# Precipitable Water Vapor Retrieval and Analysis by Multiple Data Sources: Ground-Based GNSS, Radio Occultation, Radiosonde, Microwave Satellite, and NWP Reanalysis Data

Authors: Qin Zhang; Junhua Ye; Shuangcheng Zhang; Fei Han
Journal: Journal of Sensors
 (2018)

---

## Abstract

Precipitable water vapor (PWV) content detection is vital to heavy rain prediction; up to now, lots of different measuring methods and devices are developed to observe PWV. In general, these methods can be divided into two categories, ground-based or space-based. In this study, we analyze the advantages and disadvantages of these technologies, compare retrieved atmosphere parameters by different RO (radio occultation) observations, like FORMOSAT-3/COSMIC (Formosa Satellite-3 and Constellation Observing System for Meteorology, Ionosphere, and Climate) and FY3C (China Feng Yun 3C), and assess retrieved PWV precision with a radiosonde. Besides, we interpolate PWV from NWP (numerical weather prediction) reanalysis data for more comparison and analysis with RO. Specifically, ground-based GNSS is of high precision and continuous availability to monitor PWV distribution; in our paper, we show cases to validate and compare GNSS retrieving PWV with a radiosonde. Except GNSS PWV, we give two different radio occultation sounding results, COSMIC and FY3C, to validate the precision to monitor PWV from space in a global area. FY3C results containing Beidou (China Beidou Global Satellite Navigation System) radio occultation events need to be emphasized. So, in our study, we get the retrieved atmospheric profiles from GPS and Beidou radio occultation observations and derive atmosphere PWV by a variational retrieval method based on these data over a global area. Besides, other space-based methods, such as microwave satellite, are also useful in detecting PWV distribution situations in a global area from space; in this study, we present a case of retrieved PWV using microwave satellite observation. NWP reanalysis data ECMWF (European Centre for Medium-Range Weather Forecasts) ERA-Interim and the new-generation reanalysis data ERA5 provide global grid atmosphere parameters, like surface temperature, different-level pressures, and precipitable water. We show cases of retrieved PWV and validate the precision with radiosonde results and compare new reanalysis dataset ERA5 with ERA-Interim, finding that ERA5 can get higher precision-retrieved atmosphere parameters and PWV. In the end, from our comparison, we find that the retrieved PWV from RO (FY3C and COSMIC) and ECMWF reanalysis data (ERA-Interim and ERA5) have a high positive correlation and that almost allR2 values exceed 0.9, compare retrieved PWV with a radiosonde, and find that whether it is RO and ECMWF reanalysis data, ground-based GNSS, or microwave satellite, they all show small biases.


---

## Body

## 1. Introduction

Precipitable water vapor content change plays an important role in atmospheric water transport, energy conversion, and climate change, so it is meaningful to monitor and understand its change situation and mechanism [1, 2]. Nowadays, lots of methods and platforms can be used to detect atmospheric PWV; in general, there are two categories: ground-based and space-based devices or sensors. Ground-based methods, like radiosonde, sun photometry, and microwave radiometry, and space-based observation methods, such as COSMIC radio occultation, moderate-resolution imaging spectroradiometry (MODIS), SCIAMACHY, and atmospheric infrared sounder (AIRS), have the capability to detect PWV, for example, MODIS which provides daytime and night PWV products over land and ocean areas [3–9]. Traditional atmosphere sensing techniques such as radiosonde and microwave radiometer have shortcomings in reflecting the continuous transformations of atmospheric PWV due to their low, inhomogeneous spatial distribution and insufficient time availability. But the radiosonde has been thought to be a useful reference tool for atmospheric sounding [10]. In contrast, the ground-based GNSS technique has several advantages: it has been recognized as an efficient approach to estimate PWV change situations in nearly real time, and it can keep tracking and sensing continuously [11–16]. However, GNSS stations, for example, IGS stations, still have low spatial resolution, especially in the ocean area. Some other ground-based PWV detection technologies, like sun photometry and microwave radiometry, are used somewhere, but actually, due to the cost and observation condition limitation reasons, these devices cannot be used extensively, and it is impossible to get a continuous and wide range of observations. In these years, space-based methods, like radio occultation and microwave satellite-borne passive microwave detectors, provide availability to measure meteorology parameters and PWV distribution from space [3, 17, 18]; they have lots of advantages compared to ground observation methods. Like RO, it has high vertical resolution and RO measurements are not significantly affected by clouds and precipitation [19–22]; besides, its observation coverage is global, not a matter of land or ocean area, so radio occultation is applied to obtain PWV in a global area under all-weather conditions [23–25]. Precision statistical related studies reveal over ocean-dominated geographical areas; PWV retrieved from RO and ground-based GNSS exhibits a global mean difference of around 1 mm, a root-mean-square deviation of about 5 mm, and a correlation above 0.9 [26]. ECMWF and National Centers for Environmental Prediction (NCEP) have developed some numerical weather prediction models based on historical multiple data sources and provided a global reanalysis dataset, such as the most used ERA-Interim. ERA5 is a new-generation ECMWF-retrieved reanalysis product; it has higher temporal and spatial resolutions [27]. ECMWF analysis represents optimal humidity estimates from high-quality observations among multisatellite sounders, imagers, and conventional in situ observations through a data assimilation system [28–30]; researchers have validated global ERA-Interim integrated precipitable water using ground-based GNSS observations and related techniques, which showed a high correlation [8]. Except radio occultation, microwave satellite infrared observations can also be used to provide PWV distribution and change situations.Up to now, the radiosonde and radiometer are mostly thought to be efficient atmospheric observation techniques. In this study, we present some comparison and analysis work between these observations and reference results obtained from the radiosonde; besides, we derive PWV from the IGS GNSS station’s observations and present PWV of ATVOS measurements; what is more, we retrieve PWV from ERA-Interim and ERA5 reanalysis datasets and perform a global comparison of PWV from the collocated COSMIC and FY3C radio occultation profiles. We compare these results with the radiosonde to check their applicability and precision. In this paper, background information and methodology are described in Section2. Analysis, comparison, and discussion of PWV results from these data sources ground-based GNSS, COSMIC and FY3C RO, microwave satellite, ECMWF ERA-Interim, ERA5 and radiosonde are presented in Section 3. Conclusions are presented in Section 4.


## 2. Data Source and Methodology

### 2.1. Ground-Based GNSS

GNSS is widely used in positioning or navigation areas; in recent years, several regional or global satellite navigation systems have been built. In relation to this, lots of ground GNSS observation networks are running for multiple scientific tests, such as EUREF Permanent GNSS Network (EPN), US Constellation Observation System for Meteorology, Ionosphere, and Climate (SuomiNet), and Hong Kong Satellite Positioning Reference Station Network (SatRef). These GNSS networks provide sufficient and continuous observations for space monitoring. We all know that raw GNSS observations contain different types of errors due to signal transmission from satellite to receiver, so researchers try to use GNSS observations and related data processing algorithms to get troposphere delay, then combine meteorology parameters like pressure and temperature to get a PWV distribution situation [31]. The troposphere delay effect on the GNSS signal can be divided to a hydrostatic part and a wet part by
(1)ZTD=ZHD+ZWD,where ZTD is the zenith total delay, ZHD is the zenith hydrostatic delay, and ZWD denotes the zenith wet delay. Usually, ZHD can be calculated with high accuracy based on the Saastamoinen model, which can be expressed as
(2)ZHD=0.0022768P1−0.00266cos2φ−0.00028H,where P,∅, and H represent the station total pressure and station height, respectively. PWV is related to ZWD via a conversion factor by
(3)PWV=π∗ZWD,(4)π=106ρwRvK3/Tm+K2,,where π is the conversion factor, ρw is the density constant of liquid water (103kg/m3), Rv is the gas constant for PWV (461J·kg−1·K−1), K2′ and K3 are the atmospheric refractivity constants, where K2′=16.48K·hPa−1 and K3=3.776±0.014∗105K2·hPa−1, respectively, and Tm is the weighted mean temperature of the troposphere.Tm can be expressed as
(5)Tm=∫e/Tdz∫e/T2dz,where e is the precipitable water pressure, T is the absolute temperature, and dz is the integral path. Generally, it can be calculated by empirical equation from the Bevis formula using the ground surface temperature Ts, as follows:
(6)Tm=a+bTs.


### 2.2. Space-Based GNSS Radio Occultation

