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

Abstract

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.

1. 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 quality and poor spatial coverage [11,12,13]. Besides, these traditional techniques measure the PW distribution only at coarse scales. These limitations form the main source of error in short-term precipitation forecasts [14]. In contrast, PW estimates derived from Global Positioning System (GPS) measurements have been routinely done in the last twenty years with high internal repeatability, all-weather capability, and low cost [15,16]. The spatial resolution of PW by using GPS data depends on the number GPS receivers deployed on the Earth surfaces [17,18]. Nowadays, with the always increasing density of GPS sites, the application of ground-based GPS PW detection is strongly promoted in the meteorological field [19]. Many studies are focused on the PW estimates obtained from GPS technique. Ohtani et al. (2000) have investigated the accuracy and features of PW estimates obtained from Japanese GPS networks and compared them with radiosonde observations [20]. Kwon et al. (2008) have made a comparison between GPS PW estimates derived from a permanent ground-based GPS network on the Korean Peninsula and radiosonde measurements [21]. Choy et al. (2015) have compared PW derived from grounded-based GPS receiver with traditional radiosonde measurement in Australia [22].

GPS radio signals are delayed when they pass through the ionosphere (the upper part of the atmosphere) and the neutral atmosphere. Part of this propagation delay is caused by the atmospheric water vapor, which is mostly concentrated in the troposphere, and this delay can be estimated as a function of the temperature, pressure, and water vapor contents of the atmospheric layers [23,24]. It can be essentially reduced to the knowledge of the zenithal total delay (ZTD), itself divided into two components: zenithal hydrostatic delay (ZHD) and zenithal wet delay (non–hydrostatic) (ZWD) [25,26]. The GPS ZWD can be converted to a PW estimate by a conversion coefficient$\mathsf{\Pi }$, which is related to the weighted mean temperature of water vapor (T_m) [27,28,29]. A standard processing is to compute T_m, then $\mathsf{\Pi }$ and finally PW from ZTD minus ZHD. It is clear that, from a mathematical point of view, this standard processing defines an iterative procedure, as we need an a priori guess of the water vapor contents of the atmosphere to determine the PW estimate. In practice, T_m is expressed as a linear function of surface temperature (T_s).

The ultimate goal is to get the slope and intercept of this linear function almost site-independent. This linear relationship has been explored by many scientists, both globally and regionally, from an a priori knowledge of the water vapor contents of the atmosphere from radio soundings. In 1992, Bevis et al. [30] came up with $T_m=70.2+0.72\times T_s$ from the analysis of 8718 profiles of radiosonde launches at 13 stations in the US over two years. This relation has been widely used in North Hemisphere mid-latitude regions between 27°N and 65°N, and was revised in 1994 by Bevis et al. [31], from a nearly global distribution of about 250,000 radiosonde profiles. In 1998, Mendes et al. [32] proposed another linear T_m model by using about 32,500 radiosonde profiles over one year at 50 sites between 62°S and 83°N. A dedicated model is available for high latitudes [33]. Liou et al. [34], Bai [35], Wang et al. [36], Raju et al. [37], Emardson et al. [38] also established regional and global linear relationships between T_m and T_s [39,40,41,42].

Of course, if we have an on-site good estimate of both the PW from an external source (essentially radio soundings, but this can be also absorption lines of water vapor in the atmosphere or radiometric observations), and ZWD from a collocated GPS receiver, then $\mathsf{\Pi }$and T_m can be directly derived. If, furthermore, we have a collocated weather station giving us T_s, then we can obtain by linear regression over a time series of T_m and T_s site-dependent values of the slope and intercept. This is basically the approach we applied in this paper, as we had the luxury to have a radio sounding station close to the Geodesy Observatory of Tahiti (OGT), with an International GNSS Service (IGS)-grade GPS receiver with weather data.

In this paper, the dataset and methodology are described in Section 2. The comparisons of GPS-derived results with the Center for Orbit Determination in Europe (CODE) products [43], as well as results obtained by using the Saastamoinen model are given in Section 3. In Section 4, we compare PW values downloaded from IGRA website and PW values computed from balloon raw data. The comparisons between our all-seasons T_m model and our season-specific T_m models with respect to 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.

2. Study Area, Methodology and Datasets

2.1. 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²) 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.

2.2. 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):

$$ZTD=ZHD+ZWD$$

(1)

