Metrology Assessment of the Accuracy of Precipitable Water Vapor Estimates from GPS Data Acquisition in Tropical Areas: The Tahiti Case

High precision Global Positioning System (GPS) receivers, with the advantages of all-weather work and low cost, are now widely used to routinely monitor precipitable water (PW) vapor. They are so successful that the progressive phasing out of the costly and sparse in situ radio soundings (RS) is now a certainty. Nevertheless, the sub-daily to annual monitoring of high levels of the PW by GPS receivers in the tropics and the equatorial area still needs to be asserted in terms of metrology accuracy. This is the subject of this paper, which focuses on a tropical site located in mid-ocean (Tahiti). The metrology assessment was divided into two steps. Firstly, a GPS internal assessment, with an in-house processing based on the Bernese GNSS Software Version 5.2 and a comparison with the Center for Orbit Determination in Europe (CODE) products. Secondly, an external assessment, with a comparison with RS PW estimates. In contrast with previous works that only used PW estimates from the Integrated Global Radiosonde Archive (IGRA) website, we estimated the RS PW from the balloon raw data. This is especially important in tropical areas, where IGRA estimates only consider balloon measurements taken below approximately 5500 m. We show that, in our case, this threshold is one of the main sources of bias between GPS and RS estimates, and that the formula used to translate the GPS zenith wet delays (ZWD) to PW estimates also needs to be revisited for high level water vapor contents in the atmosphere.


Introduction
Atmospheric water vapor plays an important role in atmospheric processes including atmospheric radiation balance, water cycles, the transfer of energy, and the formation of clouds and precipitation [1][2][3]. Water vapor is also the most abundant and a critical component of the greenhouse gases driving global weather and climate changes [4][5][6]. Integrated precipitable water vapor (IPWV, or PW for short) is the total water vapor contained in an air column from the Earth's surface to the top of the atmosphere. About 45-65% of the PW is included from the surface to 850 hPa altitude (roughly 1500 m) and its distribution is highly variable in time and space [7,8]. Real-time and accurate estimates of characterization of the PW at high spatial and temporal resolution are now needed for severe weather forecasting and climate warming studies [9,10].
Many meteorological techniques that have been used for decades to retrieve PW, such as radiosondes, water vapor radiometer, and infrared sounders are expensive and may give poor data the Bevis T m model take place in Section 5. In Section 6, we compare and analyze the GPS-derived PW estimates based on our new T m models with the corresponding RS PW estimates. Conclusions are presented in Section 7.

The Study Area
The study area of this paper is the Tahiti island of French Polynesia, located on the tropical South Central Pacific Ocean (149 • 25 W, 17 • 40 S; Figure 1 [44]). The island of Tahiti (around 1000 km 2 ) includes two volcanic systems that culminate at 2241 m (Tahiti-Nui), and at 1332 m in the southeast peninsula (Tahiti-Iti) [45]. The local climate, much more humid than in the other islands of French Polynesia [46], is governed by the South Pacific Convergence Zone (SPCZ) on the large scale, and by the topography of a high volcanic island on the orographic scale, including a fresh and dry period ("dry" season, with August being the driest month with 48 mm of precipitations) when southeast winds prevail (from May to October) and a warm and rainy period ("wet" season, the wettest month of which is January, peaking at 340 mm of rain in average) when northeast trade winds are predominant (from November to April) [47]. During the wet season, from November to April, significant land-slides and flash-floods are common, generating large property damages and casualties. For our data processing, we used the GPS data from the permanent IGS station (THTI) of the OGT and the radio soundings from the nearby (2 km) Météo-France weather station. Both locations ( Figure 1) are in the suburbs of Tahiti capital city, Papeete.

