**1. Introduction**

Precipitable water vapor (PWV), typically expressed in terms of height in millimeters, is the amount of water found in a vertical column of unit cross-sectional area, extending all the way from the Earth's surface to the upper edge of the atmosphere, that is potentially available for precipitation [1–3]. Although PWV accounts for less than 5% of the volume of the whole atmosphere, it is the most variable (by more than three orders of magnitude) of the major constituents of the atmosphere. In addition, PWV contributes more than any other component of the atmosphere to the greenhouse effect, playing a crucial role in the global water and energy cycles [4–7]. Therefore, the acquisition of PWV content in the air and the study of its variation over time and space are indispensable for understanding climate change and weather processes.

Traditional PWV measurement methods include radiosonde (RS), ground-based microwave radiometers, weather radar, etc. [8–13]. However, because of their high cost, most of these methods have not been widely rolled out and used on a global scale except for RS. Radiosondes are mounted on sounding balloons to record data from different pressure layers and transmit them to the ground by radio. As the most classic means of water vapor

**Citation:** Ding, J.; Chen, J.; Tang, W.; Song, Z. Spatial–Temporal Variability of Global GNSS-Derived Precipitable Water Vapor (1994–2020) and Climate Implications. *Remote Sens.* **2022**, *14*, 3493. https://doi.org/10.3390/ rs14143493

Academic Editor: Chung-yen Kuo

Received: 22 June 2022 Accepted: 19 July 2022 Published: 21 July 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

acquisition, RS has the longest record and has been considered to have the widest and most comprehensive global coverage [14–16], so RS PWV data are among the most suitable data for global analysis. In addition, because of the accuracy and reliability of RS data, they are often used as a standard to verify the observations of other instruments [17]. However, because of the high maintenance cost, the spatial and temporal resolution of sounding instrument observations is still low. There are only more than a thousand sounding stations worldwide, and usually, only two observations are made per day, mostly at 0:00 and 12:00 UTC [18,19]. Furthermore, there is a significant lack of RS data in less developed regions and in regions where extreme weather conditions exist year-round. These factors limit the application of RS for high-spatial- and -temporal-resolution water vapor monitoring and analysis. Note that numerical weather models (NWMs) are also among main sources of PWV data, and most existing PWV studies have also been based on the evaluation and analysis of NWMs. However, NWMs do not represent direct measurement, as RS and GNSS do, and may be less reliable for areas where no or limited data assimilation observations are available, so they are not discussed in this research.

As early as 1992, Bevis et al. [4] proposed the concept of global positioning system (GPS) meteorology and experimentally verified the possibility of inversion of water vapor by GPS. Ground-based GPS can operate in all weather, is virtually impervious to weather extremes, and can operate consistently over long periods of time at very low maintenance costs, which gives the PWV data obtained by this method high degrees of continuity and integrity. With the development and improvement of global and regional satellite navigation systems such as GLONASS, the BeiDou Navigation Satellite System (BDS), and Galileo, an increasing number of scholars in related fields have conducted experiments to evaluate the reliability of this technology and to verify the accuracy of GNSS (global navigation satellite system) PWV. Rocken et al. [20] assessed PWV measured with six GPS receivers for 1 month at sites in Colorado, Kansas, and Oklahoma and concluded that GPS and water vapor radiometer (WVR) estimates agreed to 1–2 mm root mean square error (RMSE). Using one and a half years of data from the western Mediterranean region, Hasse et al. [21] obtained a standard deviation (STD) difference of 1.2 cm between the RS and GNSS zenith total delay (ZTD), which was equivalent to a PWV difference of 2 mm [22]. Zhang et al. [23] obtained an STD of 3.27 mm between GPS PWV and RS PWV from 2011–2013 in Tahiti, French Polynesia, located in the tropical southern–central Pacific Ocean. The growing number of GNSS PWV studies has opened up the possibility of using GNSS to analyze high-spatial- and -temporal-resolution PWV. However, limited by experimental materials, these studies have often been conducted only regionally at only a few to a few dozen stations [24–31], and the results obtained in this way may not be applicable on a global scale. These results have validated GNSS as a good technique for acquiring PWV but have not reflected the strengths and potential of GNSS in globalscale PWV monitoring and analysis. Fortunately, the publication of the Nevada Geodetic Laboratory (NGL) troposphere products (which is probably the most comprehensive GNSS database currently available) solved this dilemma.