The GNSS radio occultation technique provides wealthy data for monitoring the global atmosphere. RO measures the time delay in the occulted signal using the GNSS receiver on board a low Earth orbiting (LEO) satellite, which can be transformed to atmospheric bending angle and refractivity; then, the vertical profiles of atmospheric pressure, temperature, and PWV are obtained. The GNSS radio occultation sounding technique provides the other possibility to monitor air PWV distribution from space. Up to now, there are several radio occultation satellite constellations, such as CHAllenging Minisatellite Payload (CHAMP), Gravity Recovery and Climate Experiment (GRACE), COSMIC, and FY3C. In this paper, we use COSMIC and FY3C radio occultation observations to retrieve PWV, so here we introduce more about COSMIC and FY3C constellations; the COSMIC radio occultation mission had been successfully launched on April 14, 2006. Receivers were installed on board the six small FORMOSAT-3/COSMIC satellites to capture the phase and amplitude of radio waves at two GPS frequencies. Termed as the FORMOSAT-3/COSMIC mission [32], the new constellation’s primary science goal is to obtain the near real-time vertical profiles of temperature, pressure, refractivity, and PWV in the neutral atmosphere and electron density in the ionosphere with global coverage at different altitudes [3, 33]. The measurements during five years of mission existence will provide about 2500 soundings per day, thus generating an extensive body of information to support operational global weather prediction, climate change monitoring, ionospheric phenomena, space weather research, and estimations of connections of meteorological and ionospheric processes with solar activity and human impact. The FY3C satellite was launched at 03:07 UTC on September 23, 2013, from the Taiyuan Satellite Base, Shanxi province, China. The Global Navigation Satellite System Occultation Sounder (GNOS) is installed in Feng-Yun 3 (FY3) satellites. The FY3C/GNOS is capable of tracking the occultation signal of the Beidou from space for the first time.Assuming that the effects of ionization have been completely corrected and that scattering due to water droplets is negligible, refractivity is related to temperature (T), partial pressure of dry air (Pd), and precipitable water (e) through
(7)N=K1∗PdZdT+K2∗eZwT2+K3∗eZwT,where K1=77.643K/hpa, K2=3.75463∗105K2/hpa, and K3=71.2952K/hpa are empirically determined constants, and Zd and Zw are the dry air and wet air compressibility factors. The compressibility factors correct for nonideal gas effects. Zd and Zw may be up to 0.05% smaller than unity in the denser parts of the atmosphere. Below 7~8 km altitude, PWV is nonignorable; formula (7) cannot resolve atmosphere temperature and PWV simultaneously. Raw measurements of GNSS radio occultation can be used to retrieve atmospheric parameters such as temperature, pressure, and humidity. The one-dimensional variational (1D-Var) method provides an effective way to combine observations and background information that are retrieved from ECMWF forecast datasets [29, 34, 35]. The method consists of finding the most probable atmospheric state by minimizing a cost function (J):
(8)Jx=hx−y0TO+F−1hx−y0+x−xbB−1x−xb,where y0 donates the observation vector; h is the observation operator (nonlinear); xb is the background reference information; O,F, and B are the error covariance matrices of observations, observation operator, and background information, respectively; and hx is an estimate of observations that would be made with a state of atmosphere x. The minimum variance problem can be solved by Quasi-Newton iteration.As a result, we get statistically optimal pressure, temperature, and humidity profiles and their respective errors. The vertical integral of the absolute vapor mass density yields the PWV per square meter, and it can also be converted to the unit of millimeter or centimeter by dividing by the density of liquid water, 1g/cc; the following is the PWV formula:
(9)PWV=∫0PsqgdP,where N denotes the refractivity index; K1,K2, and K3 are constants; e expresses the precipitable water pressure; Zd is the dry delay; T is the temperature; Zw is the wet delay; Pd represents dry air pressure; and PWV is the precipitable water vapor: PWV (kg/m2or mm/cm).The vertical integral of the humidity, both vapor and condensed moisture, yields the column water perm2. The total atmospheric PWV contained in a vertical column of unit cross-sectional area extends from the earth’s surface to the “top” boundary of the atmosphere. It is generally expressed in terms of the height to which the water substance would stand if completely condensed and collected in a vessel of the same unit cross section. Mathematically, PWV contained in a layer bounded by pressures p1 and p2 is given by this formula:
(10)PWV=∫Z0=0Zqhzρmzdz.qhz denotes the specific humidity in elevation z, and ρmzdenotes the moist air density in elevation z. The atmosphere is approximately divided to lots of layers from surface to top, so the PWV formula can be converted to follow the equation, expressed by pressure:
(11)PWV=−1g∫pZ0=0pZqhpdz≅−1g∑i=1Nq¯h,i∆p.The average specific humidity is calculated by layer bottom–specific humidity and layer top–specific humidity:(12)q¯h,i=qh,i−1+qh,i2.COSMIC wet atmosphere profiles provide multiple levels of precipitable water pressure data, so we can use the following formula to calculate specific humidity:(13)q=ε∗ep−e∗1−ε,where ε is a constant, which equals 0.622g.kg−1, r is the mixing ratio of precipitable water, e is the precipitable water pressure, ε is the ratio of the molecular weight of precipitable water to dry air, p is the total atmospheric pressure, and q is the specific humidity.


### 2.3. NWP Reanalysis Data

ERA-Interim is a global atmospheric reanalysis data from 1979 to the present, produced by a numerical weather prediction model run at ECMWF. The horizontal resolution of ERA-Interim is0.75°×0.75°, with a temporal resolution of 6 h [36]. ERA-Interim is one of the most advanced global atmospheric reanalysis data which represents the state of the atmosphere using the 4D-Var method [37] and assimilates a number of different sources of observations such as radiosonde humidity, atmospheric infrared sounder radiance (AIRS), GNSS RO-bending angle profiles, and Special Sensor Microwave/Imager (SSM/I), ERS- (European Remote Sensing Satellite-) 1 and -2 [36]. Here, we use monthly ERA-Interim data with a horizontal resolution of 1.5°×1.5° on 37 pressure levels between 1 hPa and 1000 hPa [36]. The spatial resolution of the dataset is approximately 80 km (T255 spectral) on 60 vertical levels from the surface up to 0.1 hPa [29]. This model is one of the most advanced in operational use and capable of predicting the global atmosphere with accuracy just barely less than what is theoretically possible [37]. ERA5 will be the fifth generation of ECMWF atmospheric reanalysis of the global climate, which starts with the reanalysis produced in the 1980s, followed by ERA-15, ERA-40, and the most used ERA-Interim. The new ERA5 reanalysis will span the modern observing period from 1979 onward, with daily updates continuing forward in time. ERA5 will eventually replace ERA-Interim. ERA5 data will be at a much higher resolution than ERA-Interim; hourly analysis fields will be available at a horizontal resolution of 31 km on 137 levels, from the surface up to 0.01 hPa (around 80 km). ERA-Interim and ERA5 both provide surface and multiple pressure levels of data; users can select different types of atmosphere datasets, such as temperature, pressure, or column precipitable water profiles. In the paper, we use the total column precipitable water data from 0.75°×0.75° 6 h ERA-Interim reanalysis and 0.25°×0.25° 3 h ERA5 datasets.


### 2.4. Microwave Satellite

Based on satellite-borne passive microwave detectors, it is possible to observe the atmospheric PWV over the vast oceanic regions where only limited ship-based and buoy-based observations were previously available [38]. Microwaves are less affected by scattering and absorption from cloud, allowing for PWV measurements in most weather conditions. Unlike the infrared measurements, the passive microwave detector is able to measure the atmospheric PWV content under both cloud-free and cloudy conditions since the microwave can penetrate cloud. Lots of studies had presented that the data sources from satellite measurements can be used to retrieve PWV distribution with varying accuracy; for example, researchers use microwave satellite data to retrieve PWV in polar winter conditions and analyze the global PWV trend and its diurnal asymmetry based on GPS, radiosonde, and microwave satellite measurements. The High-Resolution Infrared Radiation Sounder (HIRS) unit in the TIROS-N Operational Vertical Sounder (TOVS) package is designed to sense the amount of PWV in the upper troposphere with high quality [11]. In this paper, we use TOVS-detected data for PWV retrieval and analysis.


### 2.5. Radiosonde

Radiosonde balloon measurements are launched globally, although with sparse coverage in lots of areas, such as over oceans or in the Southern Hemisphere. Most radiosonde locations are in land; in addition, there are also some radiosonde observations made on ships over oceans, but actually it is not easy and convenient for users. Radiosonde can get high precision atmosphere observations and has high vertical resolution, and its observation covers the range from the ground surface to an altitude of 30 km. However, the data quality varies strongly depending on the sensor type; now, many different types of sensors are used globally, and each one has its unique known and unknown biases. Besides, sensor types at different locations change over time, which can lead to artificial trends or jumps. Most stations provide twice observations in one day, one is in 00 : 00 and the other one is in 12 : 00, while some stations can also provide four times observations per day. From radiosonde observations, we can get pressure, temperature, relative humidity, wind direction, and wind speed these parameters, and then we can use formula (13) to calculate the PWV in each site.


### 2.6. PWV Horizontal and Vertical Interpolation

When we do PWV comparison and precision evaluation between different datasets, the direct issue is to match them to the same site, which means we need to do horizontal and vertical interpolation. In our paper, we use five different types of datasets: radio occultation, ECMWF reanalysis data, GNSS observation, microwave satellite, and radiosonde. All these data, in horizontal, often express location using geographic coordinates latitude/longitude. If we want to interpolate PWV to RO sites using ERA-Interim reanalysis data, firstly we need to pick up the total column PWV values of four grid points surrounding RO location and then calculate the comprehensive value using the following equations:(14)PWV=∑i=14WiPi,(15)Wi=1/di∑j=141/di.di denotes the distance grid point to the target point, Wi is the weight of each grid point PWV value, Pi is the grid point retrieving PWV, and PWV is the end target point PWV.Actually, different datasets use multiple height systems, like ERA-Interim or ERA5, geopotential height is used to express different layer observations, the GNSS station often uses geodetic height, and the radiosonde uses altitude. However, in our paper, we need the total column precipitable water vapor, which means from ground surface to top no matter what kind height systems they used.


