Investigation of Antarctic Precipitable Water Vapor Variability and Trend from 18 Year (2001 to 2018) Data of Four Reanalyses Based on Radiosonde and GNSS Observations

Abstract

Precipitable water vapor (PWV) plays a vital role in climate research, especially for Antarctica in which meteorological observations are insufficient due to the adverse climate and topography therein. Reanalysis data sets provide a great opportunity for Antarctic water vapor research. This study investigates the climatological PWV means, variability and trends over Antarctica from four reanalyses, including the fifth generation of European Centre for Medium-Range Weather Forecasts (ECMWF) Reanalysis (ERA5), the Second Modern-Era Retrospective analysis for Research and Applications (MERRA-2), Japanese 55-year Reanalysis (JRA-55) and National Centers for Environmental Prediction/Department of Energy (NCEP/DOE), in the period of 2001–2018 based on radiosonde and GNSS observations. PWV data from the ERA5, MERRA-2, JRA-55 and NCEP/DOE have been evaluated by radiosonde and GNSS observations, showing that ERA5 and MERRA-2 perform better than JRA-55 and NCEP/DOE with mean root mean square (RMS) errors below 1.2 mm. The climatological PWV mean distribution over Antarctica roughly shows a decreasing trend from west to east, with the highest content in summer and the lowest content in winter. The PWV variability is generally small over Antarctica, showing a seasonal dependence that is larger in the cold season and smaller in the warm season. PWV trends for all reanalyses at most Antarctic regions are insignificant and most reanalyses present overall drying trends from 2001 to 2018, except for ERA5 exhibiting a moistening trend. PWV trends also show seasonal and regional dependence. All reanalyses are generally consistent with radiosonde and GNSS observations in reproducing the PWV means (mean differences within 1.1 mm), variability (mean differences within 3%) and trends (mean differences within 6.4% decade⁻¹) over Antarctica, except for NCEP/DOE showing spurious variability and trends in East Antarctica. Results can help us further understand these four reanalysis PWV products and promote climate research in Antarctica.

1. Introduction

Water vapor, as one of the major greenhouse gases in the atmosphere, accounts for about 60% of the natural greenhouse effect and plays a vital role in many atmospheric processes, such as global energy and hydrological cycle, global warming and climate change [1]. Following the nonlinear Clausius–Clapeyron equation, the water vapor is mainly controlled by temperature, a 1-K rise of temperature will increase about 7% water-holding capacity of the atmosphere if the relative humidity remains nearly constant, and the increase in water vapor will further induce a positive feedback in global warming [2,3]. This strong positive feedback has a profound influence on the global climate system and affects the development of extreme climate phenomena due to rapid variations of water vapor accompanied by large changes of the latent heat [4]. Antarctica, the continent with the highest average elevation around the world, is characterized by the coldest, driest, stormiest and windiest climate with scarce precipitation in the world. The annual mean precipitation of the Antarctic continent is 55 mm and there is almost no precipitation near the South Pole. In Antarctica, the spatiotemporal variability of atmospheric water vapor impacts the precipitation pattern and thus the rates of ice accumulation over the continental interior [5]. Therefore, the monitoring and investigating of atmospheric water vapor content and its variations are significant for a better understanding of global climate and weather processes, particularly in such a sensible region as Antarctica.

Precipitable water vapor (PWV) is one of the most widely used indicators for measuring the amount of atmospheric water vapor and an important parameter in numerical weather forecasting models [6,7,8,9]. It can be acquired through various techniques and generally can be divided into two categories. The first category is observations from multiple sensors including radiosondes, Global Navigation Satellite System (GNSS), weather stations and satellite remote sensing [10,11,12,13,14]. The PWV observations obtained from this category are measured by the sensors and have the advantage of high accuracy. Among them, ground-based GNSS can detect the signal zenith delay and further estimate the PWV with its long-term stability, weather independent and high temporal resolution compared to other PWV measurements such as radiosondes and satellite remote sensing [15,16,17]. However, since the sensors from radiosondes and satellites are of high cost and poor performance under harsh climate conditions, and since the observations are often sparsely distributed and incomplete, it is difficult to estimate accurate PWV variability and trend from the observational data. Another category of PWV data is the reanalysis data set that is produced by assimilating multiple sources of observations (such as radiosonde, ground-based weather station and satellite remote sensing) into a coherent data set based on the advanced data assimilation scheme and atmospheric numerical model [18,19]. Reanalysis data sets can provide complete and homogeneous atmospheric data over multiple decades with the advantages of multivariable outputs and global coverage. They thereby can be of great value for climate research and applications. In Antarctica, meteorological observations are sparse due to the inhospitable environment; thus, studies of Antarctic climate often rely upon the reanalysis of data sets. Although unified processing schemes and models are used in the reanalysis data, biases and inhomogeneities from the assimilated observations can still contaminate the reanalysis data sets, especially for Antarctica where observations used for data assimilation are scarce or contain errors. It is therefore possible that PWV variability and trend estimates derived from the reanalysis data sets may show spurious phenomenon, and this cannot be ignored.

Many studies have been conducted to investigate global or Antarctic PWV changes. Trenberth et al. [1] performed an analysis and evaluation of global PWV means, variability and trends for the National Centers for Environmental Prediction (NCEP) reanalyses 1 and 2, and for the 40-year European Centre for Medium-Range Weather Forecasts (ECMWF) reanalysis (ERA-40) from 1988 to 2001. The reanalyses showed reasonable results over land while major problems were found over the oceans. Suparta et al. [20,21,22] used GPS technology to analyze the Antarctic PWV variability and its response to some atmospheric events including solar activities and terrestrial winds. Chen and Liu [23] analyzed the global variability and trends in PWV using ERA-Interim, NCEP/DOE (Department of Energy), radiosonde, GPS and microwave satellite data, and found that the global PWV had an evident moistening trend while the Antarctic PWV showed a drying trend in ERA-Interim from 1992 to 2014. However, the result is debatable for Antarctica because of the shortage of observational data. Parracho et al. [24] also investigated the global means, variability and trends in PWV from ERA-Interim, the Second Modern-Era Retrospective analysis for Research and Applications (MERRA-2) and GNSS for the period 1980–2016, however large uncertainties in PWV variability and trends were still seen in both reanalyses over Antarctica. Ricaud et al. [25] investigated the time evolution of humidity and temperature above Antarctic Dome C using radiosondes, microwave radiometer and reanalyses. The investigations of global water vapor distribution and changes were performed in previous studies. However, most studies are limited in providing detailed water vapor distribution and changes for Antarctica. A comprehensive investigation of Antarctic water vapor has yet to be fully explored. Moreover, the reliability of the PWV derived from reanalyses also needs to be evaluated.