GPS and Weather Data
The permanent GPS station THTI (IGS reference name) of the OGT (about 100 m altitude) consists of a TRIMBLE NETR8 receiver with an ASHTECH geodetic L1/L2 antenna. GPS satellite precise orbits and clocks as well as consistent Earth-rotation parameters provided by the Center for Orbit Determination in Europe (CODE), together with precise positioning point (PPP) approach under the Bernese GNSS Software Version 2, were used to estimate ZTD [43]. The software parameters are set, in all our processing, as follows: sampling rate is 300 s, elevation cut-off angle is 3 degrees, temporal resolution is 2 h and ionosphere correction is L3 free linear combination. We also considered the atmospheric loading and the ocean tidal loading corrections in the Bernese processing. In this study, we used the hydrostatic Vienna mapping function (VMF1) [48], with input data from the European Center for Medium-Range Weather Forecasts (ECMWF), as the a priori ZHD model, with wet and dry VMF1 mapping function [49,50].
The GPS-derived ZTD comprises ZHD and ZWD, as shown in Equation (1): And, the relationship between ZWD and PW [27] can be expressed as: where Π = 10 6 ρ·R v ( k 3 T m + k 2 ) and where ρ = 1000 kg/m 3 is the density of liquid water, R v = 461.495 J/Kg·K is the specific gas constant of water vapor, k 2 = 22.1 K/mb and k 3 = 3.739·10 5 K 2 /mb are physical constants. The coefficient Π is a conversion coefficient, in practice a mean value of about 0.15 ± 20%. P v is the partial pressure of water vapor, T is the absolute temperature in Kelvin. Both P v and T are pointwise values taken along the vertical column. In our case, T m (in Kelvin) is linearly related to T s , as given by Equation (5): where a and b are the intercept and slope of this linear function. The estimation of ZWD (Equation (1)) from GPS-derived ZTD essentially requires empirical models of the ZHD, like the Saastamoinen model, the accuracy of which is approximately 1-2 mm [51] for the ZHD estimation.
The Saastamoinen model is defined as: where f (λ, H) = 1 − 0.00266· cos(2λ) − 0.00028·H, λ (rad) is the station latitude, H (km) is the geoid height in kilometers and P s (hPa) is the surface pressure. For sub-diurnal studies of the water vapor contents from GNSS data, the sampling rate is usually 30 min or less, but in our case, as we were primarily interested in comparisons with the CODE products and radio soundings in regard to a seasonal time frame, a 2 h sampling is perfectly adapted. where 1000 kg/m is the density of liquid water, 461.495 J/Kg • K is the specific gas constant of water vapor, 22.1 K/mb and 3.739 • 10 K /mb are physical constants. The coefficient Π is a conversion coefficient, in practice a mean value of about 0.15 ± 20%. Pv is the partial pressure of water vapor, T is the absolute temperature in Kelvin. Both Pv and T are pointwise values taken along the vertical column. In our case, Tm (in Kelvin) is linearly related to Ts, as given by Equation (5): where a and b are the intercept and slope of this linear function.
The estimation of ZWD (Equation (1)) from GPS-derived ZTD essentially requires empirical models of the ZHD, like the Saastamoinen model, the accuracy of which is approximately 1-2 mm [51] for the ZHD estimation.
The Saastamoinen model is defined as: For sub-diurnal studies of the water vapor contents from GNSS data, the sampling rate is usually 30 min or less, but in our case, as we were primarily interested in comparisons with the CODE products and radio soundings in regard to a seasonal time frame, a 2 h sampling is perfectly adapted. In our study, the meteorological parameters (the surface pressure Ps and the surface temperature Ts) are from the gridded VMF1 (vmf1_g) data of the European Center for Medium Range Weather Forecasts (ECMWF, http://ggosatm.hg.tuwien.ac.at). They are based on a 2° 2.5°global grid (2 degrees sampling from North to South and 2.5 degrees sampling from West to In our study, the meteorological parameters (the surface pressure P s and the surface temperature T s ) are from the gridded VMF1 (vmf1_g) data of the European Center for Medium Range Weather Forecasts (ECMWF, http://ggosatm.hg.tuwien.ac.at). They are based on a 2 • × 2.5 • global grid (2 degrees sampling from North to South and 2.5 degrees sampling from West to East), with a 6-hourly temporal resolution (at 0h00, 6h00, 12h00 and 18h00 UTC). These gridded data have facilitated some high accuracy geodetic applications over the past years [52]. To validate the vmf1_g, we compared the  vmf1_g downloaded values to the corresponding values measured at the OGT site (weather station of  the TAH1 GPS receiver, ftp://garner.ucsd.edu, collocated within a few meters from THTI).
To validate the accuracy of the estimates of the ZTD and ZWD computed through the Bernese software with the VMF1 mapping function, we compared them with the CODE troposphere products (ftp://ftp.aiub.unibe.ch) and with the ZTD and ZWD derived from the Saastamoinen model (SAAS ZWD) relative to the THTI station with on-site weather data. The CODE products are also calculated by using the Bernese 5.2 software, but with a relative positioning technique and ionosphere-free linear combination of L1 and L2, their sampling rate is 180 s. The a priori ZHD model and mapping functions and the other processing parameters of CODE troposphere products are, from their documentation [43], the same as the ones used for our GPS data processing.

Radiosonde Level 1 and Level 2 Data
As mentioned above, a radio sounding balloon station is located close to our THTI GPS receiver in Tahiti Island ( Figure 1). The horizontal distance is about 2 km and the height difference is around 100 m [44]. Balloons are launched once a day (0:00 UT) or twice a day (0:00 and 23:00 UT or 0:00 and 12:00 UT). They include pressure, temperature, relative humidity velocity and wind direction sensors.
In our study, the RS PW values coming from this station are relative to two different processing levels of the same raw data. The level 1 PW values are retrieved from the Integrated Global Radiosonde Archive (IGRA: ftp://ftp.ncdc.noaa.gov/) managed by the National Oceanic and Atmospheric Administration (NOAA) [53,54]. This archive is highly compressed, based on the characteristic points of data profiles, from the surface to 500 hPa (about 5500 m). This geopotential altitude is often referred to as a steering level, because the weather systems beneath, near to the Earth's surface, roughly move in the same direction as the winds at the 500 hPa level. We computed in our observatory the level 2 data from the 1-second balloon raw data that we retrieved from the Météo-France in-house archive. They consist of pressure, temperature, relative humidity and GPS positioning from the surface to the maximum altitude of the balloon, i.e., 25,500 m (around 1 h 20 min of ascent). When the balloon explodes, its lateral drift can reach up to 20 to 50 km, depending on the direction and intensity of the wind. Macpherson [55] found that we can neglect the space-time drift of radiosonde balloons during their ascent for mesoscale numerical weather prediction models. Specifically, while the impact of balloon drift can indeed be detected, it is small compared with the impact from the observation when assimilated into the numerical weather prediction model with neglect of balloon drift. We computed these level 2 PW values by applying the procedure outlined in [8,56] as: with and e s = f ·6.1121·e (18.729− t 227.3 )·t t+257.87 (9) where h is the geopotential altitude, ρ is the density of liquid water, ρ w is the density of water vapor, RH is the relative humidity, R v is the specific gas constant of water vapor, T = t + 273.15, t is the surface temperature in Celsius, e s is the saturation vapor pressure in hPa, f = 1.0007 + 3.46·P s ·10 −6 is the enhancement factor, and P s is the surface water vapor pressure.

Construction of the T m Models
As we have two series of GPS (ZWD_GPS) and RS (PW_RS) data coming side by side, we can construct time series of "instantaneous" Π according to Equation (10), which is the reversed version of Equation (2): Then the time series of Π values can be linearly fitted with respect to the ground temperature T s for a given duration of time (Equation (5)). We can then distinguish many subcases, depending on how PW_RS is exactly computed (choice of mapping function and local tailoring of the chosen mapping function, including temperature and pressure) and which level (1 or 2) is chosen for PW_RS. We can also distinguish between wet and dry seasons for the fit to the local ground temperature T s .
For the meteorological parameters entering the mapping functions, we used the data from the VMF1/ECWMF archive and local data from the nearby local station GPS station TAH1 (which also provides the T s time series).

Definitions
In the following sections, we will use the usual statistical tools of Bias, root mean square (RMS) and standard deviation (STD) between two time series. To avoid any misunderstanding, we precise the definitions here: where i = 1, 2, 3, . . . , n. is the number of data point, S and R represent any time series with the same time sampling, with R taken as a reference. We will also do time averaging, in the form [45]: where S[j] is the value of S corresponding to the j-th epoch. We will always use an averaging over one month, t = MJD (Modified Julian Day), ∆t = t j − t i , i, j = 1, 2, 3, . . . , n, if |∆t| > T, we set v ij = 0, otherwise, v ij = 1.

Comparisons of Our GPS-Derived ZTD and ZWD Values with CODE Products
In a first step, we compared home-brew ZTD and ZWD estimates based on our Local VMF1 Processing (LVP) over three years (2011)(2012)(2013) for the THTI station with the estimates from the CODE products relative to the same station, to detect any gross differences. For this purpose, we employed the Bernese 5.2 software with the hydrostatic VMF1/ECMWF model (VMF1 dry and wet mapping function). Nowadays, the VMF1 mapping function is widely seen as the mapping function providing, from a global point of view, the most accurate and reliable geodetic results [50]. It should be noted that we used the PPP approach, but that CODE products, also processed with the help of the Bernese 5.2 software and the VMF1 dry and wet mapping function, are using a relative positioning approach and an ionosphere-free linear combination of L1 and L2, with a sampling rate of 180 s. Figure 2 and Table 1 summarize the results. The distribution of the differences between our LVP processing and the CODE products, both for the ZTD and ZWD estimates are clearly close to a Gaussian distribution with no appreciable bias and an RMS of 6 mm. We can then consider the two processing as consistent. The ZWD differences are probably coming from the choice of PPP versus relative positioning, and the different ionosphere-free linear combination and sampling rate settings. The ZTD differences are caused by the ZWD differences. providing, from a global point of view, the most accurate and reliable geodetic results [50]. It should be noted that we used the PPP approach, but that CODE products, also processed with the help of the Bernese 5.2 software and the VMF1 dry and wet mapping function, are using a relative positioning approach and an ionosphere-free linear combination of L1 and L2, with a sampling rate of 180 s. Figure 2 and Table 1 summarize the results. The distribution of the differences between our LVP processing and the CODE products, both for the ZTD and ZWD estimates are clearly close to a Gaussian distribution with no appreciable bias and an RMS of 6 mm. We can then consider the two processing as consistent. The ZWD differences are probably coming from the choice of PPP versus relative positioning, and the different ionosphere-free linear combination and sampling rate settings. The ZTD differences are caused by the ZWD differences. (ZTD) from our GPS data processing, our Local VMF1 Processing (LVP) (cyan dots) and Center for Orbit Determination in Europe (CODE) products (magenta dots), and monthly averaged estimates of the ZTD from our GPS data processing (blue dots) and CODE products (red dots), and (c) estimates of the zenith wet delays (ZWD) from our GPS data processing (cyan dots) and CODE products (magenta dots), and monthly averaged estimates of the ZWD from our GPS data processing (blue dots) and CODE products (red dots). Right column: Histogram of (b) the raw ZTD differences between CODE ZTD and LVP ZTD, and (d) the raw ZWD differences between CODE ZWD and LVP ZWD.  from our GPS data processing, our Local VMF1 Processing (LVP) (cyan dots) and Center for Orbit Determination in Europe (CODE) products (magenta dots), and monthly averaged estimates of the ZTD from our GPS data processing (blue dots) and CODE products (red dots), and (c) estimates of the zenith wet delays (ZWD) from our GPS data processing (cyan dots) and CODE products (magenta dots), and monthly averaged estimates of the ZWD from our GPS data processing (blue dots) and CODE products (red dots). Right column: Histogram of (b) the raw ZTD differences between CODE ZTD and LVP ZTD, and (d) the raw ZWD differences between CODE ZWD and LVP ZWD. Besides, we noted that some El Niño events took place from 2014 to 2016 [46]. To validate our LVP ZTD and ZWD, we also made comparisons between our results and CODE products during these three years. Table 1 shows good agreement between our LVP ZTD and ZWD and CODE products. The bias between the ZTD and ZWD time series is respectively 0.34 mm and 0.50 mm and RMS and STD are all about 7 mm. The above discussions reflect that the results of our LVP ZTD and ZWD meet the accuracy requirement and can be used in the following data processing.

Comparison of Our SAAS ZWD with Our LVP ZWD and CODE ZWD
We did an estimation of the ZWD (called here SAAS ZWD) over the period of 2014-2016 by using the LVP ZTD from the Bernese 5.2 data processing and the Saastamoinen model [51], fed with pressure data from the vmf1_g, because the weather station of the TAH1 GPS receiver stopped working in 2015. Figure 3a,b and Table 2 summarize the results with a comparison with the ZWD estimates from the previous section (VMF1 mapping functions) and CODE products. Besides, we noted that some El Niño events took place from 2014 to 2016 [46]. To validate our LVP ZTD and ZWD, we also made comparisons between our results and CODE products during these three years. Table 1 shows good agreement between our LVP ZTD and ZWD and CODE products. The bias between the ZTD and ZWD time series is respectively 0.34 mm and 0.50 mm and RMS and STD are all about 7 mm. The above discussions reflect that the results of our LVP ZTD and ZWD meet the accuracy requirement and can be used in the following data processing.

Comparison of Our SAAS ZWD with Our LVP ZWD and CODE ZWD
We did an estimation of the ZWD (called here SAAS ZWD) over the period of 2014-2016 by using the LVP ZTD from the Bernese 5.2 data processing and the Saastamoinen model [51], fed with pressure data from the vmf1_g, because the weather station of the TAH1 GPS receiver stopped working in 2015. Figure 3a,b and Table 2 summarize the results with a comparison with the ZWD estimates from the previous section (VMF1 mapping functions) and CODE products.   Table 2. Statistical summary for the comparisons between our LVP ZWD estimates and SAAS_1 ZWD estimates (LVP minus SAAS_1), CODE ZWD and SAAS_2 ZWD estimates (CODE minus SAAS_2) in terms of max, min, bias, STD and RMS, from 2014 to 2016, relative to Figure 3a,b.  Tables 1 and 2 indicate that all the processing (in-house VMF1 and on Saastamoinen) give essentially the same results (internal consistency), paving the way to a comparison of our LVP ZWD and their derived PW values with the RS PW values.

Comparison of Level 1 RS PW Values with Level 2 RS PW Values
We will now assert the internal consistency of RS data. For this purpose, we recalculate in this section the local PW data from raw balloon measurements (level 2 RS PW in the jargon of the meteorological community, that we split into two subsets: a/level 2-A RS PW values, from the surface to the 500 hPa level and b/level 2-B RS PW values from the surface to the maximum altitude of the balloon, roughly 25,500 m in the mean) and make a comparison with the corresponding    Tables 1 and 2 indicate that all the processing (in-house VMF1 and on Saastamoinen) give essentially the same results (internal consistency), paving the way to a comparison of our LVP ZWD and their derived PW values with the RS PW values.

Comparison of Level 1 RS PW Values with Level 2 RS PW Values
We will now assert the internal consistency of RS data. For this purpose, we recalculate in this section the local PW data from raw balloon measurements (level 2 RS PW in the jargon of the meteorological community, that we split into two subsets: a/level 2-A RS PW values, from the surface to the 500 hPa level and b/level 2-B RS PW values from the surface to the maximum altitude of the balloon, roughly 25,500 m in the mean) and make a comparison with the corresponding archived PW data from the IGRA website (level 1 RS PW). We emphasize that the IGRA archive contains only the level 1 data, that are highly compressed by the archive builders, which only retained the "characteristics" points of the raw data profiles, in the sense that the behavior of the data is identified as linear between these points. Besides, the archive builders limited, by construction, the maximum altitude of the data used to compute the archive from the surface to the 500 hPa level. In other words, they assumed that no PW was present beyond this level, or that this PW is irrelevant for meteorological studies. Still, in other words, the level 1 data should be, at least, consistent with the level 2-A if the compression algorithm is correct, and consistent with level 2-B if the assumption that no PW is present beyond the 500 hPa level is correct. This is the topic of the following sub-sections.

Comparison of Level 1 RS PW Values with Level 2-A RS PW Values (from the Surface to the 500 hPa Level)
In this section, we analyze and compare the reconstructed and archived PW from the surface to 500 hPa altitude with the level 2-A RS PW values we computed from the balloon raw data. We only used the 2014-2016 period, as instrumentation was changed at the end of 2013 at the Météo-France weather station. We recall that level 1 RS PW values are relative to the surface to the 500 hPa level and were downloaded from the IGRA archive, and that level 2-B RS PW values are relative to the surface to the maximum altitude of the balloon, far over the tropopause. We recomputed the level 2 data from raw balloon data stored in Météo-France archive, but here only up to the 500 hPa level (so-called level 2-A) to assert the reliability of archived level 1 PW data with regard to this threshold altitude.  Figure 4b shows the relationship between level 1 RS PW and level 2-A RS PW, and their least-squares fit is: level_1 = 0.95 × level_2-A − 2.29, with R 2 = 94.62%. Figure 4a and Table 3 show clearly that the two sets of values are consistent, with a small bias (−0.10 mm) and an RMS of 1.30 mm. this means that the compression algorithm used to save space on the IGRA archive (storage of characteristic points of data profiles) is working, but with only a 95% "performance" in our particular case.
Remote Sens. 2018, 10, x FOR PEER REVIEW 9 of 22 archived PW data from the IGRA website (level 1 RS PW). We emphasize that the IGRA archive contains only the level 1 data, that are highly compressed by the archive builders, which only retained the "characteristics" points of the raw data profiles, in the sense that the behavior of the data is identified as linear between these points. Besides, the archive builders limited, by construction, the maximum altitude of the data used to compute the archive from the surface to the 500 hPa level. In other words, they assumed that no PW was present beyond this level, or that this PW is irrelevant for meteorological studies. Still, in other words, the level 1 data should be, at least, consistent with the level 2-A if the compression algorithm is correct, and consistent with level 2-B if the assumption that no PW is present beyond the 500 hPa level is correct. This is the topic of the following sub-sections.

Comparison of Level 1 RS PW Values with Level 2-A RS PW Values (from the Surface to the 500 hPa Level)
In this section, we analyze and compare the reconstructed and archived PW from the surface to 500 hPa altitude with the level 2-A RS PW values we computed from the balloon raw data. We only used the 2014-2016 period, as instrumentation was changed at the end of 2013 at the Météo-France weather station. We recall that level 1 RS PW values are relative to the surface to the 500 hPa level and were downloaded from the IGRA archive, and that level 2-B RS PW values are relative to the surface to the maximum altitude of the balloon, far over the tropopause. We recomputed the level 2 data from raw balloon data stored in Météo-France archive, but here only up to the 500 hPa level (so-called level 2-A) to assert the reliability of archived level 1 PW data with regard to this threshold altitude.  Figure 4b shows the relationship between level 1 RS PW and level 2-A RS PW, and their least-squares fit is: level_1 = 0.95 × level_2-A − 2.29, with R 2 = 94.62%. Figure 4a and Table 3 show clearly that the two sets of values are consistent, with a small bias (−0.10 mm) and an RMS of 1.30 mm. this means that the compression algorithm used to save space on the IGRA archive (storage of characteristic points of data profiles) is working, but with only a 95% "performance" in our particular case.    Now, we compare the same sets of level 2-A data taken from the ground to 5500 m with regard to level 2-B data, i.e., data taken from the ground to the maximum altitude of the balloon (over 25,500 m).  Table 4. The bias is 2.75 mm, the STD is 2.32 mm and the RMS is 3.60 mm. In order to better understand the relationship between the balloon altitude and the PW, we map in Figure 5b the level 2 data from a characteristic balloon launch at the Météo-France site, close to our observatory. The balloon reaches its maximum altitude, quasi-linearly with time, in about one hour. It is also clear from Figure 5b that the PW contents, integrated from the surface to the current altitude of the balloon also reaches a horizontal asymptote, and that this asymptote is pretty close, but does not correspond, to the 500 hPa level. About 95% of the maximum value of the PW is acquired after 15 min of ascent, and the last 5% in the remaining part of the ascent. Figure 5c shows the relationship between level 2-A RS PW values and level 2-B RS PW values. Their least-squares linear fit is: level_2_B = 1.09 × level_2-A − 1.35, with R 2 = 94.46%. Figure 5d shows the relationship between level 1 RS PW values and level 2-A RS PW values, and their least-squares linear fit is: level_2-B = 1.14 × level_1 − 3.29, R 2 = 98.46%. Figure 5c and its linear fit indicates that, in the mean, the level 2-A values, in our case, underestimate the PW by about 9% with respect to the level 2-B values. One can argue that this difference is caused by the balloon lateral drift that can reach tens of kilometers when the balloon bursts. But the direction of the drift is highly variable as a function of the wind conditions during the balloon ascent, and this drift is essentially the cause of the scatter of the data points (about 2000 balloon launches) in Figure 5c. Therefore, our conclusion is still valid, in a certain sense. In Figure 5d we did the comparison between and level 2-B. The conclusion is the same, but in this case we add the "noise" coming from the compression algorithm used to build the level 1 data. We emphasize that all these conclusions are only relative to our site, which presents high values of the PW in the atmosphere, up to 60 mm, and sometimes even more. The level of 500 hPa was fixed decades ago by meteorologists for weather forecasting in temperate areas, and it is simply too low for tropical areas, contrary to the opinion of [57]. In tropical and equatorial areas, the water vapor layer reaches higher altitudes than in mid-latitudes areas, and up to 8 km in appreciable quantities [58].  Now, we compare the same sets of level 2-A data taken from the ground to 5500 m with regard to level 2-B data, i.e., data taken from the ground to the maximum altitude of the balloon (over 25,500 m).  Table 4. The bias is 2.75 mm, the STD is 2.32 mm and the RMS is 3.60 mm. In order to better understand the relationship between the balloon altitude and the PW, we map in Figure 5b the level 2 data from a characteristic balloon launch at the Météo-France site, close to our observatory. The balloon reaches its maximum altitude, quasi-linearly with time, in about one hour. It is also clear from Figure 5b that the PW contents, integrated from the surface to the current altitude of the balloon also reaches a horizontal asymptote, and that this asymptote is pretty close, but does not correspond, to the 500 hPa level. About 95% of the maximum value of the PW is acquired after 15 min of ascent, and the last 5% in the remaining part of the ascent. Figure 5c shows the relationship between level 2-A RS PW values and level 2-B RS PW values. Their least-squares linear fit is: level_2_B = 1.09 × level_2-A − 1.35, with R 2 = 94.46%. Figure 5d shows the relationship between level 1 RS PW values and level 2-A RS PW values, and their least-squares linear fit is: level_2-B = 1.14 × level_1 − 3.29, R 2 = 98.46%. Figure 5c and its linear fit indicates that, in the mean, the level 2-A values, in our case, underestimate the PW by about 9% with respect to the level 2-B values. One can argue that this difference is caused by the balloon lateral drift that can reach tens of kilometers when the balloon bursts. But the direction of the drift is highly variable as a function of the wind conditions during the balloon ascent, and this drift is essentially the cause of the scatter of the data points (about 2000 balloon launches) in Figure 5c. Therefore, our conclusion is still valid, in a certain sense. In Figure 5d we did the comparison between and level 2-B. The conclusion is the same, but in this case we add the "noise" coming from the compression algorithm used to build the level 1 data. We emphasize that all these conclusions are only relative to our site, which presents high values of the PW in the atmosphere, up to 60 mm, and sometimes even more. The level of 500 hPa was fixed decades ago by meteorologists for weather forecasting in temperate areas, and it is simply too low for tropical areas, contrary to the opinion of [57]. In tropical and equatorial areas, the water vapor layer reaches higher altitudes than in mid-latitudes areas, and up to 8 km in appreciable quantities [58].

Comparison of the Estimates of Our Tm Model, Bevies Tm Model and vmf1_g
The linear relationship = 70.2 + 0.72 × between Tm and Ts has been proposed by Bevis in 1992 [30], and is now used worldwide, thanks to its simplicity. Nevertheless, in 1997 Ross and Rosenfeld [59] found that this relationship weakens in the equatorial region and becomes even weaker in summer than in winter. It is clear from Figure 6 that the Tm values that we derived by applying Equation (10)

Comparison of the Estimates of Our T m Model, Bevies T m Model and vmf1_g
The linear relationship T m = 70.2 + 0.72 × T s between T m and T s has been proposed by Bevis in 1992 [30], and is now used worldwide, thanks to its simplicity. Nevertheless, in 1997 Ross and Rosenfeld [59] found that this relationship weakens in the equatorial region and becomes even weaker in summer than in winter. It is clear from Figure 6 that the T m values that we derived by applying Equation (10) differs greatly from the Bevis et al. estimates. Therefore, we decided to fit in our case local linear relationships, keeping in mind that the range of variations of our local ground temperature T s is much smaller (8 K) than the range found in the Bevis et al. model (80 K). In Figure 6a, we made a comparison among Tm estimates (Bevis et al. Tm, vmf1_g Tm and our Tm estimate from Equation (10)) over three years, from 2014 to 2016. The vmf1_g Tm estimate is calculated by using Equation (4) based on a 2°× 2.5° global grid with 6-hourly temporal spacing [48]. It is clear from Figure 6 that our Tm estimate shows an increasing uptrend with regard to the Bevis et al. and vmf1_g models, which exhibits much lower values and are almost equivalents. Figure 7 shows that the relationship between Tm and Ts in our case is still a linear one, albeit it is more "fuzzy" than in the work of Bevis et al., and that it also depends weakly on the season (dry or wet).  Table 5 shows the parameters (intercept and slope) of the linear least-squares fit to our LVP Tm corresponding to Figure 7. We applied the usual 3-sigma rejection rule for the rare outliers.  (10)) over three years, from 2014 to 2016. The vmf1_g T m estimate is calculated by using Equation (4) based on a 2 • × 2.5 • global grid with 6-hourly temporal spacing [48]. It is clear from Figure 6 that our T m estimate shows an increasing uptrend with regard to the Bevis et al. and vmf1_g models, which exhibits much lower values and are almost equivalents. Figure 7 shows that the relationship between T m and T s in our case is still a linear one, albeit it is more "fuzzy" than in the work of Bevis et al., and that it also depends weakly on the season (dry or wet). In Figure 6a, we made a comparison among Tm estimates (Bevis et al. Tm, vmf1_g Tm and our Tm estimate from Equation (10)) over three years, from 2014 to 2016. The vmf1_g Tm estimate is calculated by using Equation (4) based on a 2°× 2.5° global grid with 6-hourly temporal spacing [48]. It is clear from Figure 6 that our Tm estimate shows an increasing uptrend with regard to the Bevis et al. and vmf1_g models, which exhibits much lower values and are almost equivalents. Figure 7 shows that the relationship between Tm and Ts in our case is still a linear one, albeit it is more "fuzzy" than in the work of Bevis et al., and that it also depends weakly on the season (dry or wet).  Table 5 shows the parameters (intercept and slope) of the linear least-squares fit to our LVP Tm corresponding to Figure 7. We applied the usual 3-sigma rejection rule for the rare outliers.  Table 5 shows the parameters (intercept and slope) of the linear least-squares fit to our LVP T m corresponding to Figure 7. We applied the usual 3-sigma rejection rule for the rare outliers. A Quantile-Quantile (QQ) plot of the residuals [60] indicates that no remaining structure can be found in the residuals of the fits in Figure 8, except maybe at the ends of the QQ line (see Section 6.2, Figure 12, and conclusions in Section 7).  [60] indicates that no remaining structure can be found in the residuals of the fits in Figure 8, except maybe at the ends of the QQ line (see Section 6.2, Figure 12, and conclusions in Section 7).  Table 5). No remaining structure can be found in the residuals of the fits in Figure 8b versus Figure 8a, except maybe at the ends of the QQ line (see the conclusions for a discussion of the origin of this possible misfit). Figure 9a shows the comparison of our monthly averaged values of LVP PW based on all-seasons Tm model and season-specific (dry and wet) Tm models. The consistency is pretty good, but it is also clear in Figure 9b that the use of separate dry and wet Tm models introduces discontinuities. During the dry season (from May to October), the mean PW differences between our LVP PW estimates computed with the all-seasons linear fit and our LVP PW estimates computed by distinguishing the dry and wet seasons (a linear fit per season) are a little bit closer to zero than the corresponding differences in the wet season (from November to April) (Figure 9b). It reflects the fact that the accuracy of PW estimates based on the wet season-specific Tm model is very close to the accuracy given by using the all-seasons specific Tm model. The LVP PW values derived by using the dry-season-specific Tm model seems to be not so good. This is may be related to the fact that during the dry season (less cloudy), the larger solar radiation heating of the humidity sensor may have caused larger biases [34].  Figure 9a shows the comparison of our monthly averaged values of LVP PW based on all-seasons T m model and season-specific (dry and wet) T m models. The consistency is pretty good, but it is also clear in Figure 9b that the use of separate dry and wet T m models introduces discontinuities. During the dry season (from May to October), the mean PW differences between our LVP PW estimates computed with the all-seasons linear fit and our LVP PW estimates computed by distinguishing the dry and wet seasons (a linear fit per season) are a little bit closer to zero than the corresponding differences in the wet season (from November to April) (Figure 9b). It reflects the fact that the accuracy of PW estimates based on the wet season-specific T m model is very close to the accuracy given by using the all-seasons specific T m model. The LVP PW values derived by using the dry-season-specific T m model seems to be not so good. This is may be related to the fact that during the dry season (less cloudy), the larger solar radiation heating of the humidity sensor may have caused larger biases [34].

Comparison of Mean LVP PW Based on the Local all Seasons T m Model and Local Season-Specific T m Models
(a) (b) Figure 9. Time series from 2014 to 2016 of (a) our local mean LVP PW estimates based on all-seasons Tm specific model (blue dots, see Figure 7a) and our mean LVP PW based on season-specific Tm models (red dots, dry or wet depending on the seasons, see Figure 7b), and (b) the differences between the mean PW estimates in Figure 9a. In order to improve legibility, we smoothed the raw PW time series through an averaging window of one month (Equation (14)) to obtain the mean time series shown in Figure 9a.

Comparison of LVP PW Based on Season-Specific Tm Model with Level 2-B RS PW Values
In order to test the accuracy of the PW estimates based on our new Tm models, we made a comparison between our LVP PW estimates based on our dry and wet season-specific Tm models and level 2-B RS PW values (from the surface to the maximum altitude of the balloon) from 2014 to 2016 (Figure 10a). We can see the consistency is good. Figure 10b shows that the PW differences are close to a Gaussian distribution with no appreciable bias (0.14 mm), an RMS of 3.29 mm and an STD of 3.29 mm, slightly biased toward positive differences (Table 6).   Figure 7a) and our mean LVP PW based on season-specific T m models (red dots, dry or wet depending on the seasons, see Figure 7b), and (b) the differences between the mean PW estimates in Figure 9a. In order to improve legibility, we smoothed the raw PW time series through an averaging window of one month (Equation (14)) to obtain the mean time series shown in Figure 9a.