### 2.7. Evaluation

In the paper, we use mean bias and RMSE (root-mean-square error) to analyze and evaluate the retrieved PWV precision of different types of measurements.(16)Meanbias=∑i=1nxobs,i−xmodel,in,(17)RMSE=∑i=1nxobs,i−xmodel,i2n.


## 2.1. Ground-Based GNSS

GNSS is widely used in positioning or navigation areas; in recent years, several regional or global satellite navigation systems have been built. In relation to this, lots of ground GNSS observation networks are running for multiple scientific tests, such as EUREF Permanent GNSS Network (EPN), US Constellation Observation System for Meteorology, Ionosphere, and Climate (SuomiNet), and Hong Kong Satellite Positioning Reference Station Network (SatRef). These GNSS networks provide sufficient and continuous observations for space monitoring. We all know that raw GNSS observations contain different types of errors due to signal transmission from satellite to receiver, so researchers try to use GNSS observations and related data processing algorithms to get troposphere delay, then combine meteorology parameters like pressure and temperature to get a PWV distribution situation [31]. The troposphere delay effect on the GNSS signal can be divided to a hydrostatic part and a wet part by
(1)ZTD=ZHD+ZWD,where ZTD is the zenith total delay, ZHD is the zenith hydrostatic delay, and ZWD denotes the zenith wet delay. Usually, ZHD can be calculated with high accuracy based on the Saastamoinen model, which can be expressed as
(2)ZHD=0.0022768P1−0.00266cos2φ−0.00028H,where P,∅, and H represent the station total pressure and station height, respectively. PWV is related to ZWD via a conversion factor by
(3)PWV=π∗ZWD,(4)π=106ρwRvK3/Tm+K2,,where π is the conversion factor, ρw is the density constant of liquid water (103kg/m3), Rv is the gas constant for PWV (461J·kg−1·K−1), K2′ and K3 are the atmospheric refractivity constants, where K2′=16.48K·hPa−1 and K3=3.776±0.014∗105K2·hPa−1, respectively, and Tm is the weighted mean temperature of the troposphere.Tm can be expressed as
(5)Tm=∫e/Tdz∫e/T2dz,where e is the precipitable water pressure, T is the absolute temperature, and dz is the integral path. Generally, it can be calculated by empirical equation from the Bevis formula using the ground surface temperature Ts, as follows:
(6)Tm=a+bTs.


## 2.2. Space-Based GNSS Radio Occultation

The GNSS radio occultation technique provides wealthy data for monitoring the global atmosphere. RO measures the time delay in the occulted signal using the GNSS receiver on board a low Earth orbiting (LEO) satellite, which can be transformed to atmospheric bending angle and refractivity; then, the vertical profiles of atmospheric pressure, temperature, and PWV are obtained. The GNSS radio occultation sounding technique provides the other possibility to monitor air PWV distribution from space. Up to now, there are several radio occultation satellite constellations, such as CHAllenging Minisatellite Payload (CHAMP), Gravity Recovery and Climate Experiment (GRACE), COSMIC, and FY3C. In this paper, we use COSMIC and FY3C radio occultation observations to retrieve PWV, so here we introduce more about COSMIC and FY3C constellations; the COSMIC radio occultation mission had been successfully launched on April 14, 2006. Receivers were installed on board the six small FORMOSAT-3/COSMIC satellites to capture the phase and amplitude of radio waves at two GPS frequencies. Termed as the FORMOSAT-3/COSMIC mission [32], the new constellation’s primary science goal is to obtain the near real-time vertical profiles of temperature, pressure, refractivity, and PWV in the neutral atmosphere and electron density in the ionosphere with global coverage at different altitudes [3, 33]. The measurements during five years of mission existence will provide about 2500 soundings per day, thus generating an extensive body of information to support operational global weather prediction, climate change monitoring, ionospheric phenomena, space weather research, and estimations of connections of meteorological and ionospheric processes with solar activity and human impact. The FY3C satellite was launched at 03:07 UTC on September 23, 2013, from the Taiyuan Satellite Base, Shanxi province, China. The Global Navigation Satellite System Occultation Sounder (GNOS) is installed in Feng-Yun 3 (FY3) satellites. The FY3C/GNOS is capable of tracking the occultation signal of the Beidou from space for the first time.Assuming that the effects of ionization have been completely corrected and that scattering due to water droplets is negligible, refractivity is related to temperature (T), partial pressure of dry air (Pd), and precipitable water (e) through
(7)N=K1∗PdZdT+K2∗eZwT2+K3∗eZwT,where K1=77.643K/hpa, K2=3.75463∗105K2/hpa, and K3=71.2952K/hpa are empirically determined constants, and Zd and Zw are the dry air and wet air compressibility factors. The compressibility factors correct for nonideal gas effects. Zd and Zw may be up to 0.05% smaller than unity in the denser parts of the atmosphere. Below 7~8 km altitude, PWV is nonignorable; formula (7) cannot resolve atmosphere temperature and PWV simultaneously. Raw measurements of GNSS radio occultation can be used to retrieve atmospheric parameters such as temperature, pressure, and humidity. The one-dimensional variational (1D-Var) method provides an effective way to combine observations and background information that are retrieved from ECMWF forecast datasets [29, 34, 35]. The method consists of finding the most probable atmospheric state by minimizing a cost function (J):
(8)Jx=hx−y0TO+F−1hx−y0+x−xbB−1x−xb,where y0 donates the observation vector; h is the observation operator (nonlinear); xb is the background reference information; O,F, and B are the error covariance matrices of observations, observation operator, and background information, respectively; and hx is an estimate of observations that would be made with a state of atmosphere x. The minimum variance problem can be solved by Quasi-Newton iteration.As a result, we get statistically optimal pressure, temperature, and humidity profiles and their respective errors. The vertical integral of the absolute vapor mass density yields the PWV per square meter, and it can also be converted to the unit of millimeter or centimeter by dividing by the density of liquid water, 1g/cc; the following is the PWV formula:
(9)PWV=∫0PsqgdP,where N denotes the refractivity index; K1,K2, and K3 are constants; e expresses the precipitable water pressure; Zd is the dry delay; T is the temperature; Zw is the wet delay; Pd represents dry air pressure; and PWV is the precipitable water vapor: PWV (kg/m2or mm/cm).The vertical integral of the humidity, both vapor and condensed moisture, yields the column water perm2. The total atmospheric PWV contained in a vertical column of unit cross-sectional area extends from the earth’s surface to the “top” boundary of the atmosphere. It is generally expressed in terms of the height to which the water substance would stand if completely condensed and collected in a vessel of the same unit cross section. Mathematically, PWV contained in a layer bounded by pressures p1 and p2 is given by this formula:
(10)PWV=∫Z0=0Zqhzρmzdz.qhz denotes the specific humidity in elevation z, and ρmzdenotes the moist air density in elevation z. The atmosphere is approximately divided to lots of layers from surface to top, so the PWV formula can be converted to follow the equation, expressed by pressure:
(11)PWV=−1g∫pZ0=0pZqhpdz≅−1g∑i=1Nq¯h,i∆p.The average specific humidity is calculated by layer bottom–specific humidity and layer top–specific humidity:(12)q¯h,i=qh,i−1+qh,i2.COSMIC wet atmosphere profiles provide multiple levels of precipitable water pressure data, so we can use the following formula to calculate specific humidity:(13)q=ε∗ep−e∗1−ε,where ε is a constant, which equals 0.622g.kg−1, r is the mixing ratio of precipitable water, e is the precipitable water pressure, ε is the ratio of the molecular weight of precipitable water to dry air, p is the total atmospheric pressure, and q is the specific humidity.


## 2.3. NWP Reanalysis Data

ERA-Interim is a global atmospheric reanalysis data from 1979 to the present, produced by a numerical weather prediction model run at ECMWF. The horizontal resolution of ERA-Interim is0.75°×0.75°, with a temporal resolution of 6 h [36]. ERA-Interim is one of the most advanced global atmospheric reanalysis data which represents the state of the atmosphere using the 4D-Var method [37] and assimilates a number of different sources of observations such as radiosonde humidity, atmospheric infrared sounder radiance (AIRS), GNSS RO-bending angle profiles, and Special Sensor Microwave/Imager (SSM/I), ERS- (European Remote Sensing Satellite-) 1 and -2 [36]. Here, we use monthly ERA-Interim data with a horizontal resolution of 1.5°×1.5° on 37 pressure levels between 1 hPa and 1000 hPa [36]. The spatial resolution of the dataset is approximately 80 km (T255 spectral) on 60 vertical levels from the surface up to 0.1 hPa [29]. This model is one of the most advanced in operational use and capable of predicting the global atmosphere with accuracy just barely less than what is theoretically possible [37]. ERA5 will be the fifth generation of ECMWF atmospheric reanalysis of the global climate, which starts with the reanalysis produced in the 1980s, followed by ERA-15, ERA-40, and the most used ERA-Interim. The new ERA5 reanalysis will span the modern observing period from 1979 onward, with daily updates continuing forward in time. ERA5 will eventually replace ERA-Interim. ERA5 data will be at a much higher resolution than ERA-Interim; hourly analysis fields will be available at a horizontal resolution of 31 km on 137 levels, from the surface up to 0.01 hPa (around 80 km). ERA-Interim and ERA5 both provide surface and multiple pressure levels of data; users can select different types of atmosphere datasets, such as temperature, pressure, or column precipitable water profiles. In the paper, we use the total column precipitable water data from 0.75°×0.75° 6 h ERA-Interim reanalysis and 0.25°×0.25° 3 h ERA5 datasets.