As the latest released reanalysis data set from ECMWF, ERA5 has the highest temporal resolution of 1 h and highest horizontal spatial resolution of 0.25° × 0.25° (lon. × lat.) with respect to other older reanalyses including MERRA-2, Japanese 55-year Reanalysis (JRA-55) and NCEP reanalysis. ERA5 benefits from many improvements, including assimilation of more recent observations and the adoption of more advanced data assimilation methods [26]. In this study, we first verify the ERA5, MERRA-2, JRA-55 and NCEP/DOE-derived PWV results with respect to radiosonde and GNSS observations. Then, once the reanalysis-derived PWV accuracy has been verified, a comprehensive investigation of the Antarctic PWV means, variability and trends for the period 2001–2018 are presented using reanalyses, GNSS and radiosonde data. Data and methods used in this paper are described in Section 2. Section 3 shows the evaluation of PWV estimations derived from reanalyses using radiosonde and GNSS observations. The means, variability and linear trends of PWV found in the reanalyses, radiosonde and GNSS data over Antarctica from 2001 to 2018 are carried out in Section 4. Section 5 discusses the results and concludes the paper.

2. Data and Methods

2.1. Four Reanalysis Products

Four global reanalysis products are used in this study, including ERA5, MERRA-2, JRA-55 and NCEP/DOE reanalysis. The basic information of these four reanalysis products is summarized in

Table 1. Summary of four reanalysis products used in this study.

Reanalysis	Horizontal Resolution (Lon. × Lat.)	Temporal Resolution	Assimilation System	Source
ERA5	0.25° × 0.25°	1 h	4D-VAR	ECMWF
MERRA-2	0.625° × 0.5°	1 h	3D-VAR	NASA GES DISC
JRA-55	1.25° × 1.25°	6 h	4D-VAR	JMA/CRIEPI
NCEP/DOE	2.5° × 2.5°	6 h	3D-VAR	NCEP/DOE

.

ERA5 is the fifth generation ECMWF atmospheric reanalysis data set for the global climate and weather, covering the period from 1 January 1979 to present. The ERA5 reanalysis is the successor of ERA-Interim reanalysis and assimilates as many observations as possible in the upper air and near surface with the four-dimensional variational analysis (4D-Var) data assimilation scheme in the Cycle 41R2 model of the ECMWF’s Integrated Forecast System (IFS), which is expected to be extensively used in the future [26,27,28].

MERRA-2 is the latest global atmospheric reanalysis data set produced by NASA’s Global Modeling and Assimilation Office (GMAO) from 1980 to present. It uses an upgraded version of the Goddard Earth Observing System Model, Version 5 (GEOS-5) data assimilation system [29] and the Gridpoint Statistical Interpolation (GSI) analysis scheme [30,31], which replaces the original MERRA reanalysis data set.

JRA-55 is the second Japanese global atmospheric reanalysis project conducted by the Japan Meteorological Agency (JMA) for the period from 1958 onward. Compared to its predecessor JRA-25, JRA-55 employs a new data assimilation and prediction system (DA) which improves many deficiencies found in JRA-25. These improvements include implementing higher spatial resolution (TL319L60), a new radiation scheme in the forecast model, the 4D-Var assimilation scheme with Variational Bias Correction (VarBC) and introduction of greenhouse gases with time-varying concentrations [32,33].

The NCEP/DOE Reanalysis-2 is using a state-of-the-art analysis/forecast system to perform data assimilation using past data from 1979 up to present in the Atmospheric Model Intercomparison Project (AMIP-II) to correct several human processing errors and update parameters of physical processes in its predecessor, NCEP/NCAR (National Center for Atmospheric Research) reanalysis [34].

These four reanalysis data sets have adopted advanced numerical climate models and assimilation methods to combine multi-source observations which minimized bias, and it is expected that PWV variability and trend estimates are more realistic. In this study, these four reanalysis data over Antarctica (60–90°S area) for the period 2001−2018 are used.

2.2. PWV Retrieval from Radiosonde Profiles

Radiosonde stations can directly measure vertical meteorological parameters by a sounding balloon at 00:00 and 12:00 UTC every day at various heights, containing mainly temperature, pressure and relative humidity. The radiosonde-derived PWV can be calculated from an integral of the meteorological profiles as follows [35]:

$$\mathrm{PWV}={\displaystyle {\int }_{P_i}^{P_{i+1}}\frac{q}{{\rho }_{w}g}}dP,$$

(1)

where P is the atmospheric pressure, P_i and P_i+1 are the pressure at the lower and upper layers, respectively (unit: Pa). q refers to the specific humidity (g∙kg⁻¹), ρ_w denotes the liquid water density (g∙cm⁻³) and g is the acceleration of gravity (cm∙s⁻²). The accuracy of PWV estimations derived from radiosonde has been demonstrated to be high with the uncertainty of 5~8% based on different radiosonde types which constantly serve as a standard to confirm other PWV measurement results [11,12]. In this study, the PWV estimations derived from eleven radiosonde sites are considered to assess the reliability of the reanalysis and GNSS PWV results and to derive the Antarctic water vapor variability and trends. For the quality control of radiosonde-derived PWV time series at each station, estimations with deviation from the mean value of estimations outside five times the standard deviation are rejected [36]. The geographical distribution of the eleven radiosonde sites is shown in Figure 1 (red triangles).

