## 2 Data and Method  
Radiosonde data used in this paper are from the National Climatic Data Center (NCDC, Asheville, NC, USA) Integrated Global Radiosonde Data (IGRA) version 2 [20]. IGRA consists of radiosonde and pilot balloon observations over 2000 globally distributed stations and over 60 years (~1957-2016) of records. The records include various meteorological elements such as temperature, relative humidity, pressure, wind speed, wind direction at standard levels (including the surface level and 1000, 925, 825, 700, 500 hPa, etc.) and at significant levels. Satellite-based global positioning system radio occultation soundings from the FORMOSAT-3/COSMIC (Formosa Satellite Mission 3/Constellation Observing System for Meteorology, Ionosphere, and Climate) have long-term, all-weather, global operational measurements of water vapor profiles, which provided additional observations over the oceans and the Southern Hemisphere for the most recent ten years (2007-2016) [2]. The FORMOSAT-3/COSMIC mission is a collaborative project held by the University Corporation for Atmospheric Research in the USA and the National Space Organization (NSPO) in Taiwan. The GPS radio-occultation receivers on the COSMIC satellites measure the phase delay of radio waves from satellites when occulted by the Earth's atmosphere. Profiles of atmospheric refractivity can be derived from the bending angles of radio wave trajectories in the ionosphere, stratosphere, and troposphere. The GPS RO retrieved wet profiles from the COSMIC Data Analysis and Archive Center consisting of temperature, pressure and analyzed refractivity up to 40 km in each profile and containing over 1000 daily profiles distributed globally with high vertical resolution (100 m). The raw data from the RO measurements can be processed to derive bending angles, which can be converted to the atmospheric refractivity (\(N\)) under the spherical symmetry assumption [21]. In the upper and middle troposphere, the refractivity is mainly used to obtain temperature. In the lower troposphere, the information of refractivity is mainly used to retrieve water vapor, since the atmospheric refractivity is more sensitive to water vapor variability than to the temperature [22]. Following [PERSON] et al. [19], the water vapor pressure can be calculated by the relationship among \(N\), temperature (\(T\)) in Kelvin, air pressure (\(P\)) in hPa and water vapor pressure (\(e\)) in hPa:  
\[N=77.6\frac{P}{T}+3.73\times 10^{5}\frac{e}{T^{2}} \tag{1}\]  
The specific humidity is easily calculated by water vapor pressure and air pressure using the Equation (2) [23]. (\(\epsilon=0.622\))  
\[q=\frac{\epsilon e}{p-(1-\epsilon)e} \tag{2}\]To calculate PW value of one single wet profile, the specific humidity of each height is integrated from the surface to the retrieved maximum height, using the Equation (3). In the Equation (3), \(W\) is the PW value in mm, \(\rho\) is the density of liquid water, \(P\) is air pressure, \(P_{s}\) is the surface air pressure, \(g\) is the gravitational acceleration value. \[w=\frac{1}{\rho g}\int_{0}^{P_{s}}qdp=-\frac{1}{g}\sum_{P_{s}}^{0}q\Delta p \tag{3}\]  
The existence of super-refraction in the lower troposphere might lead to negative refractivity [12], which might cause the biases of temperature and humidity and void data from the surface to the lowest perigee height. In this study, the void data of layers is compensated with the same value of the lowest perigee and, to reduce the uncertainty, profiles with a lowest perigee height higher than 1.0 km are discarded. It is noteworthy that, since the lowest ray perigee heights often happen in the drier atmosphere, some sampling errors may be introduced. However, [PERSON] et al. [19] proved that this bias in PW will be less than 2 mm on average. The global PW values are then averaged into regular grid points with the resolution of \(5^{\circ}\) latitude by \(5^{\circ}\) longitude bins.