## 2.4. Microwave Satellite

Based on satellite-borne passive microwave detectors, it is possible to observe the atmospheric PWV over the vast oceanic regions where only limited ship-based and buoy-based observations were previously available [38]. Microwaves are less affected by scattering and absorption from cloud, allowing for PWV measurements in most weather conditions. Unlike the infrared measurements, the passive microwave detector is able to measure the atmospheric PWV content under both cloud-free and cloudy conditions since the microwave can penetrate cloud. Lots of studies had presented that the data sources from satellite measurements can be used to retrieve PWV distribution with varying accuracy; for example, researchers use microwave satellite data to retrieve PWV in polar winter conditions and analyze the global PWV trend and its diurnal asymmetry based on GPS, radiosonde, and microwave satellite measurements. The High-Resolution Infrared Radiation Sounder (HIRS) unit in the TIROS-N Operational Vertical Sounder (TOVS) package is designed to sense the amount of PWV in the upper troposphere with high quality [11]. In this paper, we use TOVS-detected data for PWV retrieval and analysis.


## 2.5. Radiosonde

Radiosonde balloon measurements are launched globally, although with sparse coverage in lots of areas, such as over oceans or in the Southern Hemisphere. Most radiosonde locations are in land; in addition, there are also some radiosonde observations made on ships over oceans, but actually it is not easy and convenient for users. Radiosonde can get high precision atmosphere observations and has high vertical resolution, and its observation covers the range from the ground surface to an altitude of 30 km. However, the data quality varies strongly depending on the sensor type; now, many different types of sensors are used globally, and each one has its unique known and unknown biases. Besides, sensor types at different locations change over time, which can lead to artificial trends or jumps. Most stations provide twice observations in one day, one is in 00 : 00 and the other one is in 12 : 00, while some stations can also provide four times observations per day. From radiosonde observations, we can get pressure, temperature, relative humidity, wind direction, and wind speed these parameters, and then we can use formula (13) to calculate the PWV in each site.


## 2.6. PWV Horizontal and Vertical Interpolation

When we do PWV comparison and precision evaluation between different datasets, the direct issue is to match them to the same site, which means we need to do horizontal and vertical interpolation. In our paper, we use five different types of datasets: radio occultation, ECMWF reanalysis data, GNSS observation, microwave satellite, and radiosonde. All these data, in horizontal, often express location using geographic coordinates latitude/longitude. If we want to interpolate PWV to RO sites using ERA-Interim reanalysis data, firstly we need to pick up the total column PWV values of four grid points surrounding RO location and then calculate the comprehensive value using the following equations:(14)PWV=∑i=14WiPi,(15)Wi=1/di∑j=141/di.di denotes the distance grid point to the target point, Wi is the weight of each grid point PWV value, Pi is the grid point retrieving PWV, and PWV is the end target point PWV.Actually, different datasets use multiple height systems, like ERA-Interim or ERA5, geopotential height is used to express different layer observations, the GNSS station often uses geodetic height, and the radiosonde uses altitude. However, in our paper, we need the total column precipitable water vapor, which means from ground surface to top no matter what kind height systems they used.


## 2.7. Evaluation

In the paper, we use mean bias and RMSE (root-mean-square error) to analyze and evaluate the retrieved PWV precision of different types of measurements.(16)Meanbias=∑i=1nxobs,i−xmodel,in,(17)RMSE=∑i=1nxobs,i−xmodel,i2n.


## 3. Results and Discussion

In this paper, mostly we focus on analyzing PWV retrieved by multiple data sources. We show cases and analyze the precision of ground-based GNSS and space-based monitoring methods, like microwave satellites and radio occultation. We try to use the retrieved atmospheric profiles of FY3C and COSMIC radio occultation soundings on the day of February 4, 2016 for PWV calculation and analysis. Besides, we also interpolate the “total column precipitable water vapor” data from ECMWF ERA-Interim and ERA5 for comparison and validation. Same as lots of related research, we still use the retrieved PWV from ground-based radiosonde as the reference, for precision analysis and validation. The following sections will show more details and analysis.


### 3.1. GNSS Retrieved PWV Comparison and Analysis

Lots of studies have been done for PWV content and distribution monitoring by the ground-based GNSS network. Especially in recent several years, the precise point positioning (PPP) technique developed rapidly; some researchers try to use the near real-time satellite orbit, clock, and ionosphere streams, for example the International Geodetic Service (IGS) RT-PP product, to monitor atmosphere zenith total delay, then combine atmospheric temperature and pressure observations to convert to PWV. In this paper, we just give a case to show the statistical dependence of retrieved PWV of IGS stations by GNSS and the nearby radiosonde’s observations [33, 39–41].We select the GNSS station SVTL for our test as well as its coordinates (latitude:60.53deg, longitude: 29.78deg, elevation: 60.64m) and location in Europe; the nearby radiosonde station ID is 26063, its coordinates are latitude: 59.95deg, longitude: 30.70deg, and elevation: 78.00m. The two types of observations are completely in 2017, so we choose the whole year’s observations to do analysis.GNSS can get high-frequency ZTD nearly real-time, which depends on receiver antenna and hardware process ability, which means it has a possibility of getting high-frequency PWV. In this case, the SVTL station provides retrieved PWV 30 min each. In the whole year of 2017, it gets 15,410 observation values; however, in each day, the radiosonde just measures twice, one at 00 : 00 and the other at 12 : 00 a.m., so each day we can get only two PWV values in the 26063 station. Totally, we get 730 PWV values in 2017, so it is easy to find GNSS providing much more PWV observations compared to the radiosonde.In the paper, so as to compare and analyze the retrieved PWV bias and accuracy by the two types of observations, we select each day’s 00 : 00 and 12 : 00 observations from GNSS PWV to do statistical analysis. The following is the distribution map of matched PWV values in each month of 2017.From Figures1 and 2, we can see in the whole year of 2017 GNSS retrieving PWV and radiosonde retrieving PWV; their correlation is up to 0.95. The max PWV of GNSS and RS are both produced between July and September, and the local season is in the summer. Table 1 lists the statistical information min bias, max bias, mean bias, and mean RMSE of matched RS retrieving PWV and GNSS retrieving PWV.Figure 1
The graphic shows SVTL GNSS station and the nearby radiosonde station 26063 PWV time series distributions of the whole year of 2017.Figure 2
PWV correlation analysis of the ground-based GNSS (SVTL) and the matched radiosonde (26063) in the whole year of 2017.Table 1
Precision statistic of retrieved PWV from ground-based GNSS station (STVL) compared with the nearby radiosonde.
Min bias (mm)
Max bias (mm)
Mean bias (mm)
RMSE (mm)
RS_PWV-GPS_PWV
−0.01
15.15
0.23
2.41From Table1, we can easily see that the mean bias is only 0.23 mm and the mean RMSE is 2.41 mm; they are very small, which means GNSS-retrieved PWV values have high precision. It is a useful method to monitor PWV situation changes.


### 3.2. Microwave Satellite Inferred PWV Analysis

We use the day January 1, 2014 (00 : 00), ATOVS (The Advanced TIROS Operational Vertical Sounder) data (Figure3) and radiosonde data to investigate and analyze PWV content distribution status and interpolate ATOVS-detected PWV in 589 radiosonde locations, so we can get the statistical results. In Table 2, the min bias is 0.003 mm, the max bias is 12.068 mm, the mean bias is −0.348 mm, and the mean RMS is 3.461 mm.Figure 3
Distribution of the PWV detected by microwave satellite on January 1, 2014 (00 : 00).Table 2
Precision statistic of retrieved PWV from the ATOVS product compared with radiosondes.
Min bias (mm)
Max bias (mm)
Mean bias (mm)
RMSE (mm)
RS_PWV-Sat_PWV
0.003
12.068
−0.348
3.461


### 3.3. PWV Retrieval from RO and NWP Reanalysis Data

#### 3.3.1. Distributions of COSMIC and FY3C RO Events

On February 4, 2016, COSMIC produced 891 radio occultation events and FY3C had 410. Figure4 shows the distribution map of these radio occultation event locations; we can see that these occultation event locations nearly cover the entire world from the North Pole to the South Pole.Figure 4
The two graphics show the locations of COSMIC radio occultation and FY3C radio occultation events produced on February 4, 2016.