2.3. PWV Retrieval from GNSS

The International GNSS Service (IGS) center is operating a worldwide GNSS network with more than 400 GNSS reference stations and provides regularly a 5 min sampling ZTD (zenith total delay) tropospheric product. The accuracy of IGS ZTD estimates has been proven to be better than 5 mm [37,38,39]. To ensure the reliability of the IGS ZTD estimates, the ZTD data are screened using a quality control scheme described by Bock et al. [40], and the abnormal observation data are removed. After screening, the ZTD data from nine GNSS stations over Antarctica for the period 2001–2018 (the selection retained those with the time series more than ten years) were validated to be used in this study. The 5 min ZTD data are averaged in 1-h bins. The geographical distribution of the nine GNSS stations is shown in Figure 1 (yellow circles).

The conversion of ZWD (zenith wet delay) to PWV is through the following formula [41]:

$$\mathrm{PWV}=\mathrm{\Pi }\cdot \mathrm{ZWD}=\frac{{10}^6}{{\rho }_w{R}_v[(k_2^{\prime }+k_3/T_m)]}\cdot \mathrm{ZWD}$$

(2)

where Π is the water vapor conversion factor, a function of atmospheric weighted mean temperature (T_m); R_v denotes the specific vapor gas constant (461.495 J∙kg⁻¹∙K⁻¹); ρ_w refers to the density of liquid water (1 × 10³ kg/m³); k’₂ and k₃ are the atmospheric refractivity parameters given by Bevis et al. [42]; ZWD is obtained by subtracting the zenith hydrostatic delay (ZHD) from the ZTD and the ZHD (in m) is calculated by the Saastamoinen model [43]:

$$\mathrm{ZHD}=0.0022768\times \frac{P_s}{1-0.00266\cdot \cos (2\phi )-2.8\times {10}^{-7}\times h_o}$$

(3)

where φ represents the latitude of the GNSS station (unit: radians) and h_o is the height of the GNSS station (unit: m); P_s is the surface pressure of the GNSS station (unit: hPa). In this study, P_s at GNSS site is obtained from ERA5 data through spatial interpolation of four surrounding grid points, and the detailed interpolation process is described in Section 2.4. T_m is directly obtained from the Global Grid T_m (GGTm) model proposed by Huang et al. [44]. The GGTm model considers the latitude and altitude variations in T_m estimates and shows better performance than the widely used GPT2w model [45] in T_m estimates, especially in the regions with large undulating terrain such as Antarctic areas [44]. Therefore, the GGTm model is used to obtain the T_m estimates for each GNSS station. It should be noted that negative GNSS-derived PWV results were removed for quality assurance. Negative GNSS-derived PWV values are due to a poor performance of the underlying ZHD models (Saastamoinen model) over Antarctica [46]. In Antarctica, the PWV contents are far smaller than in other parts of the world and the GNSS-derived PWV must be considered with great care even if they are positive. In this work, the quality of the PWV retrieval from GNSS is verified using radiosonde observations. Then, the GNSS-derived PWV is adopted as a reference value to validate the PWV results estimated from the reanalysis data and to derive the PWV variability and trends.

2.4. Spatial Adjustments for Meteorological Data

Since the GNSS data are referred on the geodetic height system, both radiosonde and reanalysis data adopt the geopotential height system, and it is crucial to unify the different height systems when two types of data are compared. In this study, we employed the scheme described in Wang et al. [47] to unify the different height systems. To adjust the reanalysis data to each station, we first identified the four surrounding grid points in horizontal. Then, a vertical adjustment is implemented to adjust the grid points to the height of the station. Finally, the grid points are interpolated bi-linearly to the location of the station. For the vertical adjustment, different meteorological parameters depend on different vertical adjustment methods. The vertical adjustment formula for pressure is as follows [48]:

$$P=P_0{\left[1-0.0000226(h-h_0)\right]}^{5.225},$$

(4)

where P and P₀ are the surface pressure at GNSS station and grid point, respectively; h and h₀ represent the height of the GNSS station and grid point, respectively. For the PWV vertical adjustment, an empirical equation derived by Leckner et al. [49] and Kouba [50] is widely used to adjust for the height difference:

$${\mathrm{PWV}}_{h_1}={\mathrm{PWV}}_{h_2}\cdot \exp \left(\beta \cdot \left(h_1-h_2\right)\right),$$

(5)

where ${\mathrm{PWV}}_{h_1}$ and ${\mathrm{PWV}}_{h_2}$ represent the PWVs corresponding to the heights of h₁ and h₂ (unit: km), respectively; β denotes the empirical decay coefficient of −0.5 mm/km given by Kouba [50]. The advantage of this equation is that it can be applied directly to the PWV data without requiring any extra meteorological data. However, the β value given by Kouba is obtained from the average of the gridded—site-specific VMF1 (Vienna Mapping Function 1) differences of wet delay for IGS station KOKB [50]. Moreover, the β may have some seasonal variations due to the PWV seasonal changes. Therefore, the β value should be different in Antarctica and the parameter should be recalculated. In this work, we use the method proposed by Huang et al. [51] to develop a PWV vertical correction model for Antarctica using the PWV derived from ERA5 data.

To investigate whether the PWV and its corresponding height would follow the exponential decay equation shown in Equation (5) in Antarctica, the daily ERA5 gridded PWV data in 2018 and corresponding gridded geopotential height data are used to carry out a regression analysis to explore the relationship between PWV and geopotential height over Antarctica. After experiments with MATLAB cftools, we found that the exponential equation could fit well the relation, and the regression analysis result is shown in Figure 2.

In Figure 2, the regression analysis shows that a significant exponential decay correlation between the PWV and the corresponding height with the correlation coefficient (R) of −0.85, and the exponential decay equation, can fit well PWV versus geopotential height with RMS error of 1.1 mm, indicating that PWV is highly related to elevation in Antarctica, and their relation can be determined as Equation (5). The daily β values are obtained using daily ERA5-derived PWV and the corresponding geopotential height data over Antarctica from 2001 to 2018, and the seasonal variations of the β are investigated. In addition, the fast Fourier transform (FFT) algorithm is performed to detect the characteristic periodicity of β. The results are shown in Figure 3.