And, the relationship between ZWD and PW [27] can be expressed as:

$$PW=\mathsf{\Pi }·ZWD$$

(2)

where

$$\mathsf{\Pi }=\frac{{10}^6}{\rho ·R_v(\frac{k_3}{T_m}+k_2^{\prime })}$$

(3)

and

$$T_m=\frac{{\displaystyle \int \frac{P_v}{T}dz}}{{\displaystyle \int \frac{P_v}{T^2}dz}}$$

(4)

where $\rho =1000\mathrm{kg}/{\mathrm{m}}^3$ is the density of liquid water, $R_v=461.495\mathrm{J}/\mathrm{Kg}·\mathrm{K}$ is the specific gas constant of water vapor, $k_2^{\prime }=22.1\mathrm{K}/\mathrm{mb}$ and $k_3=3.739·{10}^5{\mathrm{K}}^2/\mathrm{mb}$ are physical constants. The coefficient $\mathsf{\Pi }$ 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):

$$T_m=a+b·T_s$$

(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:

$$ZHD=\frac{2.2768P_s}{f(\lambda ,H)}$$

(6)

where $f\left(\lambda ,\mathsf{{\rm H}}\right)=1-0.00266·\cos \left(2\lambda \right)-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.

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^{°}\times {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.unibe.aiub.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.

2.3. 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:

$$PW=\frac{1}{\rho }{\displaystyle \int {\rho }_w}dh$$

(7)

with

$${\rho }_w=\frac{RH}{100·R_v·T}·e_s$$

(8)

and

$$e_s=f·6.1121·e^{\frac{(18.729-\frac{t}{227.3})·t}{t+257.87}}$$

(9)

where $h$ is the geopotential altitude, $\rho$is the density of liquid water, ${\rho }_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.

2.4. 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” $\mathsf{\Pi }$according to Equation (10), which is the reversed version of Equation (2):

$$\frac{PW\_RS}{ZWD\_GPS}=\mathsf{\Pi }$$

(10)

Then the time series of $\mathsf{\Pi }$ 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).

2.5. 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:

$$Bias=\frac{1}{n}{\displaystyle \sum _{i=1}^n(S(i)-R(i))}$$

(11)

$$STD=\sqrt{\frac{1}{n}{\displaystyle \sum _{i=1}^n{((S(i)-R(i))-Bias)}^2}}$$

(12)

$$RMS=\sqrt{\frac{1}{n}{\displaystyle \sum _{i=1}^n{(S(i)-R(i))}^2}}$$

(13)

where $i=1,2,3,\dots ,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]:

$$\overline{S}[i]=\frac{{\displaystyle \sum _{j=1}^{n}S[j]*v_{ij}}}{{\displaystyle \sum _{j=1}^n{v}_{ij}}}$$

(14)

where $S\left[j\right]$ is the value of S corresponding to the j-th epoch. We will always use an averaging over one month, $t=\mathrm{MJD}$ (Modified Julian Day), $\Delta t=t_j-t_i$, $i,j=1,2,3,\dots ,n$, if $\left|\Delta t\right|>T$, we set $v_{ij}=0$, otherwise,$v_{ij}=1$.

3. Comparisons of ZTD and ZWD Values from Our GPS Data Processing, CODE Products and Saastamoinen Model

3.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–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. Statistical summary for the comparison between our LVP ZTD and ZWD and CODE products (LVP minus CODE) in terms of max, min, bias, STD, and RMS, from 2011 to 2013, relative to Figure 2b,d.

Data Differences	Max (mm)	Min (mm)	Bias (mm)	STD (mm)	RMS (mm)	Data Points
ZTD	19.91	−19.80	−0.24	6.10	6.11	12,702
ZWD	20.17	−19.64	−0.16	6.11	6.12	12,702

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.

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.

3.2. 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. 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.

Data Differences	Max (mm)	Min (mm)	Bias (mm)	STD (mm)	RMS (mm)	Data Points
LVP-SAAS_1	4.95	−10.00	0.85	1.21	1.48	4178
CODE-SAAS_2	4.72	−3.57	0.69	1.16	1.35	4178

summarize the results with a comparison with the ZWD estimates from the previous section (VMF1 mapping functions) and CODE products.

Table 1 and Table 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.

4. 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.

4.1. 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 4a shows the time series of level 1 RS PW and level 2-A RS PW and the corresponding monthly averaged values from 2014 to 2017. 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² = 94.62%. Figure 4a and