#### 3.3.2. Retrieved Atmosphere Parameters of FY3C and COSMIC

In Figure5, we select two nearby FY3C and COSMIC radio occultation events as an example for comparison analysis; (a), (b), and (c) are the retrieved specific humidity profile, pressure profile, and temperature profile, respectively, of FY3C radio occultation on the site latitude:24.99°,longitude:−11.68° at 11 : 51 : 03 of the day February 4, 2016, so we find the nearest COSMIC radio occultation event regardless of the location but also the event time which occurred latitude:26.95°,longitude:−11.35° at 14 : 10 : 41 of February 4, 2016; (d), (e), and (f) list the COSMIC-retrieved specific humidity profile, pressure profile, and temperature profile, respectively. From (a) and (d), we can see that radio occultation detected specific humidity whether it is FY3C or COSMIC; in nearly ground height, the value exceeds 4k/kg and the height is between 2 km and 7 km. COSMIC detected that specific humidity has a bigger jump than FY3C, above 10 km, both of which detected a specific humidity reduction to 0, which means no PWV exists in the upper layer. From (b) and (e), the graphics show FY3C and COSMIC retrieving pressure from the ground surface to the upper atmosphere; FY3C shows a pressure of 0~60 km and COSMIC gives 0~40 km. We can see that above 20 km, the air pressure becomes very small; from the ground surface to 20 km, the pressure degrades rapidly. The max value of pressure is located at the ground surface, both showing no more than 1000 Mb; besides, the pressure variation tendencies of both pressure data are in substantial agreement. (c) shows FY3C retrieving temperature from the surface to 60 km; from the graphic, we can see that the highest temperature is in the ground surface, and then from the surface to 20 km, the temperature becomes lower and lower, with the lowest temperature being nearly 20 km in height. However, from 20 km, the temperature (205 K) increases, reaching 255 K in the 42 km height layer. From 42 km, the temperature decreases again until reaching 220 K in the 60 km layer. The same situation occurs in graphic (f); from the ground surface to 20 km, the temperature decreases, especially between the ground surface and the 12 km height layer, in which the temperature variation trend is nearly linearly decreasing; from 20 km to 40 km, the temperature becomes linearly increasing, although FY3C and COSMIC events have more than a 2 h gap, but as can be easily seen, their variation trends are basically the same.Figure 5
The above group graphics show the retrieved atmosphere profiles from FY3C and COSMIC radio occultation observations on February 4, 2016.
(a)
(b)
(c)
(d)
(e)
(f)


#### 3.3.3. Evaluation of COSMIC and FY3C Retrieved PWV

The paper uses radiosonde measurements as the reference standard to do analysis; due to this reason, radiosonde observations can only be measured at around 11 : 00 and 23 : 00. We use the following three rules to collect nearby matched RO and RS observations for comparison and do statistical analysis. In the first rule, we set latitude and longitude differences between radiosonde observation sites and radio occultation event which produced sites below1°; for the temporal resolution, we limit it to a ±2h difference. For example, if we use 11 : 00 radiosonde observations, the related radio occultation events must be produced at 9 : 00~13 : 00. That is to say, all of these observations we selected are produced in the time range 9 : 00~13 : 00 and 21 : 00~01 : 00 on February 4, 2016. For the second and third standards, we change the search range from latitude-longitude difference 1°/1° to 2°/2° and 3°/3° individually.On February 4, 2016, COSMIC produced 891 radio occultation events and FY3C had 410. We use the above three standards to collect matched data from FY3C, COSMIC radio occultations, and radiosonde observations in the global area. Statistical results are listed in Table3, whether it is FY3C- or COSMIC-retrieved PWV; compared to nearby radiosonde observations, their RMSE values are no more than 10 mm. From our statistical data, we find that RMSEs at 23 : 00 are all larger than 11 : 00, whether it is FY3C or COSMIC; it is not sure if RO measurements have big biases at 23 : 00 or RS produces large bias. In the future, we will do more deep research about this phenomenon.Table 3
Retrieved PWV precision statistic between FY3C, COSMIC, and radiosonde.
RO
Latitude/longitude difference (deg)
RMSE (mm)
11 : 00
23 : 00
FY3C
1/1
0.34
6.64
2/2
1.09
5.28
3/3
1.26
4.65
COSMIC
1/1
3.03
6.03
2/2
4.30
5.32
3/3
3.29
6.92


#### 3.3.4. PWV Retrieved from ECMWF ERA-Interim and ERA5

Figure6 shows the retrieved PWV distribution in the global area from ECMWF ERA-Interim and ERA5 in four periods (00 : 00, 06 : 00, 12 : 00, and 18 : 00); we can see that surrounding the equator, the PWV content exceeds 40 mm, and the PWV distribution decreases along the latitude from the equator to the South Pole to the North Pole.Figure 6
These graphics show the retrieved PWV in global by ECMWF ERA-Interim and ERA5 reanalysis data on February 4, 2016.


#### 3.3.5. Evaluation of ERA-Interim and ERA5-Retrieved PWV

In Figure7, it shows the retrieved PWV statistical coefficients between different data sources; we present the statistical analysis using two groups’ observations; respectively, one is in 00 : 00 and the other one uses 12 : 00 observations, so from these graphics, we can see that graph (a) shows the retrieved PWV coefficient between ERA5 and radiosonde at 00 : 00, with R2 equal to 0.9444. From (d), at 12 : 00, R2 equals 0.9708, so we can see that the retrieved PWV by ERA5 and radiosonde has a higher linear correlation. In (b) and (e), they show a retrieved PWV relationship of ERA-Interim and ERA5 at 00 : 00 and 12 : 00, and their R2 values are 0.9833 and 0.9325. Respectively, in theory both methods can get nearly the same results; ERA5 is little superior compared to ERA-Interim, due to the reason that it uses more data sources and new technology for assimilation analysis. (c) and (f) express the relationship of PWV between ERA-Interim and radiosonde, and their R2 values are 0.8699 and 0.9539, respectively, which means they have a higher linear correlation.Figure 7
The PWV correlation statistical of the ERA-Interim, ERA5, and RS. (a, b, c) 00 : 00 observations and (d, e, f) 12 : 00 observations.
(a)
(b)
(c)
(d)
(e)
(f)


#### 3.3.6. PWV Intercomparison between RO and ERA-Interim and ERA5

In Figure8, these graphics show the retrieved PWV distribution along the latitude: (a) shows obtained COSMIC-retrieved atmospheric profiles, (b) shows interpolated PWV in the site COSMIC radio occultation events that occurred based on ERA-Interim data, (c) shows interpolation in the same sites from ERA5, (d) shows obtained FY3C radio occultation–retrieved atmospheric profiles, and (e)-(f) show interpolated PWV in FY3C radio occultation events producing sites from ERA-Interim and ERA5 data. From the above graphics, we can easily see that the PWV has symmetrical distribution along the equator; from latitude:0° to the North Pole or South Pole, the retrieved PWV decreases to nearly 0 in the South or North Pole. The biggest PWV is detected in the equator area which exceeds 60g/m2. In (a), (b), and (c), these graphics show that they have nearly the same variation tendency in these sites by different data sources; (d), (e), and (f) have the same situations.Figure 8
The six graphics show the retrieved PWV in global along the latitude by COSMIC radio occultation (a), FY3C radio occultation (d), ERA-Interim (c, f), and ERA5 (b, e) reanalysis data on February 4, 2016. (a, b, c) Retrieved PWV from COSMIC RO and interpolated total column water vapor from ECMWF ERA-Interim and ERA5 reanalysis data on COSMIC radio occultation sites. (d, e, f) Retrieved PWV from FY3C RO and interpolated total column water vapor from ECMWF ERA-Interim and ERA5 reanalysis data on FY3C radio occultation sites.
(a)
(b)
(c)
(d)
(e)
(f)We calculate the mean bias and mean RMSE of COSMIC and FY3C-retrieved PWV with interpolated PWV from ERA-Interim and ERA5 datasets; these results are listed in Table4.Table 4
Retrieved PWV precision statistic between FY3C, COSMIC, and ERA-Interim, ERA5.
Mean bias (mm)
Mean RMSE (mm)
COSMIC RO sites
COSMIC-Interim
−2.3
4.5
COSMIC-ERA5
−2.5
4.8
Interim-ERA5
−0.1
2.0
FY3C RO sites
FY3C-Interim
0.9
3.6
FY3C-ERA5
1.0
3.8
Interim-ERA5
0.1
1.9As can be seen from Table4, the mean bias of COSMIC-retrieved PWV and ERA-Interim is −2.3 mm, RMSE is 4.5 mm, compared to COSMIC-ERA5 results, the statistical value has only −0.2 mm bias, and RMSE has 0.3 mm bias; the same situation occurs in FY3C observations, but on the whole, FY3C-retrieved PWV has smaller mean bias values and mean RMSE than COSMIC-retrieved results. COSMIC has 891 observations, and FY3C contains 410 observations; whether in COSMIC radio occultation sites or FY3C, ERA-Interim and ERA5 have stable statistical results, and the mean bias is around 0 and RMSE is approximately 2.0 mm.In the paper, we also analyze the bias distribution of radio occultation and ECMWF reanalysis data; as can be easily seen in Figure9, the most departures of COSMIC-retrieved PWV and ERA-Interim or ERA5 are no more than 10 cm, and with the same situation to FY3C, most time departures are below 10 mm. For the two data sources of ECMWF, the departures between ERA-Interim and ERA5 mostly are no more than 5 mm.Figure 9
The above graphics show the departures of the retrieved PWV by different observations.
(a)
(b)
(c)
(d)
(e)
(f)