As shown in Figure 3, β clearly shows a strong correlation with time and seasonal variations, therefore, β can be modeled as:

$$\begin{array}{ll}\beta (\mathrm{DOY})=& {\beta }_0+{\beta }_1\cos (\frac{\mathrm{DOY}}{365.25}2\pi )+{\beta }_2\sin (\frac{\mathrm{DOY}}{365.25}2\pi )\\ & +{\beta }_3\cos (\frac{\mathrm{DOY}}{365.25}4\pi )+{\beta }_4\sin (\frac{\mathrm{DOY}}{365.25}4\pi )\end{array}$$

(6)

where DOY is the day of the year; β₀ is the mean value of the β; (β₁, β₂) and (β₃, β₄) are the annual and semiannual periodicity coefficients of the β, respectively. These model coefficients can be estimated by applying the least squares method to the daily β from 2001 to 2018 based on Equation (6), and the model coefficient values of β₀, β₁, β₂, β₃ and β₄ are −0.799, 0.058, 0.003, 0.035 and −0.003, respectively. When using this regional PWV vertical correction model for Antarctica, only the PWV, height, target height and DOY are required. The PWV values at the target height can be calculated using the model through Equations (5) and (6).

2.5. Method for Linear Trend Estimation

PWV linear trends in time were estimated using the Theil–Sen method [52]. The Theil–Sen method is a non-parametric method that estimates the trend by calculating the median of slopes of lines through all pairs of points in the time series. This method provides a more robust slope estimate than the least squares method since it is insensitive to extreme or outlier values, particularly for relatively small data sets, and it can also reduce the effects of autocorrelation in time series [53]. The Theil–Sen method is exploited to estimate the PWV trends are as:

$$\mathrm{trend}=\underset{1\le i<j\le n}{\mathrm{median}}\left(\frac{y_j-y_i}{t_j-t_i}\right),$$

(7)

where y_j and y_i are the monthly PWV anomaly values at time t_j and t_i (j > i), respectively. The monthly PWV anomalies are obtained by removing the monthly PWV climatology (mean value of monthly PWV over the whole period from 2001 to 2018) from the monthly PWV. To understand the PWV change intuitively, the percentage trend P_trend (unit: % decade⁻¹) is used in this work [54]:

$$P_{\text{trend}}=1000\cdot \frac{\mathrm{trend}}{\overline{y}},$$

(8)

where $\overline{y}$ is the annual average PWV value. The statistical significance of trends is assessed using a modified Mann–Kendall trend test [55] with a 10% significance level.

In the analysis of the PWV time series, the daily PWV values are obtained by averaging hourly PWV data, the monthly PWV values are generated by averaging daily PWV values and then the annual PWV values are calculated by averaging the monthly PWV values. In addition, the seasonal PWV values are computed for December-January-February (DJF), March-April-May (MAM), June-July-August (JJA) and September-October-November (SON) by averaging the monthly PWV values in each season. To ensure that the computed PWV values can represent the monthly, seasonal and annual means, the monthly PWV is computed from the daily PWV when at least half of the expected monthly data (i.e., at least nearly 15 days of data) are available, the seasonal PWV is calculated from the monthly PWV if at least 2 months of data are available and the annual PWV is generated from the monthly PWV when at least 10 out of 12 months of data are available.

3. Assessment of PWV from Reanalyses by Radiosonde and GNSS Observations

The reliability of the reanalyses is required to be assessed before they are further applied to estimate PWV variability and trends. As mentioned above, radiosonde and GNSS techniques can measure the PWV content with high accuracy. Therefore, we carried out an evaluation to evaluate the ERA5, MERRA-2, JRA-55 and NCEP/DOE-derived PWV using radiosonde and GNSS observations. First, the PWV derived from GNSS observations were evaluated using collocated radiosonde stations. After verifying the GNSS PWV results with respect to radiosonde, the four reanalysis PWV products were evaluated using 18 years of concurrent estimates at radiosonde and GNSS stations. A quality control method [36] is implemented in each PWV comparison: PWV differences with deviation from the mean value larger than three times the standard deviation are rejected as gross errors.

3.1. Comparison of GNSS-Derived PWV with Radiosonde Observations

As independent measurements, radiosonde observations are constantly used as a standard to confirm the results of the GNSS-derived PWV [11]. The quality of the GNSS-derived PWV needs to be verified using radiosonde observations before using GNSS to assess the reanalysis data. To minimize the impacts caused by distance and height differences between the GNSS and radiosonde stations, five collocated sites were selected for the PWV comparison, and the distances and height differences between the GNSS and collocated radiosonde stations are all within 2 km and 100 m, respectively. Comparisons of GNSS-derived PWV with radiosonde-derived PWV at the five collocated sites during the period from 2001 to 2018 are shown in Figure 4, and statistical comparison results are shown in

Table 2. Statistical results of the comparison between GNSS-derived PWV and radiosonde-derived PWV at the five collocated sites from 2001 to 2018.

Site	Bias (mm)	RMS (mm)	R
CAS1	1.18	1.35	0.945
DAV1	0.72	0.94	0.945
MAW1	1.54	1.68	0.923
MCM4	0.73	0.97	0.93
SYOG	0.73	1	0.933

.

According to Figure 4, the GNSS-derived PWV at all five stations show good agreement with the corresponding radiosonde estimates. Compared with the radiosonde-derived PWV, the GNSS PWV shows a strong correlation with the radiosonde PWV at the five collocated stations with a mean R value of 0.935, and all the R values are larger than 0.92. In Table 2, the mean biases between the radiosonde-derived and GNSS-derived PWV at the five collocated stations are all positive with the range of 0.72–1.54 mm, indicating that the GNSS-derived PWV is prone to slight overestimation in Antarctica. Unfortunately, the reason for the bias remains to be determined. Potential factors contributing to the bias may be the antenna phase center variation at GNSS stations, which directly impacts the GNSS PWV estimates [56], or uncertainty in radiosonde observations [57] or the position offsets between the radiosonde and GNSS stations [58]. In terms of RMS error, the overall mean RMS error at the five sites is 1.12 mm with the range of 0.94–1.68 mm, which shows high accuracy and is consistent with the findings of previous studies [12,15]. Therefore, the GNSS-derived PWV could be used as a reference to assess the PWV estimates derived from reanalyses.