In this study, we examined the RS dataset from Integrated Global Radiosonde Archive Version 2 (IGRA2) and the GNSS products from the NGL. Then, the methods of PWV retrieval from these two observation techniques were introduced. To evaluate the reliability of the NGL products, the IGRA2 PWV for the 1994–2020 period was used as a reference, and the mean deviation (bias) and RMSE of the NGL PWV were calculated. Finally, the PWV data derived from more than 14,000 NGL GNSS stations over 20 years were analyzed in terms of temporal and spatial variation characteristics, anomalies in PWV were analyzed for correlations with global temperature and sea level data, and some valuable conclusions were drawn that are discussed and summarized herein.

#### **2. PWV Data and Analysis Methods**

#### *2.1. PWV Retrieval*

As PWV is not a raw observation but needs to be derived by certain methods and models, which are not exactly the same for different databases, it is necessary to give a brief introduction to the method of PWV retrieval.

#### 2.1.1. Radiosonde PWV Retrieval

The RS provides pressure level data from the ground to the top of the atmosphere, including temperature, pressure, relative humidity, geopotential height, dewpoint depression, wind direction and wind speed. RS-derived PWV is described as follows [32,33]:

$$\varepsilon = 6.112 \exp(\frac{17.6 T\_d}{T\_d + 243.15}) \tag{1}$$

$$PWV = \frac{1}{\mathcal{g}} \int\_{p1}^{p2} \frac{0.622e}{p-e} dp \tag{2}$$

where *Td* is the dewpoint temperature (◦C); *e* is the vapor pressure (hPa); *p* is the atmospheric pressure (hPa); *p*<sup>1</sup> and *p*<sup>2</sup> are the atmospheric pressures at the lower and upper layers (hPa), respectively; and *g* is the acceleration of gravity (m/s2).

#### 2.1.2. GNSS PWV Retrieval

ZTD can be divided into two components, known as zenith hydrostatic delay (ZHD) and zenith wet delay (ZWD). PWV can be derived from ZWD by the following formula [34]:

$$PWV = \Pi \times ZWD \tag{3}$$