## 3.1. GNSS Retrieved PWV Comparison and Analysis

Lots of studies have been done for PWV content and distribution monitoring by the ground-based GNSS network. Especially in recent several years, the precise point positioning (PPP) technique developed rapidly; some researchers try to use the near real-time satellite orbit, clock, and ionosphere streams, for example the International Geodetic Service (IGS) RT-PP product, to monitor atmosphere zenith total delay, then combine atmospheric temperature and pressure observations to convert to PWV. In this paper, we just give a case to show the statistical dependence of retrieved PWV of IGS stations by GNSS and the nearby radiosonde’s observations [33, 39–41].We select the GNSS station SVTL for our test as well as its coordinates (latitude:60.53deg, longitude: 29.78deg, elevation: 60.64m) and location in Europe; the nearby radiosonde station ID is 26063, its coordinates are latitude: 59.95deg, longitude: 30.70deg, and elevation: 78.00m. The two types of observations are completely in 2017, so we choose the whole year’s observations to do analysis.GNSS can get high-frequency ZTD nearly real-time, which depends on receiver antenna and hardware process ability, which means it has a possibility of getting high-frequency PWV. In this case, the SVTL station provides retrieved PWV 30 min each. In the whole year of 2017, it gets 15,410 observation values; however, in each day, the radiosonde just measures twice, one at 00 : 00 and the other at 12 : 00 a.m., so each day we can get only two PWV values in the 26063 station. Totally, we get 730 PWV values in 2017, so it is easy to find GNSS providing much more PWV observations compared to the radiosonde.In the paper, so as to compare and analyze the retrieved PWV bias and accuracy by the two types of observations, we select each day’s 00 : 00 and 12 : 00 observations from GNSS PWV to do statistical analysis. The following is the distribution map of matched PWV values in each month of 2017.From Figures1 and 2, we can see in the whole year of 2017 GNSS retrieving PWV and radiosonde retrieving PWV; their correlation is up to 0.95. The max PWV of GNSS and RS are both produced between July and September, and the local season is in the summer. Table 1 lists the statistical information min bias, max bias, mean bias, and mean RMSE of matched RS retrieving PWV and GNSS retrieving PWV.Figure 1
The graphic shows SVTL GNSS station and the nearby radiosonde station 26063 PWV time series distributions of the whole year of 2017.Figure 2
PWV correlation analysis of the ground-based GNSS (SVTL) and the matched radiosonde (26063) in the whole year of 2017.Table 1
Precision statistic of retrieved PWV from ground-based GNSS station (STVL) compared with the nearby radiosonde.
Min bias (mm)
Max bias (mm)
Mean bias (mm)
RMSE (mm)
RS_PWV-GPS_PWV
−0.01
15.15
0.23
2.41From Table1, we can easily see that the mean bias is only 0.23 mm and the mean RMSE is 2.41 mm; they are very small, which means GNSS-retrieved PWV values have high precision. It is a useful method to monitor PWV situation changes.


## 3.2. Microwave Satellite Inferred PWV Analysis

We use the day January 1, 2014 (00 : 00), ATOVS (The Advanced TIROS Operational Vertical Sounder) data (Figure3) and radiosonde data to investigate and analyze PWV content distribution status and interpolate ATOVS-detected PWV in 589 radiosonde locations, so we can get the statistical results. In Table 2, the min bias is 0.003 mm, the max bias is 12.068 mm, the mean bias is −0.348 mm, and the mean RMS is 3.461 mm.Figure 3
Distribution of the PWV detected by microwave satellite on January 1, 2014 (00 : 00).Table 2
Precision statistic of retrieved PWV from the ATOVS product compared with radiosondes.
Min bias (mm)
Max bias (mm)
Mean bias (mm)
RMSE (mm)
RS_PWV-Sat_PWV
0.003
12.068
−0.348
3.461


## 3.3. PWV Retrieval from RO and NWP Reanalysis Data

### 3.3.1. Distributions of COSMIC and FY3C RO Events

On February 4, 2016, COSMIC produced 891 radio occultation events and FY3C had 410. Figure4 shows the distribution map of these radio occultation event locations; we can see that these occultation event locations nearly cover the entire world from the North Pole to the South Pole.Figure 4
The two graphics show the locations of COSMIC radio occultation and FY3C radio occultation events produced on February 4, 2016.


### 3.3.2. Retrieved Atmosphere Parameters of FY3C and COSMIC

In Figure5, we select two nearby FY3C and COSMIC radio occultation events as an example for comparison analysis; (a), (b), and (c) are the retrieved specific humidity profile, pressure profile, and temperature profile, respectively, of FY3C radio occultation on the site latitude:24.99°,longitude:−11.68° at 11 : 51 : 03 of the day February 4, 2016, so we find the nearest COSMIC radio occultation event regardless of the location but also the event time which occurred latitude:26.95°,longitude:−11.35° at 14 : 10 : 41 of February 4, 2016; (d), (e), and (f) list the COSMIC-retrieved specific humidity profile, pressure profile, and temperature profile, respectively. From (a) and (d), we can see that radio occultation detected specific humidity whether it is FY3C or COSMIC; in nearly ground height, the value exceeds 4k/kg and the height is between 2 km and 7 km. COSMIC detected that specific humidity has a bigger jump than FY3C, above 10 km, both of which detected a specific humidity reduction to 0, which means no PWV exists in the upper layer. From (b) and (e), the graphics show FY3C and COSMIC retrieving pressure from the ground surface to the upper atmosphere; FY3C shows a pressure of 0~60 km and COSMIC gives 0~40 km. We can see that above 20 km, the air pressure becomes very small; from the ground surface to 20 km, the pressure degrades rapidly. The max value of pressure is located at the ground surface, both showing no more than 1000 Mb; besides, the pressure variation tendencies of both pressure data are in substantial agreement. (c) shows FY3C retrieving temperature from the surface to 60 km; from the graphic, we can see that the highest temperature is in the ground surface, and then from the surface to 20 km, the temperature becomes lower and lower, with the lowest temperature being nearly 20 km in height. However, from 20 km, the temperature (205 K) increases, reaching 255 K in the 42 km height layer. From 42 km, the temperature decreases again until reaching 220 K in the 60 km layer. The same situation occurs in graphic (f); from the ground surface to 20 km, the temperature decreases, especially between the ground surface and the 12 km height layer, in which the temperature variation trend is nearly linearly decreasing; from 20 km to 40 km, the temperature becomes linearly increasing, although FY3C and COSMIC events have more than a 2 h gap, but as can be easily seen, their variation trends are basically the same.Figure 5
The above group graphics show the retrieved atmosphere profiles from FY3C and COSMIC radio occultation observations on February 4, 2016.
(a)
(b)
(c)
(d)
(e)
(f)


### 3.3.3. Evaluation of COSMIC and FY3C Retrieved PWV

The paper uses radiosonde measurements as the reference standard to do analysis; due to this reason, radiosonde observations can only be measured at around 11 : 00 and 23 : 00. We use the following three rules to collect nearby matched RO and RS observations for comparison and do statistical analysis. In the first rule, we set latitude and longitude differences between radiosonde observation sites and radio occultation event which produced sites below1°; for the temporal resolution, we limit it to a ±2h difference. For example, if we use 11 : 00 radiosonde observations, the related radio occultation events must be produced at 9 : 00~13 : 00. That is to say, all of these observations we selected are produced in the time range 9 : 00~13 : 00 and 21 : 00~01 : 00 on February 4, 2016. For the second and third standards, we change the search range from latitude-longitude difference 1°/1° to 2°/2° and 3°/3° individually.On February 4, 2016, COSMIC produced 891 radio occultation events and FY3C had 410. We use the above three standards to collect matched data from FY3C, COSMIC radio occultations, and radiosonde observations in the global area. Statistical results are listed in Table3, whether it is FY3C- or COSMIC-retrieved PWV; compared to nearby radiosonde observations, their RMSE values are no more than 10 mm. From our statistical data, we find that RMSEs at 23 : 00 are all larger than 11 : 00, whether it is FY3C or COSMIC; it is not sure if RO measurements have big biases at 23 : 00 or RS produces large bias. In the future, we will do more deep research about this phenomenon.Table 3
Retrieved PWV precision statistic between FY3C, COSMIC, and radiosonde.
RO
Latitude/longitude difference (deg)
RMSE (mm)
11 : 00
23 : 00
FY3C
1/1
0.34
6.64
2/2
1.09
5.28
3/3
1.26
4.65
COSMIC
1/1
3.03
6.03
2/2
4.30
5.32
3/3
3.29
6.92