3.2. Verification of Reanalysis-Derived PWV by Radiosonde and GNSS Observations

The accuracies of the four reanalysis-derived PWV are evaluated using concurrent PWV estimates from radiosonde and GNSS observations.

Table 3. Statistical results of the assessment of PWV derived from four reanalyses by radiosonde- and GNSS-derived PWV from 2001 to 2018.

Comparison		Bias (mm)	RMS (mm)	Relative RMS	R
Radiosonde versus	ERA5	−0.1	0.5	15.45%	0.97
	MERRA-2	0.14	0.65	20.56%	0.95
	JRA-55	0.09	0.85	26.54%	0.91
	NCEP/DOE	0.69	1.11	34.88%	0.95
GNSS versus	ERA5	−0.88	1.2	26%	0.96
	MERRA-2	−0.52	1.15	25.08%	0.94
	JRA-55	−0.66	1.34	28.83%	0.92
	NCEP/DOE	0.09	1.59	35.31%	0.88

gives the bias, RMS, relative RMS and R of the differences between ERA5/MERRA-2/JRA-55/NCEP-derived PWVs and radiosonde/GNSS-measured PWVs during the period from 2001 to 2018. The spatial distributions of the bias, RMS and R of the four reanalysis products validated by radiosonde and GNSS stations are depicted in Figure 5 and Figure 6, respectively.

In the evaluation with respect to radiosonde, ERA5 gets the highest accuracy with respect to MERRA-2, JRA-55 and NCEP/DOE with mean bias, RMS, relative RMS and R values of −0.1 mm, 0.5 mm, 15.45% and 0.97, respectively. However, NCEP/DOE shows a worse accuracy than the other three reanalyses with mean bias, RMS, relative RMS and R values of 0.69 mm, 1.11 mm, 34.88% and 0.95, respectively. The R values between the reanalysis-derived PWV and radiosonde-derived PWV for eleven stations vary between 1 and 0.8, which indicates that the PWVs derived from reanalysis and radiosonde have strong correlations at these sites. As shown in Figure 5, compared with radiosonde, the R values of ERA5 and MERRA-2 are slightly higher than JRA-55 and NCEP/DOE among the radiosonde stations (see Figure 5i–l). Positive biases are observed in most sites except for ERA5, which are nearly 0 mm level. The RMS errors of the four reanalyses are all below 2 mm at all 11 radiosonde stations, and ERA5 shows a best performance with respect to other reanalyses, while NCEP/DOE performs the worst. All reanalyses seem to have low biases and RMS errors; however, NCEP/DOE has a larger relative RMS error of 34.88%. It is better to refer to relative RMS errors instead of considering biases and RMS errors only, especially in Antarctica where the PWV contents are very low. The sites with larger RMS errors in the four reanalyses are mainly located in coastal East Antarctica; the water vapor and its uncertainty are relatively large under the influence of the ocean.

In the evaluation with respect to GNSS, the PWV accuracies of ERA5 and MERRA-2 are higher than JRA-55 and NCEP/DOE among the four reanalyses in accordance with the mean RMS and relative RMS values (see Table 3). Nevertheless, the biases, RMS and relative RMS errors for all four reanalyses are larger than those with respect to radiosonde. The main reason is attributed to the fact that radiosonde observations are already assimilated to the reanalyses [26,51]. Therefore, it is normal to find a better agreement as radiosonde data were assimilated in reanalysis models. GNSS-observed PWV are independent of reanalysis models and constitute a robust validation data set. Compared with the GNSS estimates, the PWVs derived from four reanalyses show strong correlations with the GNSS-derived PWV at all nine sites with the R values ranging from 0.76 to 1. In Figure 6, all reanalyses show significant negative biases at most stations except for NCEP/DOE, indicating that the reanalysis-derived PWV is slightly lower than the GNSS-derived PWV in Antarctica. The RMS errors of the four reanalyses are all below 2.2 mm at all nine GNSS stations and the relative RMS errors are all below 36%. At most GNSS sites, NCEP/DOE has the poorest performance among the four reanalyses with higher RMS values, while the other three reanalyses show lower RMS values.

For an each-year comparison, Figure 7 presents the variations of bias and RMS error with year. In the comparison between radiosonde and reanalysis data, the biases and RMS errors of four reanalyses are small and steady without large fluctuations, except for JRA-55 which has a sudden increase in RMS errors from 2017 to 2018. The RMS errors of JRA-55 in 2017–2018 are nearly 1 mm larger than those of the other years, which are caused by the large differences between JRA-55 and radiosonde observations at several radiosonde stations. This issue in JRA-55 requires to be further investigated once more in situ observations are available in Antarctica. Compared with radiosonde, the biases and RMS errors have larger fluctuations in the comparison between GNSS and reanalysis data. Such issues may be caused by the changes in GNSS equipment (such as the receiver type, and the antenna brand and model) and the varying processing strategies (such as the software, and the mapping function and input parameter) in IGS ZTD product over time [59]. Even so, the fluctuations are not large, with RMS values fluctuating around 0.5 mm. It is noteworthy that the higher the spatial resolution of reanalyses is, the higher the accuracy of reanalyses is. This can be explicable by the fact that the better performance of the reanalysis models is attributed not only to the spatial resolution, but also to underlying better physical models and assimilation methods [24,51].

4. PWV Means, Variability and Trend

4.1. PWV Means and Variability