Table 3. Statistical summary for the comparison between level 1 RS PW estimates and level 2-A RS PW estimates (level 1 minus level 2-A) in terms of max, min, bias, STD, and RMS, from 2014 to 2016.

PW Differences	Max (mm)	Min (mm)	Bias (mm)	STD (mm)	RMS (mm)	Data Points
Level 1-2-A	6.21	−5.22	−0.10	1.30	1.30	2151

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.

4.2. Comparison of Level 2-A RS PW Values (from the Surface to the 500 hPa Level) with Level 2-B RS PW Values (from the Surface to the Maximum Altitude of the Balloon)

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). Figure 5a map level-2A versus level 2-B RS PW values over three years, with basically two launches of balloons per day. The level 2-B RS PW values are always larger than the corresponding level 2-A values. The statistic of their differences can be found in

Table 4. Summary for the comparison between the level 2-A RS PW estimates from the surface to the 500 hPa level and level 2-B RS PW estimates from the surface to the maximum altitude of the balloon (the latter minus the former) in terms of max, min, bias, STD and RMS, from 2014 to 2016, relative to Figure 5b.

PW Differences	Max (mm)	Min (mm)	Bias (mm)	STD (mm)	RMS (mm)	Data Points
Level 2	14.19	−19.24	2.75	2.32	3.60	2179

. 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² = 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² = 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].

5. All-Seasons T_m Model, Dry and Wet Season-Specific T_m Models for Our Site

5.1. 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\times 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 T_m estimates (Bevis et al. T_m, vmf1_g T_m and our T_m estimate from Equation (10)) over three years, from 2014 to 2016. The vmf1_g T_m estimate is calculated by using Equation (4) based on a $2°\times 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).

Table 5. Intercept and slope of our new on-site T_m models for all-seasons and the two different tropical seasons (dry and wet). Number of data points fitted: 2060, number of data points deleted by applying the 3-sigma rule: 33.

Tm Model	Intercept	Sigma	Slope	Sigma
All seasons	−476.44	145.55	2.68	0.49
Dry season	−535.69	252.91	2.89	0.85
Wet season	−979.71	226.86	4.36	0.76

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).

5.2. Comparison of Mean LVP PW Based on the Local all Seasons T_m Model and Local Season-Specific T_m Models

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].

6. Comparison of LVP PW Values with the Level 2-B RS PW Values

6.1. 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. Statistical summary of the comparison between our LVP PW estimates based on our season-specific T_m models and level 2-B RS PW (LVP PW minus RS PW) in terms of max, min, bias, STD and RMS, from 2014 to 2016, relative to Figure 10b.

PW Differences	Max (mm)	Min (mm)	Bias (mm)	STD (mm)	RMS (mm)	Data Points
LVP-level 2-B RS	19.29	−14.91	0.14	3.29	3.29	2060

).

6.2. Comparison of Our LVP PW Values Based on Our all Seasons T_m Model with Level 2-B RS PW Values

Figure 11a shows the comparison of LVP PW estimates based on the all-seasons specific T_m model with the level 2-B RS PW estimates for the 2014–2016 period. 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. Statistical summary for the comparison between our LVP PW estimates based on our all-seasons T_m model and level 2 RS PW values (LVP minus RS) in terms of max, min, bias, STD, and RMS, from 2014 to 2016, relative to Figure 11b.

PW Differences	Max (mm)	Min (mm)	Bias (mm)	STD (mm)	RMS (mm)	Data Points
LVP-RS	53.87	−27.50	0.04	4.29	4.29	2060

). 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 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), and in this case the all-seasons T_m model is a good choice to estimate the LVP PW values during the dry season.

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.

6.3. 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. Statistical summary for the comparison between monthly averaged surface pressure and temperatures T_s from the gridded VMF1 files and from TAH1 weather station (vmf1_g minus TAH1) in terms of max, min, bias, STD, and RMS, from 2011 to 2014, relative to Figure 13a,c.

Mean Differences	Max	Min	Bias	STD	RMS	Data Points
Pressure (hPa)	−0.21	−1.20	−0.46	0.14	0.48	4882
Ts (K)	0.15	−1.51	−0.73	0.36	0.81	4882

). 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.

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 $\mathsf{\Pi }$ of the Equation (3).

7. 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 $\mathsf{\Pi }$ 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.