### 3.3.4. PWV Retrieved from ECMWF ERA-Interim and ERA5

Figure6 shows the retrieved PWV distribution in the global area from ECMWF ERA-Interim and ERA5 in four periods (00 : 00, 06 : 00, 12 : 00, and 18 : 00); we can see that surrounding the equator, the PWV content exceeds 40 mm, and the PWV distribution decreases along the latitude from the equator to the South Pole to the North Pole.Figure 6
These graphics show the retrieved PWV in global by ECMWF ERA-Interim and ERA5 reanalysis data on February 4, 2016.


### 3.3.5. Evaluation of ERA-Interim and ERA5-Retrieved PWV

In Figure7, it shows the retrieved PWV statistical coefficients between different data sources; we present the statistical analysis using two groups’ observations; respectively, one is in 00 : 00 and the other one uses 12 : 00 observations, so from these graphics, we can see that graph (a) shows the retrieved PWV coefficient between ERA5 and radiosonde at 00 : 00, with R2 equal to 0.9444. From (d), at 12 : 00, R2 equals 0.9708, so we can see that the retrieved PWV by ERA5 and radiosonde has a higher linear correlation. In (b) and (e), they show a retrieved PWV relationship of ERA-Interim and ERA5 at 00 : 00 and 12 : 00, and their R2 values are 0.9833 and 0.9325. Respectively, in theory both methods can get nearly the same results; ERA5 is little superior compared to ERA-Interim, due to the reason that it uses more data sources and new technology for assimilation analysis. (c) and (f) express the relationship of PWV between ERA-Interim and radiosonde, and their R2 values are 0.8699 and 0.9539, respectively, which means they have a higher linear correlation.Figure 7
The PWV correlation statistical of the ERA-Interim, ERA5, and RS. (a, b, c) 00 : 00 observations and (d, e, f) 12 : 00 observations.
(a)
(b)
(c)
(d)
(e)
(f)


### 3.3.6. PWV Intercomparison between RO and ERA-Interim and ERA5

In Figure8, these graphics show the retrieved PWV distribution along the latitude: (a) shows obtained COSMIC-retrieved atmospheric profiles, (b) shows interpolated PWV in the site COSMIC radio occultation events that occurred based on ERA-Interim data, (c) shows interpolation in the same sites from ERA5, (d) shows obtained FY3C radio occultation–retrieved atmospheric profiles, and (e)-(f) show interpolated PWV in FY3C radio occultation events producing sites from ERA-Interim and ERA5 data. From the above graphics, we can easily see that the PWV has symmetrical distribution along the equator; from latitude:0° to the North Pole or South Pole, the retrieved PWV decreases to nearly 0 in the South or North Pole. The biggest PWV is detected in the equator area which exceeds 60g/m2. In (a), (b), and (c), these graphics show that they have nearly the same variation tendency in these sites by different data sources; (d), (e), and (f) have the same situations.Figure 8
The six graphics show the retrieved PWV in global along the latitude by COSMIC radio occultation (a), FY3C radio occultation (d), ERA-Interim (c, f), and ERA5 (b, e) reanalysis data on February 4, 2016. (a, b, c) Retrieved PWV from COSMIC RO and interpolated total column water vapor from ECMWF ERA-Interim and ERA5 reanalysis data on COSMIC radio occultation sites. (d, e, f) Retrieved PWV from FY3C RO and interpolated total column water vapor from ECMWF ERA-Interim and ERA5 reanalysis data on FY3C radio occultation sites.
(a)
(b)
(c)
(d)
(e)
(f)We calculate the mean bias and mean RMSE of COSMIC and FY3C-retrieved PWV with interpolated PWV from ERA-Interim and ERA5 datasets; these results are listed in Table4.Table 4
Retrieved PWV precision statistic between FY3C, COSMIC, and ERA-Interim, ERA5.
Mean bias (mm)
Mean RMSE (mm)
COSMIC RO sites
COSMIC-Interim
−2.3
4.5
COSMIC-ERA5
−2.5
4.8
Interim-ERA5
−0.1
2.0
FY3C RO sites
FY3C-Interim
0.9
3.6
FY3C-ERA5
1.0
3.8
Interim-ERA5
0.1
1.9As can be seen from Table4, the mean bias of COSMIC-retrieved PWV and ERA-Interim is −2.3 mm, RMSE is 4.5 mm, compared to COSMIC-ERA5 results, the statistical value has only −0.2 mm bias, and RMSE has 0.3 mm bias; the same situation occurs in FY3C observations, but on the whole, FY3C-retrieved PWV has smaller mean bias values and mean RMSE than COSMIC-retrieved results. COSMIC has 891 observations, and FY3C contains 410 observations; whether in COSMIC radio occultation sites or FY3C, ERA-Interim and ERA5 have stable statistical results, and the mean bias is around 0 and RMSE is approximately 2.0 mm.In the paper, we also analyze the bias distribution of radio occultation and ECMWF reanalysis data; as can be easily seen in Figure9, the most departures of COSMIC-retrieved PWV and ERA-Interim or ERA5 are no more than 10 cm, and with the same situation to FY3C, most time departures are below 10 mm. For the two data sources of ECMWF, the departures between ERA-Interim and ERA5 mostly are no more than 5 mm.Figure 9
The above graphics show the departures of the retrieved PWV by different observations.
(a)
(b)
(c)
(d)
(e)
(f)


## 3.3.1. Distributions of COSMIC and FY3C RO Events

On February 4, 2016, COSMIC produced 891 radio occultation events and FY3C had 410. Figure4 shows the distribution map of these radio occultation event locations; we can see that these occultation event locations nearly cover the entire world from the North Pole to the South Pole.Figure 4
The two graphics show the locations of COSMIC radio occultation and FY3C radio occultation events produced on February 4, 2016.


## 3.3.2. Retrieved Atmosphere Parameters of FY3C and COSMIC

In Figure5, we select two nearby FY3C and COSMIC radio occultation events as an example for comparison analysis; (a), (b), and (c) are the retrieved specific humidity profile, pressure profile, and temperature profile, respectively, of FY3C radio occultation on the site latitude:24.99°,longitude:−11.68° at 11 : 51 : 03 of the day February 4, 2016, so we find the nearest COSMIC radio occultation event regardless of the location but also the event time which occurred latitude:26.95°,longitude:−11.35° at 14 : 10 : 41 of February 4, 2016; (d), (e), and (f) list the COSMIC-retrieved specific humidity profile, pressure profile, and temperature profile, respectively. From (a) and (d), we can see that radio occultation detected specific humidity whether it is FY3C or COSMIC; in nearly ground height, the value exceeds 4k/kg and the height is between 2 km and 7 km. COSMIC detected that specific humidity has a bigger jump than FY3C, above 10 km, both of which detected a specific humidity reduction to 0, which means no PWV exists in the upper layer. From (b) and (e), the graphics show FY3C and COSMIC retrieving pressure from the ground surface to the upper atmosphere; FY3C shows a pressure of 0~60 km and COSMIC gives 0~40 km. We can see that above 20 km, the air pressure becomes very small; from the ground surface to 20 km, the pressure degrades rapidly. The max value of pressure is located at the ground surface, both showing no more than 1000 Mb; besides, the pressure variation tendencies of both pressure data are in substantial agreement. (c) shows FY3C retrieving temperature from the surface to 60 km; from the graphic, we can see that the highest temperature is in the ground surface, and then from the surface to 20 km, the temperature becomes lower and lower, with the lowest temperature being nearly 20 km in height. However, from 20 km, the temperature (205 K) increases, reaching 255 K in the 42 km height layer. From 42 km, the temperature decreases again until reaching 220 K in the 60 km layer. The same situation occurs in graphic (f); from the ground surface to 20 km, the temperature decreases, especially between the ground surface and the 12 km height layer, in which the temperature variation trend is nearly linearly decreasing; from 20 km to 40 km, the temperature becomes linearly increasing, although FY3C and COSMIC events have more than a 2 h gap, but as can be easily seen, their variation trends are basically the same.Figure 5
The above group graphics show the retrieved atmosphere profiles from FY3C and COSMIC radio occultation observations on February 4, 2016.
(a)
(b)
(c)
(d)
(e)
(f)


## 3.3.3. Evaluation of COSMIC and FY3C Retrieved PWV