Comparison of LVP PW Based on Season-Specific T m Model with Level 2-B RS PW Values
In order to test the accuracy of the PW estimates based on our new T m models, we made a comparison between our LVP PW estimates based on our dry and wet season-specific T m models and level 2-B RS PW values (from the surface to the maximum altitude of the balloon) from 2014 to 2016 (Figure 10a). We can see the consistency is good. Figure 10b shows that the PW differences are close to a Gaussian distribution with no appreciable bias (0.14 mm), an RMS of 3.29 mm and an STD of 3.29 mm, slightly biased toward positive differences (Table 6).  Figure 7a) and our mean LVP PW based on season-specific Tm models (red dots, dry or wet depending on the seasons, see Figure 7b), and (b) the differences between the mean PW estimates in Figure 9a. In order to improve legibility, we smoothed the raw PW time series through an averaging window of one month (Equation (14)) to obtain the mean time series shown in Figure 9a.

Comparison of LVP PW Based on Season-Specific Tm Model with Level 2-B RS PW Values
In order to test the accuracy of the PW estimates based on our new Tm models, we made a comparison between our LVP PW estimates based on our dry and wet season-specific Tm models and level 2-B RS PW values (from the surface to the maximum altitude of the balloon) from 2014 to 2016 (Figure 10a). We can see the consistency is good. Figure 10b shows that the PW differences are close to a Gaussian distribution with no appreciable bias (0.14 mm), an RMS of 3.29 mm and an STD of 3.29 mm, slightly biased toward positive differences (Table 6).    We can see that their consistency is pretty good. Figure 11b shows that the PW differences are clearly close to a Gaussian distribution with no appreciable bias (0.04 mm), an RMS of 4.29 mm and an STD of 4.29 mm (Table 7). In term of bias, our all-seasons T m model is better than our season-specific T m model, but it is worse in terms of variance and extreme values.   We can see that their consistency is pretty good. Figure 11b shows that the PW differences are clearly close to a Gaussian distribution with no appreciable bias (0.04 mm), an RMS of 4.29 mm and an STD of 4.29 mm (Table  7). In term of bias, our all-seasons Tm model is better than our season-specific Tm model, but it is worse in terms of variance and extreme values.   We also made a comparison of our monthly averaged LVP PW values based on our all-seasons Tm model and season-specific (dry and wet) Tm models with the monthly averaged values of level 2-B RS PW from 2014 to 2016 (Figure 12a). The consistency is good (Figure 12a), but the differences (Figure 12b), albeit centered around zero, clearly indicate that a seasonal signal is still at play (see Figure 8b). It further reflects that the accuracy of our PW estimates based on our wet season-specific Tm model is better than the PW estimates based on our dry season-specific Tm model (Figure 9b), and in this case the all-seasons Tm model is a good choice to estimate the LVP PW values during the dry season.  We also made a comparison of our monthly averaged LVP PW values based on our all-seasons T m model and season-specific (dry and wet) T m models with the monthly averaged values of level 2-B RS PW from 2014 to 2016 (Figure 12a). The consistency is good (Figure 12a), but the differences (Figure 12b), albeit centered around zero, clearly indicate that a seasonal signal is still at play (see Figure 8b). It further reflects that the accuracy of our PW estimates based on our wet season-specific T m model is better than the PW estimates based on our dry season-specific T m model (Figure 9b  From the above comparisons of our LVP PW values with level 2-B RS PW values, we can conclude that our Tm models according to the all-seasons, dry season and wet season are all reliable for accurate PW estimation in the tropical region of Tahiti. It further reflects that in our case the 500 hPa (5500 m) threshold is one of the main sources of bias between GPS PW and RS PW estimates.

Correlation of the Time Series of PW with the Fluctuations of Pressure and Temperature
We will now consider the correlations of our LVP PW estimates with weather data. We have two sources of weather data (temperature and pressure): The local TAH1 weather station (collocated with the THTI station), but, unfortunately, this weather station was stopped in mid-2015 and replaced by a new one by the end of 2017. The second one is the gridded VMF1 (vmf1_g) files (ECMWF grid weather data). Figure 13a,c shows the comparison between the surface pressure and temperature downloaded from the gridded VMF1 files and the ones recorded at the TAH1 weather station, from 2011 to 2014. We note that the monthly averaged values of TAH1 recorded pressure and temperature is always higher than the vmf1_g values by 0.46 hPa and 0.73 K ( Table 8). The vmf1_g documentation states that their records are relative to the site, and provides the station orthometric height. The observed bias of 0.46 hPa and 0.73 K may be due therefore to an instrumental bias on the sensor of the TAH1 weather station, but it is more likely related to how exactly the vmf1_g estimates are computed, and to the complex interactions between the lower atmosphere and the orography of the Tahiti Island. It cannot be due to the 100 m elevation of the station, as in this case the bias should have had the opposite value. According to these comparisons (Figure 13a,c and Table 8), we think that the surface pressure and temperature from gridded VMF1 files (ECMWF grid weather data) are reliable. Figure  13b,d shows the fluctuations of surface pressure and temperature downloaded from the gridded VMF1 files from 2014 to 2016. We used these surface pressure and temperature to calculate our ZHD values with regard to the Saastamoinen model (Equation (6)) and to model our new linear relationships (Equation (5)) between Tm and Ts for estimating the GPS PW values in Tahiti. From the above comparisons of our LVP PW values with level 2-B RS PW values, we can conclude that our T m models according to the all-seasons, dry season and wet season are all reliable for accurate PW estimation in the tropical region of Tahiti. It further reflects that in our case the 500 hPa (5500 m) threshold is one of the main sources of bias between GPS PW and RS PW estimates.

Correlation of the Time Series of PW with the Fluctuations of Pressure and Temperature
We will now consider the correlations of our LVP PW estimates with weather data. We have two sources of weather data (temperature and pressure): The local TAH1 weather station (collocated with the THTI station), but, unfortunately, this weather station was stopped in mid-2015 and replaced by a new one by the end of 2017. The second one is the gridded VMF1 (vmf1_g) files (ECMWF grid weather data). Figure 13a,c shows the comparison between the surface pressure and temperature downloaded from the gridded VMF1 files and the ones recorded at the TAH1 weather station, from 2011 to 2014. We note that the monthly averaged values of TAH1 recorded pressure and temperature is always higher than the vmf1_g values by 0.46 hPa and 0.73 K ( Table 8). The vmf1_g documentation states that their records are relative to the site, and provides the station orthometric height. The observed bias of 0.46 hPa and 0.73 K may be due therefore to an instrumental bias on the sensor of the TAH1 weather station, but it is more likely related to how exactly the vmf1_g estimates are computed, and to the complex interactions between the lower atmosphere and the orography of the Tahiti Island. It cannot be due to the 100 m elevation of the station, as in this case the bias should have had the opposite value. According to these comparisons (Figure 13a,c and Table 8), we think that the surface pressure and temperature from gridded VMF1 files (ECMWF grid weather data) are reliable. Figure 13b,d shows the fluctuations of surface pressure and temperature downloaded from the gridded VMF1 files from 2014 to 2016. We used these surface pressure and temperature to calculate our ZHD values with regard to the Saastamoinen model (Equation (6)) and to model our new linear relationships (Equation (5)) between T m and T s for estimating the GPS PW values in Tahiti. weather station (cyan dots, with some gaps) and surface pressure as downloaded from the vmf1_g (red dots), and monthly averaged surface pressure recorded by the TAH1 weather station (blue dots, with some gaps) and monthly averaged surface pressure as downloaded from the vmf1_g (magenta dots), and (c) surface temperature Ts from the TAH1 weather station (cyan dots, with some gaps) and monthly averaged temperature Ts from the vmf1_g (red dots), and monthly averaged surface temperature Ts from the TAH1 weather station (green dots, with some gaps) and monthly averaged temperature Ts from the vmf1_g (magenta dots). Right column: Time series from 2014 to 2016 of: (b) surface pressure from the vmf1_g (red dots), and monthly averaged surface pressure from the vmf1_g (magenta dots), and (d) surface temperature Ts from the vmf1_g (red dots), and monthly averaged surface temperature Ts from the vmf1_g (magenta dots). The variations of our mean LVP PW estimates (Figure 12a) are very clearly anti-correlated with the pressure variations ( Figure 13a) and weakly correlated with the temperature variations ( Figure  13c). This indicates that the culprit is probably the gridded VMF1 files (ECMWF grid weather data) that over-smooths the pressure and temperature that are entered in the VMF1 functions with respect to the true local values, and maybe, for high PW values, over-smooths the conversion factor Π of the Equation (3). weather station (cyan dots, with some gaps) and surface pressure as downloaded from the vmf1_g (red dots), and monthly averaged surface pressure recorded by the TAH1 weather station (blue dots, with some gaps) and monthly averaged surface pressure as downloaded from the vmf1_g (magenta dots), and (c) surface temperature T s from the TAH1 weather station (cyan dots, with some gaps) and monthly averaged temperature T s from the vmf1_g (red dots), and monthly averaged surface temperature T s from the TAH1 weather station (green dots, with some gaps) and monthly averaged temperature T s from the vmf1_g (magenta dots). Right column: Time series from 2014 to 2016 of: (b) surface pressure from the vmf1_g (red dots), and monthly averaged surface pressure from the vmf1_g (magenta dots), and (d) surface temperature T s from the vmf1_g (red dots), and monthly averaged surface temperature T s from the vmf1_g (magenta dots). The variations of our mean LVP PW estimates (Figure 12a) are very clearly anti-correlated with the pressure variations ( Figure 13a) and weakly correlated with the temperature variations (Figure 13c). This indicates that the culprit is probably the gridded VMF1 files (ECMWF grid weather data) that over-smooths the pressure and temperature that are entered in the VMF1 functions with respect to the true local values, and maybe, for high PW values, over-smooths the conversion factor Π of the Equation (3).

Conclusions and Outlook
In this research, we investigated, from a metrological point of view, the reliable retrieval of Integrated Precipitable Water Vapor (IPWV, or PW for short), from GPS data, in the tropical island of Tahiti. For this purpose, we did a cross-check between our GPS-derived estimates with the radio soundings data estimates from a nearby site, but also checked the internal consistency of the radio soundings estimates and of the GPS data estimates, in contrast to many studies that assume that the radio sounding data taken from the IGRA archive are "the" reference. In a first step, we checked the internal consistency of the modeling of the ZTD and ZWD estimates from our GPS data processing with respect to the CODE products and the Saastamoinen model. In a second step, we asserted the internal reliability of RS data. The first dataset we used, in this second step, is the archived PW data from the IGRA website (level 1 RS PW from the surface to the 500 hPa level). The second dataset, from the Météo-France archive, is the PW data from raw balloon measurements close to our site, sub-divided into two subsets: a/level 2-A RS PW, from the surface to the 500 hPa level and b/level 2-B RS PW from the surface to the maximum altitude of the balloon. Our analysis clearly indicates that the assumption of taking the IGRA archive as a reference can lead to an underestimation of the water vapor presents in the atmosphere in tropical areas (up to 9% in our case), and that only the RS level 2-B can be compared with the GPS-derived estimates in a consistent way. In a third step, we developed new models of the mean temperature T m of the atmosphere with regard to its water vapor contents, based on a combination of level 2-B RS PW values and GPS ZWD values over Tahiti. The new T m models based on all-seasons and different seasons (dry and wet) were derived by using meteorological data from the gridded VMF1 files. The results show that our all seasons T m model and season-specific (dry and wet) T m models are more reliable in our site (Tahiti) than the Bevis et al. model. In a fourth and last step, we compared our GPS PW estimates based on our new T m models with level 2-B RS PW values. The results show that the PW derived with a "good" T m model from GPS data is highly accurate, and as a final consequence confirms that the threshold of 500 hPa (5500 m) is the main source of bias between the GPS PW estimates and RS PW estimates taken from the IGRA archive in tropical areas.
Historically, the algorithms that are now widely employed to derive PW estimates from GNSS propagation delays can be traced down to the works of Owens in 1967 [61] on the optical refractive index of air, then Thayer in 1974 [62] on improved formulas, followed by Davis et al. in 1985 [27] who introduced the concept of the mean temperature T m of the atmosphere with respect to its water vapor contents. The last essential brick to the edifice was laid down by Bevis et al. in 1992 [30] in the form of a linear relationship from the surface temperature T s at the GNSS receiver site to the mean temperature T m of the overlying column of the atmosphere and the introduction of the proportionality constant Π between the PW estimate and the ZWD estimate. On the GNSS side, Marini in 1972 [63] pioneered the introduction of the mapping functions, which are now culminating with mapping functions tailored to each site and to each epoch with the use of the VMF1 family of mapping functions and the online vmf1_g; the introduction of PPP orbits and precise clock information hammered down these sources of incertitude on the propagation delays [48,50]. A survey of mapping functions can be found in [64]. Nevertheless, the analysis done in this paper reveals some drawbacks, which can be either corrected and/or investigated.
The first one, as mentioned above, is linked to the use, by most of the authors, of the IGRA database that assume an air column limited to the 500 hPa level, i.e., roughly 5500 m. As the scale height of the water vapor is about 2000 m worldwide [65] with regard to a scale height of 8000 m for the whole atmosphere, that sounds reasonable, but this measurement threshold is probably too low in tropical and equatorial areas, were most of the water vapor exchanges between the ocean and the atmosphere take place, as evidenced by our analysis of Section 4. Secondly, the definition of T m , as given by Davis et al. [27] is also derived under the hypothesis that water vapor contents are not so high [27]. In the paper of Bevis et al. [30], it is also apparent in their plot of T m versus T s that the scattering of these couples of values is larger for higher values of T s [30]. This is evidenced by the study done in this paper in Section 5, where we found different linear relationships between T m and T s , with a good correspondence with the radio soundings taken up to the maximum altitude of the balloons, always roughly around 25,500 m. We have no indication in the paper of Bevis et al. of the altitude threshold of the measurements for the determination of the PW estimates, but it was likely also 5500 m. We have to stress that our linear relationships were determined only over a three-year period, and besides, that they are relative to only an 8 K range between the extrema of temperatures (294 K to 302 K) found on the small island of Tahiti, which is surrounded by the thermal buffer of the deep South Pacific ocean, redarding the 80 K range of T s explored by Bevis et al. This, nevertheless, points to the fact that all the chain of state equations of the atmosphere and its optical properties must be revisited for high water vapor contents, not only by adding more GNSS and radio soundings observations in the tropical areas, but also by going back to the laboratory to check again the behavior of the atmosphere for high water vapor contents and maybe by redefining the concept of T m . This is of course far beyond the scope of this study. It is our opinion that trying to derive better empirical-only relations between T m and T s , for example by adding a seasonal variation (probably related to the lack of fit at the ends of the QQ plot line in Figure 8), as done by Yao et al. [29] and Bai [35], is going nowhere without a firm understanding of the physics of the atmosphere for high water vapor contents.