$$\Pi = \frac{10^6}{\rho R\_v [(k\_3/T\_m) + k\_2']} \tag{4}$$

where *ρ* is the density of liquid water, *Rv* is the gas constant of water vapor, *Tm* denotes the weighted mean temperature of the atmosphere, and *k* <sup>2</sup> and *k*<sup>3</sup> are constants determined experimentally. NGL products obtained *Tm* from VMF1 (Vienna Mapping Function 1) gridded NWM data [35].

It should be noted that since the vast majority of GNSS stations are not equipped with meteorological sensors, ZHD is usually obtained from a model in GNSS data processing. However, this may result in model errors being absorbed into the ZWD and then into the PWV. The NGL solution used the VMF1 grid products as the ZHD input. The VMF1 was computed from the NWM by ray-tracing methods. When obtaining the site ZHD from the VMF1 grid products, the grid point ZHD was first converted to the corresponding height pressure, the pressure was lifted from grid height to site height, and finally, the lifted pressure was converted to ZHD again [36].

Another way to obtain high-accuracy ZHD is with reanalysis data. Reanalysis data give the best estimates where no data are recorded, in between stations and observation times. It is an "optimal" reflection of the atmospheric conditions, but not a true reflection. According to experiments by Fernandes et al. [37] and Urquhart et al. [38], the bias of the ZHD obtained from the VMF1 ZHD and the in-situ measurement was approximately 3 mm, which translated into a PWV of less than 0.5 mm. Furthermore, according to the experiments of Urquhart et al. [38] and Wang et al. [39], the bias of the ZHD difference between VMF1 and reanalysis data into PWV was around 0.3–0.4 mm. The results of these studies showed that there was no significant difference between using VMF1 and reanalysis data to obtain ZHD in regard to the accuracy of the derived PWV.

#### *2.2. RS and GNSS PWV*

In this section, the RS and GNSS data used and the data preprocessing methods used in this study are introduced.

#### 2.2.1. IGRA2 RS PWV

The Integrated Global Radiosonde Archive (IGRA) is an RS dataset maintained, archived, and distributed by the National Oceanic and Atmospheric Administration (NOAA) National Centers for Environmental Information (NCEI), formerly National Climatic Data Center (NCDC) [40]. IGRA2 is the second version of IGRA and consists of quality-controlled RS observations at stations across all continents. Its data are drawn from more than 30 different sources. The earliest year of data is 1905, and the data are updated on a daily basis. In this study, the IGRA2 RS PWV data from 1994 to 2020 were selected, covering the same period as the GNSS PWV data. To ensure the collocation of the sounding stations used for reference, the following two constraints were applied: (i) the distance between the two types of stations was less than 40 km, and the height difference was less than 100 m (existing research on this restriction has proven it to be sufficiently stringent, feasible, and effective [12]); (ii) the common available data covered at least one year. After the selection, a total of 1534 IGRA2 stations were available.

#### 2.2.2. NGL GNSS PWV

At the end of 2017, the Nevada Geodetic Lab announced the new public availability of over 34 million station-days of their tropospheric products (total zenith delay, north gradient, and east gradient every 5 min since 1996) from over 16,000 stations. The data processing strategy was precise point positioning (PPP). In March 2020, the products were updated with an earlier starting year of 1994, more than 19,000 stations, over 43 million station-days, and the addition of ZWD, *Tm*, and PWV data [35,41]. Table 1 summarizes the properties of NGL products, including the adopted models, strategies, and products used. In addition, more detailed information can be obtained from NGL's solution strategy file: http://geodesy.unr.edu/gps/ngl.acn.txt (accessed on 8 August 2021). The NGL products provided the corresponding STD of the elements listed in the table, with the exception of *Tm*, as this was provided by the VMF1 gridded numerical weather model. In addition, although the database is updated weekly, the data delay is approximately three weeks. Although it has been open for only a short time, NGL's GNSS products surpass common databases, such as the International GNSS Service (IGS) and Crustal Movement Observation Network of China (CMONOC), in terms of data volume and timespan. In particular, for the number of stations, NGL surpasses the first two by two orders of magnitude. In addition, the coverage and density of GNSS are wider and denser than those of RS in a large number of less developed and bipolar regions, which makes the NGL tropospheric products possibly the most comprehensive and extensive PWV database currently available.


**Table 1.** Adopted models and strategies of Nevada Geodetic Laboratory products.

The number of valid NGL stations, the number of total stations, and the total number of products in the database were counted for each day from 1994 to 2020 and are displayed in Figure 1. Figure 1a shows the variation in these quantities with time, and Figure 1b shows the effective time for each station (blue is the background color, each yellow horizontal line represents a station, and stations are in alphabetical order from top to bottom). The following information is relevant to Figure 1: (i) on 4 December 2004, the number of valid stations dropped abruptly from approximately 2000 to 102 and recovered immediately afterwards, probably because of a network outage; (ii) the sudden increase in the total number of stations on 1 January 2009 was due to the addition of a large number of stations from Japan to the database, with numbers ranging from J001 to J999. In addition, starting in 2009, the number of stations entered a period of rapid growth, reaching approximately 1000 new stations added per year; (iii) stations J001–J999, located in Japan, and P001–P821, on the west coast of the United States, had good observation quality, with fewer missing data and interruptions than other stations.

**Figure 1.** Variability in the numbers of Nevada Geodetic Laboratory (NGL) stations and products. (**a**) The number of valid NGL stations (red solid line), the total number of stations (red dotted line), and the total number of products in the database data (blue solid line); (**b**) the effective time for each station (each yellow horizontal line represents a station, and station names are in alphabetical order from top to bottom).

Information on the distribution of the stations is shown in Figure 2. The colors of the dots indicate the ellipsoid heights of the corresponding stations. Per the map, the majority of stations were below 1 km in ellipsoid height, and the stations with higher ellipsoid heights were mainly concentrated on the west coast of the American continent, which is where the Cordillera de los Andes is located. In addition, there were a number of stations in the Antarctic and Greenland, which are valuable for studying PWV changes in both polar regions. It should be noted that because of the large number of stations, there are inevitably some overlaps in the display. For better visualization, important information is displayed first as needed.

#### *2.3. Product Validation and Analysis Methods*

In this section, a comparison between NGL GNSS PWV and IGRA2 RS PWV is presented, and the extraction of seasonal and interannual variation signals of PWV is introduced.

**Figure 2.** Distribution of global NGL stations. Colors represent the ellipsoids of the corresponding sites. The blue pentagrams represent the stations used for power spectral density (PSD) analysis.

#### 2.3.1. Comparison with RS PWV

A total of 1534 pairs of NGL and IGRA2 stations were screened, each pair within 40 km on the ellipsoid surface and within 100 m in height, with a common data validity of more than one year (pairs of stations that did not have a full year of common observations because of missing data were excluded). The bias values, relative bias, RMSE values, and normalized RMSE values of these stations are presented in Figure 3. The biases of most of the stations were positive, and stations with negative bias values (GNSS PWV less than RS PWV) are mainly concentrated on the west coast of America and inland on the Eurasian continent. The RMSE value showed a negative correlation that decreased with increasing latitude, and for relative bias and normalized RMSE, many large values appeared in the polar regions, which indicated that the difference between the two PWV acquisition methods was larger in the high-latitude regions with low PWV content than in other regions.

**Figure 3.** Bias (**a**), relative bias (**c**) (in % relative to mean PWV), root mean square error (RMSE) (**b**), and normalized RMSE (**d**) of 1534 pairs of NGL global navigation satellite system (GNSS) and Integrated Global Radiosonde Archive Version 2 (IGRA2) RS stations.

The statistical results of the data are shown in Figure 4. By fitting a normal distribution to the frequency histogram shown in Figure 4, we revealed that for NGL GNSS PWV, the global average bias was ~0.72 ± 1.29 mm, the mean RMSE is ~2.56 ± 1.13 mm, the global average relative bias was ~0.03%, and the normalized mean RMSE was ~0.14%. Based on these results, we screened NGL stations and eliminated stations with RMSEs lower than 0.14%. It should be noted that the comparison here was only to analyze how the NGL products differed from the mainstream PWV products (there have been a large number of studies based on the RS products) and not to consider whether the RS products had higher accuracy. In fact, although NCEI scientists applied a comprehensive set of quality control procedures, IGRA data still suffered from jumps due to changes in instrumentation, etc., and the unhomogenized RS data even affected most of the reanalysis data [12]. In the GNSS solution, however, the tropospheric delay and the station coordinates were estimated together, and the errors caused by these issues were mainly absorbed by the station coordinates and hardly affected the tropospheric delay. In addition, since GNSS satellites are more than 20,000 km from the ground, they cover the full length of the water vapor path and do not need to be imputed to the top atmosphere in the same way as RS data.

**Figure 4.** Frequency histograms and normal fittings of bias (**a**), RMSE (**b**), relative bias (**c**), and normalized RMSE (**d**) of 1534 stations during 1994–2020.

2.3.2. Acquisition of PWV Spatial–Temporal Distribution Features

Significant annual periodic terms can be found in the PWV time series. To obtain all the periodic terms in the PWV series, 26 globally distributed stations (the blue pentagrams in Figure 2) were selected, which were well distributed and from which the available data were more than 20 years old. The power spectral densities (PSDs) of the PWV time series from these stations were obtained and are plotted in Figure 5. In Figure 5, the gray line is the PSD of these 26 stations, and the red line is the average PSD of the 26 stations. This figure shows that in addition to the annual cycle term, there was a significant semiannual cycle term in the PWV series. Theoretically, high-temporal-resolution GNSS data would allow smaller-period features to be extracted. However, some scholars have claimed that smaller-period temporal characteristics, such as diurnal or semidiurnal cycle characteristics, can be extracted from the reanalysis data [44]. Such features were not obtained from GNSS PWV, probably because these period amplitudes were very small (less than PWV accuracy) and were masked by observation noise.

**Figure 5.** Time series of the power spectral density (PSD) of PWV. Each gray line represents the PSD of a station, and the red line represents the average PSD value.

Trend fitting and seasonal cycle fittingwere performed usingEquations (5) and (6), respectively:

$$PVV(t) = a + b \cdot t \tag{5}$$

where *a* is the intercept, *b* is the growth rate of the PWV, and *t* is the decimal year;

$$PVV(day) = A\_1 \cos(2\pi \frac{day + P\_1}{365.25}) + A\_2 \cos(4\pi \frac{day + P\_2}{365.25}) + \mathcal{C} \tag{6}$$

where *A*<sup>1</sup> and *P*<sup>1</sup> are the annual amplitude and phase of the PWV, respectively; *A*<sup>2</sup> and *P*<sup>2</sup> are the semiannual amplitude and phase, respectively; *C* is a constant term that represents the annual average value of PWV; and *doy* is the day of the year.

#### **3. Results**

#### *3.1. Spatial Analysis*

PWV data were extracted from the ZWD while retaining most of the spatial characteristics of the ZWD. We counted all PWV data according to the latitudes and ellipsoidal heights of the stations from which they were collected. Since the temporal resolution of GNSS PWV was 5 min, this led to an overwhelming amount of data that was hard to process, so we calculated first the daily mean value and STD of the PWV and then the average value of each station's means and STDs for the 1994–2020 period. Subsequently, the PWV means and STDs were calculated by group according to latitude and ellipsoid height, and the results of the means and STDs for each group are shown in Figure 6.

**Figure 6.** Distribution of mean PWV and PWV standard deviation (STD) with latitude and ellipsoid height. (**a**) PWV mean with latitude, (**b**) PWV STD with latitude, (**c**) PWV mean with ellipsoid height, (**d**) PWV STD with ellipsoid height. Light gray dots indicate each station, red dots indicate group means, and dark gray bars represent group STD. Note: the groups for latitude and ellipsoid height are in intervals of 1◦ and 100 m, respectively.

Figure 6a shows that: (i) PWV mean values were symmetrically distributed in the Northern and Southern Hemispheres and negatively correlated with latitude, and their relationship with latitude could be fitted by a segmented linear function in which the absolute slope value for 0–35 degrees latitude was greater than that for 35–90 degrees; (ii) STD in the latitude groups decreased with increasing latitude, which indicated that the fluctuation of PWV mean values in different regions at the same latitude increased with decreasing latitude. Figure 6b shows that: (i) the mean value of the PWV daily STD decreased with increasing latitude in the regions with latitudes greater than 35◦, but there was no significant correlation in the regions with latitudes less than 35◦ which indicated that the PWV fluctuations reached the limit in low- and mid-latitude regions; (ii) in the group where the maximum value was, i.e., at approximately 35◦, the degree of variation in fluctuations was the greatest, which also indicated that the variation in PWV was the greatest among the geographical areas around this group.

From Figure 6c,d: (i) the mean PWV and STD were exponentially negatively correlated with the ellipsoid height of the stations, the PWV values varied greatly between stations in the low ellipsoid height region, and the difference decreased as the height increased; (ii) the rate of decrease in PWV with increasing ellipsoid height decreased with increasing ellipsoid height, because water vapor is mainly concentrated on the Earth's surface.

### *3.2. Seasonal Cycle*

The seasonal cycle was the most intuitive and significant temporal feature exhibited by the PWV time series. We fit the PWV using the annual periodicity + semiannual periodicity function of (11), and the fitting residuals and fitting parameters are shown in Figure 7 (Figures A1–A4 in Appendix A show the PWV fit of four stations at different latitudes). It should be noted that, considering that accurate periodicity can be extracted only from data with a sufficient timespan, we excluded stations with less than two years of data, and the fitting results for the remaining >14,000 stations are shown in Figure 7.

**Figure 7.** Global distribution of PWV fitting parameters and residuals. (**a**) Annual amplitude *A*1; (**b**) annual phase *P*1; (**c**) semiannual amplitude *A*2; (**d**) semiannual phase *P*2; (**e**) constant term *C*; (**f**) RMSE of the fitting residuals.

Figures 7a and 7b show the annual amplitude A1 and phase P1, respectively, and Figures 7c and 7d display the semiannual amplitude A2 and phase P2, respectively. The annual amplitude was between 0 and 22 mm and showed an overall pattern of increasing from high latitudes to low latitudes. However, the maximum value appeared not in the equatorial region but near the Tropic of Cancer, especially in the coastal areas of the region. This was because these regions have sufficient water vapor and a large temperature difference during the year, while the equatorial region has a small temperature difference throughout the year, and the fluctuation of water vapor content is relatively small.

The annual phase can reflect the time of the year when the PWV content is highest and the time when it is lowest. As shown in Figure 7b, the phase sign was opposite in the Northern and Southern Hemispheres. The phase value was concentrated in the periods between late May and early June in the Northern Hemisphere and in the periods between late November and early December in the Southern Hemisphere, and the semiannual cycle phase was similar to the annual phase. To obtain the distribution characteristics of the phases, statistical plots of the frequency distribution of the phase values of all the stations are shown in Figure 8. Figure 8 shows that both the annual and semiannual phases showed the phenomenon of clustering towards different values in the Northern and Southern Hemispheres. The mean value of the annual phase in the Northern Hemisphere stations was ~152.1, and the mean value of the semiannual phase was ~61.6, while the mean value of the annual phase in the Southern Hemisphere was ~−34, and the mean value of the semiannual phase was ~−40. The conversion led to the conclusion that PWV levels in the Northern Hemisphere peaked around August 1 (the 213th day of the year) and reached a trough around January 31 each year, while PWV levels in the Southern Hemisphere peaked around February 3 (the 34th day of the year) and reached a trough around August 5 each year. This indicated that the water vapor content in the atmosphere peaked approximately 40 days after the solar radiation reached its maximum (on the summer solstice) during the year. These results may be more informative for global water and energy cycle studies than some studies of global temperature through surface sensor data [45–47].

**Figure 8.** Statistical histogram and normal fitting of the annual phase (blue) and semiannual phase (red) in the Northern and Southern Hemispheres. The phase in the Southern Hemisphere was less than zero, the phase in the Northern Hemisphere was greater than zero, and the subscripts N and S represent the Northern and Southern Hemispheres, respectively.

The semiannual amplitude ranged from 0 to 8 mm, which was much smaller than the annual amplitude, and the decreasing pattern with increasing latitude was not as significant as that of the annual amplitude. The semiannual amplitude was significantly larger on the southeast coast of China–Japan and in the California axis region of the west coast of the United States than in the other regions.

The constant term C (Figure 7e) represents the annual average value of PWV, from which it was found that the annual average value of PWV ranged from 0 to 50 mm and that the values decreased with increasing latitude, from more than 40 mm in the tropical regions to less than 10 mm in the polar regions. Figure 7f shows the RMSE value of the fitted residuals, from which it was found that the polar regions had small RMSE values, larger values of RMSE appeared on the west coasts of the Pacific and Atlantic Oceans, and the global average RMSE was ~5.72 ± 1.89 mm. This indicated that these regions had greater daily variation around the seasonal means. These were also the areas where severe weather, such as typhoons and rainstorms, often occurs.

#### *3.3. Trend Analysis*

The trend signal reflects changes in the time series over long periods of time, and we used (5) to obtain a linear trend of PWV time series for over 5000 NGL stations between 2009 and 2020. All of the selected stations had more than ten years of available data (stations with slight data missing were complemented by the fitting value of (11), while stations with serious missing data were eliminated). After calculating the interannual variation in the annual mean value of each station, we found that the PWV trend had no regional distribution characteristics. Considering the large variation in PWV values between seasons, we grouped the statistics according to four seasons: spring (March– April–May for the Northern Hemisphere (NH) and September–October–November for the Southern Hemisphere (SH)), summer (June–July–August for the NH and December– January–February for the SH), autumn (September–October–November for the NH and March–April–May for the SH), and winter (December–January–February for the NH and June–July–August for the SH). The results are presented in Figure 9. In addition, stations with positive values in all four seasons and stations with negative values in all four seasons are shown in Figures 9e and 9f, respectively.

**Figure 9.** Spatial pattern of seasonal PWV trends during 1994–2020. (**a**) March–April–May (MAM), (**b**) June–July–August (JJA), (**c**) September–October–November (SON), (**d**) December–January– February (DJF), (**e**) stations with positive values in all four seasons, (**f**) stations with negative values in all four seasons.

From Figure 9a,d, we found that more stations had positive than negative values in all four seasons. However, in regard to the stations with positive (e) and negative (f) values in each of the four seasons, we found that the number of stations with positive PWV growth rates was much larger than the number with negative growth rates. In addition, the following patterns were found: (i) the PWV growth rates in Antarctica, the tropical region of South America, and Alaska were almost all positive, indicating that the PWV values in these regions were increasing year by year; (ii) a large number of stations in the Canada and South Africa regions had negative PWV growth rates for all four seasons, while very few stations had positive values, indicating that the PWV values in these regions were decreasing year by year.

Figure 10 shows the frequency distribution of the PWV growth rate of all stations in the four seasons. From the figure, it was found that all four seasons conformed to the

normal distribution, and by fitting the normal distribution, we obtained mathematical expectations of the PWV growth rate of ~0.68 ± 0.92 mm/decade, ~0.43 ± 1.14 mm/decade, ~0.80 ± 1.36 mm/decade, and ~0.64 ± 1.02 mm/decade in spring, summer, autumn, and winter, respectively. This result was numerically consistent with studies using data obtained from other data sources, including reanalysis data [48], RS data [49], MODIS (moderate resolution imaging spectroradiometer) data [50], etc. This indicated that: (i) the global PWV content showed a yearly increasing trend in all four seasons, which was most significant in autumn and least significant in summer, and (ii) the dispersion (standard deviation) of the PWV growth rate among stations increased gradually in spring, summer, and autumn but was lowest in winter.

**Figure 10.** Seasonal frequency histograms of PWV trends during 1994–2020 for (**a**) MAM, (**b**) JJA, (**c**) SON, and (**d**) DJF.

In addition, it was found from the results that the PWV growth rate had large uncertainty and that its spatial distribution was not continuous in some regions. This phenomenon was also found in some studies of long-term trends in PWV obtained using other methods [50,51]. Combining these studies, the following possible reasons for this phenomenon are listed: (i) PWV is highly correlated with land cover type [52], and the density of discrete GNSS stations is still low; (ii) differences in analysis periods and spatial coverage could be responsible for these discrepancies [49].

Table 2 shows the percentages of stations with positive and negative PWV growth rates in each of the four seasons to the total number of stations. From the table, it was found that in each season, the proportion of positive values was approximately triple to quadruple that of negative values, and the proportion of stations with positive values in all seasons was more than 20 times the proportion of stations with negative values in all seasons. This indicated that PWV was increasing year by year on a global scale.


**Table 2.** Percentages of positive and negative values for the four seasons and all seasons.

The relationships between the PWV growth rate and both latitude and ellipsoid height were likewise analyzed, with the latitude and ellipsoid height grouped in the same way as above, and the results are presented in Figure 11. It should be noted that the increase in PWV did not mean that the evaporation in this area had also increased, but it meant that increased water vapor was concentrated in these areas. The results in Figure 11 showed that the PWV growth rate was not simply positively or negatively correlated with latitude and ellipsoid height. It was found that the increased water vapor was more concentrated near the equator, the Tropic of Cancer, and 60◦ latitude, while from the ellipsoid height, the increased water vapor was more concentrated in the area within 1 km from the ground (these areas all had growth rates higher than 0.5 mm/decade). This suggested that these regions are more likely to suffer drastic weather changes due to uneven increases in water vapor.

**Figure 11.** Distribution of PWV growth rates with latitude and ellipsoid height. (**a**) PWV growth rate with latitude; (**b**) PWV growth rate with ellipsoid height. The error bars are STDs of PWVs of stations in the same latitude/height grouping area. Note: the groups for latitude and ellipsoid height are in intervals of 5◦ and 100 m, respectively.

Under a warming climate, increasing temperatures exacerbate the evaporation of water, a process that in turn is enhanced by increasing levels of water vapor, the world's most dominant greenhouse gas. To establish the link between the growth rates of PWV and temperature, etc. and to analyze the correlation between them, we obtained global temperature anomalies and sea height variations from Berkeley Earth [53] and NASA's Goddard Institute for Space Studies (GISS) [54] for the last two decades and compared PWV anomalies with them. After data processing, the results of the comparison are shown in Figure 12, where the red points represent the global annual values of anomalies (variations), the red error bars are the yearly fluctuations in temperature anomalies and sea height variations, and the blue error bars are the yearly fluctuation of PWV anomalies. The calculation of the PWV anomalies followed that of temperature anomalies and sea height variations, with all being calculated as the differences relative to the first value. Note that the data used here were annual averages of the global data, which were calculated from monthly averages, and the yearly fluctuation is the STD of these data. The results shown in Figure 12 showed strong positive correlations between PWV anomalies and both temperature anomalies and sea level variations, with correlation coefficients of 0.81. This

result was similar to that obtained by Kwon et al. [50] with a correlation coefficient of 0.8 using MODIS PWV compared with temperature anomalies. In addition, the results of the linear fit indicated that a 1 ◦C temperature increase corresponded to a PWV increase of ~2.075 ± 0.765 mm, and a 1 mm sea level height increase corresponded to a PWV increase of ~0.015 ± 0.005 mm.

**Figure 12.** Comparison of PWV anomalies with temperature anomalies (**a**) and sea height variations (**b**). The red dots are the annual averages, the red error bars are the yearly fluctuations in temperature anomalies and sea height variations, and the blue error bars are the yearly fluctuations in PWV anomalies. k represents the slope of the linear fit, and R represents the correlation coefficient.

#### **4. Discussion**

In this contribution, the spatial and temporal distribution characteristics of GNSSderived PWV for the last two decades from more than 14,000 stations distributed globally were extracted and analyzed. In the spatial analysis, we found that the annual STDs of PWV tended to be flat in the middle- and low-latitude regions, especially at low latitudes, which may be due to the saturation of water vapor content in this region. In the seasonal analysis, the amplitude and phase of each station was extracted, which could be used as model values for PWV recovery, as were the spatial distribution characteristics of each station. In the trend analysis, the interannual variability of each station was extracted, and correlations with the interannual variability of temperature and sea level were obtained, which can be used as additional supporting evidence for global warming.

Compared with existing studies on PWV (which have mostly been based on reanalysis data or RS data), this study did not have an advantage in terms of the timespan of the data because of the short development time of GNSS. Furthermore, it is difficult to obtain long-term GNSS PWV signatures for marine areas, since ground-based GNSS stations are usually fixed to bedrock.

Finally, although there are already many stations in the NGL database, they are still not evenly distributed spatially and are mainly concentrated in developed regions. Therefore, greater efforts need to be made to build stations, especially in northern–central Africa, the Asian continent, and the polar regions. With the development of low-cost, small GNSS receivers, the acquisition of GNSS-based PWV will become more and more convenient. The data records of GNSS PWV will become longer, and the coverage will become wider. GNSS will become one of the most mainstream ways to acquire PWV.

#### **5. Conclusions**

Reliable temporal and spatial variability information on water vapor is crucial to understanding climate processes and the global water and energy cycles. This study used PWV data derived from IGRA2 RS to evaluate NGL GNSS PWV products in 1994–2020. A total of 1534 RS sites were selected as the reference. The comparison showed a very high level of agreement between these two types of PWV data. Then, the temporal and

spatial distribution characteristics of PWV from 1994–2020 for over 14,000 NGL stations worldwide were analyzed and discussed. In addition, long-term PWV anomalies were analyzed in comparison with temperature anomalies and sea surface height variability data. The following conclusions were obtained:


**Author Contributions:** Conceptualization, J.D.; writing—original draft preparation, J.D.; writing—review and editing, J.C., J.D., W.T. and Z.S.; visualization, J.D.; validation, J.D. and J.C.; supervision, J.C.; funding acquisition, J.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Program of Shanghai Academic/Technology Research Leader (No.20XD1404500); the National Natural Science Foundation of China (No.11673050); the Key Program of Special Development funds of Zhangjiang National Innovation Demonstration Zone (Grant No. ZJ2018-ZD−009); the National Key R&D Program of China (No. 2018 YFB0504300); and the Key R&D Program of Guangdong province (No. 2018 B030325001).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The radiosonde data were provided by the National Oceanic and Atmospheric Administration's (NOAA's) NCEI (National Centers for Environmental Information) and are available from ftp://ftp.ncdc.noaa.gov/pub/data/igra/ (accessed on 8 August 2021). The GNSS PWV data were provided by the Nevada Geodetic Lab at http://geodesy.unr.edu/gps\_ timeseries/trop/ (accessed on 17 August 2021). The global temperature and sea level data were provided by Berkeley Earth and NASA's Goddard Institute for Space Studies (GISS) at http:// berkeleyearth.org/data/ (accessed on 9 November 2021) and https://climate.nasa.gov/vital-signs/ (accessed on 9 November 2021), respectively.

**Acknowledgments:** The authors gratefully acknowledge the NOAA's NCEI, the Nevada Geodetic Lab, Berkeley Earth, and NASA's Goddard Institute for Space Studies for providing the datasets.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A**

**Figure A1.** PWV time series (blue dots) and its fitting (red line) from the GNSS station TIXI.

**Figure A2.** PWV time series (blue dots) and its fitting (red line) from the GNSS station PTBB.

**Figure A3.** PWV time series (blue dots) and its fitting (red line) from the GNSS station CAGS.

**Figure A4.** PWV time series (blue dots) and its fitting (red line) from the GNSS station BRAZ.

### **References**