The paper uses radiosonde measurements as the reference standard to do analysis; due to this reason, radiosonde observations can only be measured at around 11 : 00 and 23 : 00. We use the following three rules to collect nearby matched RO and RS observations for comparison and do statistical analysis. In the first rule, we set latitude and longitude differences between radiosonde observation sites and radio occultation event which produced sites below1°; for the temporal resolution, we limit it to a ±2h difference. For example, if we use 11 : 00 radiosonde observations, the related radio occultation events must be produced at 9 : 00~13 : 00. That is to say, all of these observations we selected are produced in the time range 9 : 00~13 : 00 and 21 : 00~01 : 00 on February 4, 2016. For the second and third standards, we change the search range from latitude-longitude difference 1°/1° to 2°/2° and 3°/3° individually.On February 4, 2016, COSMIC produced 891 radio occultation events and FY3C had 410. We use the above three standards to collect matched data from FY3C, COSMIC radio occultations, and radiosonde observations in the global area. Statistical results are listed in Table3, whether it is FY3C- or COSMIC-retrieved PWV; compared to nearby radiosonde observations, their RMSE values are no more than 10 mm. From our statistical data, we find that RMSEs at 23 : 00 are all larger than 11 : 00, whether it is FY3C or COSMIC; it is not sure if RO measurements have big biases at 23 : 00 or RS produces large bias. In the future, we will do more deep research about this phenomenon.Table 3
Retrieved PWV precision statistic between FY3C, COSMIC, and radiosonde.
RO
Latitude/longitude difference (deg)
RMSE (mm)
11 : 00
23 : 00
FY3C
1/1
0.34
6.64
2/2
1.09
5.28
3/3
1.26
4.65
COSMIC
1/1
3.03
6.03
2/2
4.30
5.32
3/3
3.29
6.92


## 3.3.4. PWV Retrieved from ECMWF ERA-Interim and ERA5

Figure6 shows the retrieved PWV distribution in the global area from ECMWF ERA-Interim and ERA5 in four periods (00 : 00, 06 : 00, 12 : 00, and 18 : 00); we can see that surrounding the equator, the PWV content exceeds 40 mm, and the PWV distribution decreases along the latitude from the equator to the South Pole to the North Pole.Figure 6
These graphics show the retrieved PWV in global by ECMWF ERA-Interim and ERA5 reanalysis data on February 4, 2016.


## 3.3.5. Evaluation of ERA-Interim and ERA5-Retrieved PWV

In Figure7, it shows the retrieved PWV statistical coefficients between different data sources; we present the statistical analysis using two groups’ observations; respectively, one is in 00 : 00 and the other one uses 12 : 00 observations, so from these graphics, we can see that graph (a) shows the retrieved PWV coefficient between ERA5 and radiosonde at 00 : 00, with R2 equal to 0.9444. From (d), at 12 : 00, R2 equals 0.9708, so we can see that the retrieved PWV by ERA5 and radiosonde has a higher linear correlation. In (b) and (e), they show a retrieved PWV relationship of ERA-Interim and ERA5 at 00 : 00 and 12 : 00, and their R2 values are 0.9833 and 0.9325. Respectively, in theory both methods can get nearly the same results; ERA5 is little superior compared to ERA-Interim, due to the reason that it uses more data sources and new technology for assimilation analysis. (c) and (f) express the relationship of PWV between ERA-Interim and radiosonde, and their R2 values are 0.8699 and 0.9539, respectively, which means they have a higher linear correlation.Figure 7
The PWV correlation statistical of the ERA-Interim, ERA5, and RS. (a, b, c) 00 : 00 observations and (d, e, f) 12 : 00 observations.
(a)
(b)
(c)
(d)
(e)
(f)


## 3.3.6. PWV Intercomparison between RO and ERA-Interim and ERA5

In Figure8, these graphics show the retrieved PWV distribution along the latitude: (a) shows obtained COSMIC-retrieved atmospheric profiles, (b) shows interpolated PWV in the site COSMIC radio occultation events that occurred based on ERA-Interim data, (c) shows interpolation in the same sites from ERA5, (d) shows obtained FY3C radio occultation–retrieved atmospheric profiles, and (e)-(f) show interpolated PWV in FY3C radio occultation events producing sites from ERA-Interim and ERA5 data. From the above graphics, we can easily see that the PWV has symmetrical distribution along the equator; from latitude:0° to the North Pole or South Pole, the retrieved PWV decreases to nearly 0 in the South or North Pole. The biggest PWV is detected in the equator area which exceeds 60g/m2. In (a), (b), and (c), these graphics show that they have nearly the same variation tendency in these sites by different data sources; (d), (e), and (f) have the same situations.Figure 8
The six graphics show the retrieved PWV in global along the latitude by COSMIC radio occultation (a), FY3C radio occultation (d), ERA-Interim (c, f), and ERA5 (b, e) reanalysis data on February 4, 2016. (a, b, c) Retrieved PWV from COSMIC RO and interpolated total column water vapor from ECMWF ERA-Interim and ERA5 reanalysis data on COSMIC radio occultation sites. (d, e, f) Retrieved PWV from FY3C RO and interpolated total column water vapor from ECMWF ERA-Interim and ERA5 reanalysis data on FY3C radio occultation sites.
(a)
(b)
(c)
(d)
(e)
(f)We calculate the mean bias and mean RMSE of COSMIC and FY3C-retrieved PWV with interpolated PWV from ERA-Interim and ERA5 datasets; these results are listed in Table4.Table 4
Retrieved PWV precision statistic between FY3C, COSMIC, and ERA-Interim, ERA5.
Mean bias (mm)
Mean RMSE (mm)
COSMIC RO sites
COSMIC-Interim
−2.3
4.5
COSMIC-ERA5
−2.5
4.8
Interim-ERA5
−0.1
2.0
FY3C RO sites
FY3C-Interim
0.9
3.6
FY3C-ERA5
1.0
3.8
Interim-ERA5
0.1
1.9As can be seen from Table4, the mean bias of COSMIC-retrieved PWV and ERA-Interim is −2.3 mm, RMSE is 4.5 mm, compared to COSMIC-ERA5 results, the statistical value has only −0.2 mm bias, and RMSE has 0.3 mm bias; the same situation occurs in FY3C observations, but on the whole, FY3C-retrieved PWV has smaller mean bias values and mean RMSE than COSMIC-retrieved results. COSMIC has 891 observations, and FY3C contains 410 observations; whether in COSMIC radio occultation sites or FY3C, ERA-Interim and ERA5 have stable statistical results, and the mean bias is around 0 and RMSE is approximately 2.0 mm.In the paper, we also analyze the bias distribution of radio occultation and ECMWF reanalysis data; as can be easily seen in Figure9, the most departures of COSMIC-retrieved PWV and ERA-Interim or ERA5 are no more than 10 cm, and with the same situation to FY3C, most time departures are below 10 mm. For the two data sources of ECMWF, the departures between ERA-Interim and ERA5 mostly are no more than 5 mm.Figure 9
The above graphics show the departures of the retrieved PWV by different observations.
(a)
(b)
(c)
(d)
(e)
(f)


## 4. Conclusions

In this paper, generally, we analyze and compare retrieved PWV from ground-based GNSS and RS and space-based observation methods, such as microwave satellite, radio occultation and ECMWF reanalysis data, ERA-Interim, and ERA5. This section will show more detailed conclusions:(1)
We analyze the ground-based GNSS retrieved PWV and compare its results with radiosonde; GNSS can not only get high-frequency PWV values, but still they are of high precision, and the mean RMSE is about 2.41 mm, so it is real useful and meaningful to monitor PWV distribution and change situations by the ground-based GNSS technique
(2)
We analyze the PWV from microwave satellite ATOVS and interpolate it to global 589 radiosonde sites, compare the detected PWV of ATOVS with radiosonde, and get a statistical min bias of 0.003 mm, a max bias of 12.068 mm, a mean bias of −0.348 mm, and a mean RMS of 3.461 mm. It can satisfy the PWV detection requirement and should be a valuable method to get PWV distribution from space
(3)
We verify China FY3C radio occultation observations and compare its retrieved atmospheric parameters with COSMIC; we find that there is no big bias regardless of whether these are FY3C-retrieved atmospheric parameters, such as temperature, pressure, or water vapor pressure profiles. Besides, although FY3C can provide retrieved atmosphere observations above 40 km, for the retrieved PWV of FY3C, we compare it to COSMIC’s results and find that on the whole it has a smaller bias, compared to the reference radiosondes, and their average RMSE values are both no more than 10 mm. From our statistic results, whether it is FY3C or COSMIC, on the whole, the RMSE value at 11 : 00 (day) is smaller than 23 : 00’s results (night). We need do more future research to understand the reasons in detail
(4)
We also analyze the NWP reanalysis data with a traditional radiosonde, interpolate the PWV from ERA-Interim and new reanalysis dataset ERA5 to radiosonde sites, and find that on the whole they have high consistency and no big biases. However, relatively, for the PWV of ERA5, we find that it has higher precision compared to ERA-Interim data, and it is also reasonable, because ERA5 assimilates more different types of observations and has higher spatial and temporal resolution
(5)
We try to detect the retrieved PWV bias between RO and ECMWF reanalysis data, so we compare FY3C and COSMIC radio occultation results with ERA-Interim and ERA5. From the precision statistical results and correlation analysis, we know that the mean bias is 0.9 mm and the mean RMSE is 3.6 mm between FY3C results and ERA-Interim; comparing FY3C results to those of ERA5, the mean bias and mean RMSE are 1.0 mm, 3.8 mm, respectively; however, the mean bias (around 2.5 mm) and mean RMSE (around 4.5 mm) between COSMIC and ECMWF reanalysis data are both larger than those of FY3C results. But on the whole, regardless of whether it is COSMIC, FY3C, ERA-Interim, or ERA5, their PWV all have a higher correlation



---

*Source: 3428303-2018-12-24.xml*