Figure 8 shows the spatial distribution of PWV means derived from ERA5, MERRA-2, JRA-55 and NCEP/DOE reanalysis in Antarctica from 2001 to 2018. The mean PWV derived from the eleven radiosonde and seven GNSS sites (the GNSS stations VESL and PALM are excluded because only ten years of data are available) are also presented for comparison. Their PWV interannual variability are also calculated (standard deviation of the PWV series divided by its mean) and presented in Figure 9. From Figure 8, the spatial distributions of PWV means derived from ERA5, MERRA-2, JRA-55 and NCEP/DOE reanalysis data sets are consistent with the radiosonde-derived and GNSS-derived PWV, and they reproduce the spatial variability well compared to the radiosonde and GNSS data. The PWV spatial distribution patterns reproduced by these four reanalyses are generally similar, but ERA5, MERRA-2 and JRA-55 show more details of PWV spatial distribution than NCEP/DOE due to their higher spatial resolution, including sharper gradients in PWV. For example, on the eastern and western flanks of Antarctica, in the regions of steep orography (for instance, along the Ellsworth Mountains, which are the highest mountain ranges in Antarctica), NCEP/DOE fails to reproduce these details. The annual mean PWV content varies in the range of 0–11 mm in Antarctica. Larger mean PWV content occurs in the ocean around Antarctica, whereas a lower average PWV content is observed in the Antarctic continent with a range of 0–7 mm. Particularly for the Eastern Antarctic Plateau, the mean PWV content is less than 4 mm due to the higher terrain and desert regions. Figure 9 shows that the PWV variability derived from ERA5, MERRA-2 and JRA-55 reanalysis has similar patterns except for some small discrepancies, and all agree well with the radiosonde and GNSS data. In Antarctica, the PWV keeps at a low level all year round and its interannual variability is small in most regions, less than 10%. Furthermore, the interannual variability of PWV generally shows slightly stronger over the land with respect to the ocean. However, compared with other reanalyses, NCEP/DOE is observed a huge discrepancy in East Antarctica, with very strong PWV variability (exceeding 16%), which is unrealistic.

The spatial distributions of seasonal PWV means for the four seasons DJF, MAM, JJA and SON over Antarctica are shown in Figure 10, their corresponding PWV interannual variability are presented in Figure 11. Statistics of the mean absolute difference (MAD) between reanalyses and radiosonde/GNSS in seasonal PWV means and relative variability can be found in

Table 4. Statistical results of the mean absolute differences in PWV mean and relative variability for the four seasons at the eleven radiosonde and seven GNSS stations over Antarctica.

Comparison	Differences	Seasons	ERA5	MERRA-2	JRA-55	NCEP/DOE
Radiosonde	PWV mean (mm)	DJF	0.32	0.35	0.43	1.35
		MAM	0.21	0.16	0.22	0.79
		JJA	0.13	0.10	0.18	0.62
		SON	0.24	0.22	0.25	0.69
		Average	0.22	0.21	0.27	0.86
	PWV relative variability (%)	DJF	3.10	3.53	3.49	3.89
		MAM	1.80	1.84	1.90	2.43
		JJA	1.94	1.76	2.43	2.30
		SON	2.27	2.44	2.14	3.01
		Average	2.28	2.39	2.49	2.91
GNSS	PWV mean (mm)	DJF	0.98	0.57	0.73	0.59
		MAM	1.13	0.78	0.88	0.50
		JJA	1.01	0.78	0.93	0.39
		SON	1.15	0.89	1.09	0.48
		Average	1.07	0.75	0.91	0.49
	PWV relative variability (%)	DJF	1.27	1.17	1.67	1.07
		MAM	1.79	2.26	3.02	3.51
		JJA	3.56	4.16	3.95	4.37
		SON	2.56	2.17	3.22	3.15
		Average	2.30	2.44	2.97	3.03

.

For seasonal PWV means, all the reanalyses reasonably reproduce the mean PWV spatial distributions revealed by the radiosonde and GNSS observations as shown in Figure 10, with most of MAD values below 1 mm during the four seasons (see Table 4). In addition, the patterns of seasonal means between the four reanalyses are also generally similar. The mean PWV is 0−12 mm over Antarctica with small seasonal variations during the four seasons. In DJF (summer), the PWV content is higher than in other seasons, whereas it is lowest in JJA (winter). This phenomenon results from the fact that most Antarctic regions enjoy the polar day in summer and the round-the-clock sunlight increases the temperature and the water vapor pressure and accelerates the hydrological cycle and melting of surface snow and sea ice, causing an increase in water vapor. In winter, most Antarctic regions are in the polar night, and the low temperature decreases the surface evaporation and water vapor content [22]. In addition, the interseasonal variation of PWV content in East Antarctica is small. Larger PWV interseasonal variations are mainly reproduced in West Antarctica and the peripheral sea. This is mainly due to the extremely cold and dry atmosphere in East Antarctica, where the altitude is relatively high.

For the PWV interannual variability, it should be noted that outliers with negative PWV were removed from NCEP/DOE for quality assurance in the calculation of relative variability. Even so, compared with the other three reanalyses, NCEP/DOE has the largest average MAD in interannual variability (see Table 4). In addition, it also shows unrealistic PWV interannual variations in East Antarctica for the four seasons, especially for MAM, which even reached 96.92% (see Figure 11). The reason could be that few observations were available for assimilation in the NCEP/DOE model in Antarctica, and compared with other reanalyses, the satellite radiances are not assimilated in NCEP/DOE [34]. For the other three reanalyses, as shown in Figure 11, their spatial patterns are also generally similar over the four seasons except for a few small discrepancies. Strong PWV interannual variability can be found to occur mainly in JJA (winter), while relatively weak in DJF (summer). For the cold season, PWV is more sensitive to temperature changes than in the warm season since Antarctic climate is extreme cold and dry. In JJA, there are three high-value centers of more than 20% located over the Marie Byrd Land, Enderby Land and Wilkes Land in Antarctica. Among them, Marie Byrd Land has the strongest PWV interannual variability, which may be related to the West Antarctic mantle thermal anomaly beneath the Marie Byrd Land that is accelerating the glacier melting and water vapor changes [60]. These variations also reveal that more climatic events, such as sea level rise [22]. In DJF, the spatial distribution of the PWV relative variability is more uniform than in other seasons with values lower than 14%. Table 4 shows that ERA5 and MERRA-2 are in better agreement with radiosonde and GNSS in PWV variability than JRA-55 and NCEP/DOE over Antarctica, indicating that ERA5 and MERRA-2 are more realistic than other reanalyses and that the ERA5 model has the best performance with average MADs of 2.28% and 2.3% verified by radiosonde and GNSS observations, respectively. Nevertheless, a few sites show large differences compared to all reanalyses, e.g., radiosonde station 89009 and GNSS station MCM4 (see Figure 1). In general, most of the Antarctic regional features of PWV interannual variability derived from reanalyses are confirmed by radiosonde and GNSS observations. However, they are difficult to diagnose in site-sparse regions and sites with gross errors due to complex topography and the extremely cold climate in Antarctica. Such issues can be further confirmed by comparison with other PWV measurements from satellite products (such as MODIS or AIRS) or other collocated instruments (such as VLBI or DORIS) [61,62].

4.2. PWV Trends

Figure 12 shows the time series of ERA5, MERRA-2, JRA-55 and NCEP/DOE annual PWV averaged over Antarctica from 2001–2018. In addition, the time series averaged from radiosonde, GNSS and the four reanalyses data at the five collocated sites (see Figure 1) are shown in Figure 13, and their relative PWV linear trends, expressed as percentages, are listed in

Table 5. Statistical results of relative PWV trends (in % decade⁻¹) from radiosonde, GNSS, ERA5, MERRA-2, JRA-55 and NCEP/DOE data at the 5 collocated sites from 2001 to 2018. Significant trends are highlighted in bold values.

Date	Sites
	CAS1	DAV1	MAW1	MCM4	SYOG
Radiosonde	−3.24	−2.92	3.70	3.75	−1.12
GNSS	−4.25	−9.78	−4.49	0.98	4.87
ERA5	−5.21	−5.40	−0.65	2.56	−1.40
MERRA-2	−9.55	−7.44	−5.19	−6.33	−0.96
JRA-55	−13.85	−3.46	1.25	1.43	1.10
NCEP/DOE	−5.76	−2.62	−1.20	−3.12	−1.86

.

In Figure 12, most reanalyses show overall negative trends (drying) except ERA5 which shows an opposite (moistening) trend. The PWV trends estimated from ERA5, MERRA-2, JRA-55 and NCEP/DOE over the period 2001−2018 for Antarctica are 0.46% decade⁻¹, −2.31% decade⁻¹, −0.03% decade⁻¹ and −1.51% decade⁻¹, respectively, but are only statistically significant in MERRA-2 and NCEP/DOE. Such discrepancy may be caused by differences of surface–atmosphere processes, large-scale water vapor transport and data assimilation in the four reanalyses [24]. The PWV variations derived from the four reanalyses are basically the same, but NCEP/DOE annual PWV estimates are larger than the other three. The higher NCEP/DOE PWV estimates are attributed to a larger estimation over the ocean around Antarctica compared with other reanalyses (see Figure 8 and Figure 10). In the previous sections, NCEP/DOE has been shown to have a poorer performance compared with other PWV results; thus, the NCEP/DOE estimated trends are not as reliable as other trends from the three reanalyses. From Figure 13 and Table 5, significant drying trends are observed at stations CAS1 (−13.85~−4.25% decade⁻¹) and DAV1 (−9.78~−2.62% decade⁻¹) in both PWV data sets, except for radiosonde observations which are insignificant at all five collocated stations. For the other three stations, only GNSS data show a significant drying trend (−4.49% decade⁻¹) at station MAW1, and only MERRA-2 exhibits a significant drying trend (−6.33% decade⁻¹) at station MCM4. No significant drying or moistening trend is found at station SYOG in both PWV data sets. Higher annual PWV estimates are also found in NCEP/DOE and GNSS data compared with the other PWV data sets at several stations. Among the four reanalyses, ERA5 and MERRA-2 show a better consistency with radiosonde observations at all five collocated stations.

The spatial distributions of relative PWV linear trends for both reanalyses are shown in Figure 14. Figure 14 illustrates that increasing PWV trends are observed near the South Pole in both reanalyses, whereas downward trends are observed over the Wilkes Land and Ellsworth Land. Among the reanalyses, MERRA-2, JRA-55 and NCEP/DOE generally have overall negative trends in most Antarctic regions, while ERA5 generally shows an overall positive trend. Moreover, MERRA-2 and NCEP/DOE show stronger trends than ERA5 and JRA-55.

We further investigated the spatial distributions of seasonal PWV linear trends over Antarctica, and the results are shown in Figure 15. Statistics of the MAD between reanalyses and radiosonde/GNSS in PWV trends for the four seasons are listed in

Table 6. Statistical results of the mean absolute differences in relative PWV trends (in % decade⁻¹) for the four seasons at the eleven radiosonde and seven GNSS stations over Antarctica.

Comparison	Seasons	ERA5	MERRA-2	JRA-55	NCEP/DOE
Radiosonde	DJF	2.63	3.05	3.09	3.43
	MAM	1.54	1.61	1.79	2.36
	JJA	1.69	1.53	2.23	2.23
	SON	1.91	2.06	1.88	2.76
	Average	1.94	2.06	2.25	2.69
GNSS	DJF	3.44	3.49	4.11	4.27
	MAM	4.11	5.48	6.46	5.88
	JJA	8.31	8.49	9.11	9.59
	SON	2.85	5.42	5.82	4.93
	Average	4.68	5.48	6.37	6.17

. As shown in Figure 15, the PWV trends depend upon the season and region location over Antarctica. For NCEP/DOE, it presents very strong drying trends in East Antarctica in MAM and JJA, with most values greater than −38% decade⁻¹, even reaching its highest value of −75.55% decade⁻¹ in JJA, which is unrealistic with respect to the other PWV data sets. For the other three reanalyses, in DJF, Wilkes Land and the area around the Weddell Sea significant drying trends are exhibited, and JRA-55 shows a more significant drying trend than ERA5 and MERRA-2. In MAM, drying trends are mainly found in Wilkes Land and Marie Byrd Land, whilst ERA5 and JRA-5 are insignificant in Marie Byrd Land. Moreover, moistening trends are identified in Ross Ice Shelf only for EAR5 and JRA-55. In JJA, significant drying trends are observed in Dronning Maud Land, Wilkes Land and Marie Byrd Land for all three reanalyses, and MERRA-2 presents more significant drying trends than ERA5 and JRA-5. In contrast, the Antarctic Peninsula shows general moistening trends in JJA for all three reanalyses. In SON, drying trends are mainly found in West Antarctica, whilst moistening trends mainly occur over the East Antarctic Plateau.

Table 6 confirms that the PWV trends estimated from ERA5 and MERRA-2 generally agree better with radiosonde and GNSS than JRA-55 and NCEP/DOE over Antarctica. In addition, ERA5 performs best with average MADs of 1.94% decade⁻¹ and 4.68% decade⁻¹ verified by radiosonde and GNSS, respectively, whereas NCEP/DOE performs worst with average MADs of 2.69% decade⁻¹ and 6.17% decade⁻¹, respectively. Compared with radiosonde, the MADs of all four reanalyses for relative PWV trends are largest in DJF (summer), with values of 2.63%, 3.05%, 3.09% and 3.43% per decade for ERA5, MERRA-2, JRA-55 and NCEP/DOE, respectively. When compared with GNSS, the largest MADs are detected in JJA (winter): 8.31%, 8.49%, 9.11% and 9.59% per decade for ERA5, MERRA-2, JRA-55 and NCEP/DOE, respectively, which are about three times larger than those of radiosonde data. Such a seasonal discrepancy between radiosonde and GNSS may be related to how the performance of different equipment varies in different seasons, or changes in equipment. Moreover, the artificial errors of pressure provided by ERA5 introduced in GNSS PWV retrieval are relatively large in winter because of the poorer accuracy of pressure in winter with respect to summer for ERA5 [36].

It is important to point out that the PWV trends estimated from the four reanalyses and other PWV measurements are not always statistically significant and have a high uncertainty over Antarctica in most periods and regions, especially for NCEP/DOE, which has the poorest performance. In addition, only the linear PVW trends are estimated and compared in this work. The actual PWV trends are not simply linear, and they can vary dramatically in different periods and areas, thus estimations and comparisons in this section may not be valid for other different areas or different periods.

5. Discussion and Conclusions

Currently, since the accuracy of atmospheric reanalysis data sets has been improved significantly, they play a more important role in both climate and weather study. In this study, the PWV means, variability and trends over Antarctica in four reanalyses (ERA5, MERRA-2, JRA-55 and NCEP/DOE) are investigated using radiosonde and GNSS observations covering 18 years from 2001–2018.

PWV estimations derived from ERA5, MERRA-2, JRA-55 and NCEP/DOE have been evaluated by radiosonde and GNSS observations over Antarctica. The assessment results show that ERA5 and MERRA-5 present a better performance than JRA-55 and NCEP/DOE. Among the four reanalyses, ERA5 shows the best comprehensive performance whereas NCEP/DOE performs the worst. For the comparison with radiosonde, the mean RMS errors of ERA5, MERRA-2, JRA-55 and NCEP/DOE are 0.5 mm, 0.65 mm, 0.85 mm and 1.11 mm, respectively. Regarding the comparison with GNSS, ERA5, MERRA-2, JRA-55 and NCEP/DOE have mean RMS errors of 1.2 mm, 1.15 mm, 1.34 mm and 1.59 mm, respectively. The main reason for the better performance of ERA5 and MERRA-2 comes certainly from improvements in spatial resolution and data assimilation methods.

The annual and seasonal distributions in PWV means reproduced from the four reanalysis data sets both perform fairly well over Antarctica, and they are all consistent with the radiosonde and GNSS observations with mean differences less than 1.1 mm. The PWV contents over Antarctic continent are very low and less than surrounding areas, showing a decreasing trend from west to east. The highest PWV content appears in summer whilst the lowest content is observed in winter. More details and features in spatial distribution of PWV means are reproduced well by ERA5, MERRA-2 and JRA-5, while NCEP/DOE is not recommended due to its lower spatial resolution and artefacts. In terms of the PWV interannual variability, the annual PWV relative variability over Antarctica is small with values smaller 10%. The seasonal variability is strongest in the cold season and relatively weak in the warm season. All reanalyses depict interannual variability well over Antarctica with mean differences smaller than 3%, except for NCEP/DOE which shows spurious variability in East Antarctica. PWV trends for all reanalyses in most Antarctic regions in the period 2001–2018 are insignificant tested by the Mann–Kendall method with a 10% significance level, and most reanalyses show overall drying trends in annual averages except for ERA5, which exhibits a moistening trend of 0.46% decade⁻¹. PWV trends also vary with season and region over Antarctica. Most reanalyses generally agree with radiosonde and GNSS in PWV trends over Antarctica with mean differences smaller 6.4% decade⁻¹, except for NCEP/DOE which also generates spurious trends in East Antarctica. Overall, the performances of ERA5 and MERRA-2 are the best among these reanalyses, and ERA5 performs slightly better than MERRA-2. In conclusion, NCEP/DOE is not recommended for use in Antarctica.

The results drawn from this work provide a complete investigation of the water vapor climatology over Antarctica for the last 18 years, thanks to the rapid development of global reanalysis data sets, especially for the new reanalysis ERA5. It also gives the reader an understanding of the reliability of PWV derived from the four reanalysis data sets over Antarctica. However, decadal variability and trends in PWV are not indicative of the long-term atmospheric climate change since several studies [1,24,25] suggested that a longer time series (nearly three decades) may be needed to acquire full stable patterns of trends. In addition, due to the fact that few observations are available, PWV variability and trends in Antarctica still suffer from large uncertainties and more in situ observations are needed for further investigation.