IWV observations in the Caribbean Arc from a network of ground-based GNSS receivers during EURECA

Ground-based Global Navigation Satellite System (GNSS) measurements from nearly fifty stations distributed over the Caribbean Arc have been analysed for the period 1 January-29 February 2020 in the framework of the EURECA (Elucidate the Couplings Between Clouds, Convection and Circulation) field campaign. The aim of this effort is to deliver high-quality Integrated Water Vapour (IWV) estimates to investigate the moisture environment of mesoscale cloud patterns in the Tradewinds and their feedback on the large-scale circulation and energy budget. 5 This paper describes the GNSS data processing procedures and assesses the quality of the GNSS IWV retrievals from four operational streams and one reprocessed research stream which is the main data set used for offline scientific applications. The uncertainties associated with each of the data sets, including the zenith tropospheric delay (ZTD) to IWV conversion methods and auxiliary data, are quantified and discussed. The IWV estimates from the reprocessed data set are compared to the Vaisala RS41 radiosonde measurements operated from the Barbados Cloud Observatory (BCO) and to the measurements 10 from the operational radiosonde station at Grantley Adams international airport (GAIA), Bridgetown, Barbados. A significant dry bias is found in the GAIA humidity observations with respect to the BCO sondes (-2.9 kgm−2) and the GNSS results (-1.2 kgm−2). A systematic bias between the BCO sondes and GNSS is also observed (1.7 kgm−2) where the Vaisala RS41 measurements are moister than the GNSS retrievals. The IWV estimates from a collocated microwave radiometer agree with the BCO soundings after an instrumental update on 27 January, while they exhibit a dry bias compared to the soundings and 15 to GNSS before that date. IWV estimates from the ECMWF fifth generation reanalysis (ERA5) are overall close to the GAIA observations, probably due to the assimilation of these observations in the reanalysis. However, during several events where strong peaks in IWV occurred, ERA5 is shown to significantly underestimate the GNSS derived IWV peaks. Two successive peaks are observed on 22 January and 23/24 January which were associated with heavy rain and deep moist layers extending

ultra-rapid satellite clocks and orbits, and processes the measurements in 6-hourly windows shifted by one hour every hour.
For the permanent GNSS stations in the region, the ZTD estimates for the most recent hour are available within 45 min after the end of measurements and are distributed to Numerical Weather Prediction (NWP) centres via the EUMETNET-GNSS 85 water vapour program (E-GVAP). Although the BCO stations were included in this stream they are not registered as permanent stations at E-GVAP and are thus not assimilated by the NWP centres. Due to the sliding window processing approach, each hourly time slot is thus processed six times. The accuracy of the corresponding ZTD estimates depends on the time slot as will be illustrated in the next Section. The main reason is that the quality of the NRT clocks and orbits is lower in the more recent time slot because they are predicted (observations are not available in due time to constrain that part of the orbit calculation). 90 The sampling of the ZTD estimates of the NRT solution is 15 min. On the other hand, the 'daily' processing stream is updated every day before 05 UTC. It is based on a longer time windows of 24 hours and benefits from slightly improved clock and orbit products. Its accuracy is thus slightly higher. The sampling of the ZTD estimates of the daily solution is 1 h. This processing stream is not sent to E-GVAP but is rather devoted to an operational quality control of the observations from the RGP stations.
It is also used for a posteriori verification of the NRT solution. 95 The data from the BCO stations were also included in two other operational streams operated by ENSTA_B/IPGP in the framework of ACTRIS-France (Aerosol, Cloud and Trace Gases Research Infrastructure -France, https://www.actris.fr/, last access: 29 January 2021). These two streams are referred to as 'GIPSY rapid' and 'GIPSY final'. They differ only by the clock and orbit products which are released with daily and fortnightly latency, respectively. The GIPSY rapid processing is equivalent to the RGP daily processing, while the GIPSY final processing on the other hand provides a higher level of accuracy. Both a 100 streams provide ZTD estimates with a sampling of 5 minutes.
Apart from the usage of different clock and orbit products, the processing streams run by both analysis centres also differ by their software packages and related processing options as summarised in Table 1. Most importantly, IGN uses the Bernese GNSS version 5.2 software in double-difference processing mode (Dach et al., 2015) while ENSTA_B/IPGP use the GIPSY OASIS II version 6.2 software in precise point positioning (PPP) mode (Zumberge et al., 1997). Various studies have compared 105 the results and discussed the benefits and drawbacks of both approaches. When consistent clock and orbit products are used, there is generally a good level of agreement (i.e. a few mm RMS differences in ZTD estimates).
Finally, the regional network including 49 GNSS stations was reprocessed by ENSTA_B/IPGP using again GIPSY OASIS II version 6.2 software, in a very similar scheme as the final stream but with a few improved processing options. This processing is expected to be the most accurate. It is referred to as 'GIPSY repro1' in the following.

GNSS data quality checking
The quality of the data processing results can be checked with two types of information: 1) global information quantifying the accuracy of the processing in the 6-hourly or 24-hourly time window (e.g. reduced sum of squares of residuals, percentage of solved ambiguities) and 2) the formal errors of the estimated parameters (station coordinates and ZTDs).
For the operational streams, only the second type of information was available and a basic screening procedure was thus 115 applied in the form of a range check on ZTD estimates with limits [1 m, 3 m] and on for ZTD formal errors with a limit of 0.02 m.
For the GIPSY repro1 solution, a more elaborated quality check was made using both types of information for all stations.
Mean statistics over the 2-month period are reported in Table S1 of the supplemental material. The mean RMS of residuals for all 49 stations ranged between 0.009 m and 0.015 m. Five stations with values above 0.013 m may be unreliable (N240, 120 DEHA, AIRS, BOUL, HABL). These stations have also small mean percentage of fixed ambiguities (≤ 50 %) and/or significant temporal variations in both parameters. Station position repeatability was in general good both in the horizontal (≤ 0.005m) and vertical (≤ 0.015m) except for two stations (BOUL, OLVN). In a second step, a tighter range check completed with an outlier check were applied following the recommendations of Bock (2020a). The range check limits were [0.5 m, 3.0 m] on ZTD values and 0.01 m on formal errors. The outlier check limits were computed for each station based on its median µ and 125 standard deviation σ statistics over the 2-month period according to [µ − 5σ, µ + 5σ] for the ZTD values and [0, µ + 3.5σ] for the formal errors. This screening procedure rejected 0.21 % of the ZTD data due exclusively to the formal error outlier check.
The GIPSY ZTD time series are generally quite small with few outliers as confirmed by the low outlier detection rate.

GNSS ZTD to IWV conversion
The ZTD is the sum of the zenith hydrostatic delay (ZHD) and the zenith wet delay (ZWD): During the data processing, ZHD is fixed to an a priori value, either from an empirical model or from a NWP model analysis (Boehm et al., 2007), and ZWD is estimated. Errors in the a priori ZHD can propagate to the estimated ZWD but, in general, this error is small (< 1 mm) and the sum of both is very close to the real ZTD. When IWV is to be extracted, it is desirable to use a more accurate ZHD value computed from the surface air pressure: where k 1 is the dry air refractivity coefficient, R d the dry air specific gas constant, P s surface air pressure, and g m the mean acceleration due to gravity (Davis et al., 1985). IWV is derived from ZWD by applying a delay to mass conversion factor κ(T m ): κ is a semi-empirical function of the weighted mean temperature T m defined as (Bevis et al., 1992): where k 2 and k 3 are refractivity coefficients for the water molecule, R v is the specific gas constant for water vapour. T m is defined as: where ρ v (z) and T (z) and the specific mass of water vapour and the air temperature, respectively, at height z above the surface. The integrals are from the surface to the top of the atmosphere.
Practically, the GNSS IWV estimates are converted from the ZTD data using equations (1) to (4) with the auxiliary data, P s and T m , given at the position of the GNSS station with interpolation to the GNSS times. On average for the 2-month period (January-February 2020), the conversion parameter amounts to κ = 164 kg m −3 and the observed ZW D = 0.20 m results into 150 IW V = 32 kg m −2 . The sensitivity of IWV estimates to errors in the auxiliary data can be assessed by the partial derivatives: where the numerical values were computed from the average EUREC 4 A conditions. Various sets of refractivity coefficients and auxiliary data were available and used with our operational and research GNSS products and are described in the subsections below.

RGP NRT and daily results
This dataset is provided only for the BCON and BCOS stations which benefit from colocated surface air pressure and temperature measurements. The measurements come from Vaisala PTU200 Sensors connected directly to the GNSS receivers and are included in the raw GNSS data files. The ZHD estimates are therefore computed using Eq. 2 and the T m values are computed using the widely used Bevis et al. (1992) formula: where T s is the surface pressure. This formula was derived by Bevis et al. (1992) from a radiosonde data set over the USA with a RMS error of 4.7K. The refractivity coefficients used for the computation of ZHD and κ are from Thayer (1974), the thermodynamics constants are from ICAO (1993), and the g m formula is from Saastamoinen (1972).

165
This dataset is also provided only for the BCON and BCOS stations but it used the a priori values for ZHD available during the data processing. The values were obtained from the Technical University of Vienna (TUV) as part of the VMF1 tropospheric parameters computed from the ECMWF operational analysis (Boehm et al., 2006). The TUV also computes T m values using Eq. 5 from the same analysis. Both products are made available on a global grid with 2°latitude and 2.5°longitude resolution every 6 hours. The refractivity coefficients and g m are the same as for the RGP data set.

GIPSY repro1 results
For the reprocessed data set we used a more rigorous approach following Bock (2020a). Here the ZHD and T m values are computed from ERA5 pressure-level data at the four surrounding grid points and are interpolated bilinearly to the positions of the GNSS antennas. The ERA5 fields were used with the highest temporal and spatial resolutions: 1-hourly and 0.25°x 0.25°, respectively. The refractivity coefficients and thermodynamics constants were updated from Bock (2020a) to account for the global CO 2 content for January 2020 (see the Appendix). The g m formula proposed by Bosser et al. (2007) was used instead of the Saastamoinen (1972) one.
This methodology provides the most accurate IWV estimates regarding the available data at the 49 stations of the extended network and the best knowledge of uncertainties due to the various empirical formulas (T m and g m ) and auxiliary data.

180
The three data sets (RGP, GIPSY operational, and GIPSY repro1) used different conversion methods, which are associated with different uncertainties. Although repro1 is the most accurate, it is not an operational stream and covers only a limited period of time. Users may thus be interested by the near real time and extended time coverage of the two operational data sets available for stations BCON and BCOS. The consistency of the different data sets is assessed in the next Section for ZTD and IWV. Here we describe how the different error sources contribute to overall uncertainty in the IWV data.

185
Let us consider two ZTD solutions from two different processing streams, ZT D 1 and ZT D 2 , converted to IWV using different ZHD and κ data, ZHD 1,2 and κ 1,2 . The two IWV solutions, IW V 1 and IW V 2 write: The difference ∆IW V = IW V 2 − IW V 1 can be separated into systematic and random components and written as a func-190 tion of the contributions (assumed independent) of the systematic and random differences of the three parameters ∆ZT D = ZT D 2 − ZT D 1 , ∆ZHD = ZHD 2 − ZHD 1 , ∆κ = κ 2 − κ 1 as: The magnitude of the actual errors in the auxiliary data used in the different conversion methods can be appreciated from 195 Fig. 2. The Figure shows that the ZHD data computed from ERA5 reanalysis and PTU sensor agree well. The mean pressure (ZHD) difference, ERA5 -PTU, amounts to −0.13 ± 0.20 hPa (−0.30 ± 0.45 mm), which impact on IWV represents 0.045 ± 0.070 kg m −2 . ERA5 is thus a good source of surface pressure for the IWV conversion with negligible bias at this site. The ZHD estimates computed from ECMWF by TUV on the other hand have significant aliasing errors (Fig. 2). This is due to the strong atmospheric tide seen in the surface pressure that cannot be resolved by the 6-hourly analysis data (shown by crosses 200 upon the red dashed line in the lower ZHD panel).
Regarding the κ data, TUV and ERA5 agree well but the PTU-derived curve has significant errors because it uses the empirical formula given by Eq. 6 which is based on surface temperature and does not well correlate with the upper air variations that influence T m . In the case of the RGP data set, the use of Bevis et al. (1992) formula is the dominant error source in the IWV conversion process, while in the case of GIPSY rapid and final results, the use of 6-hourly ZHD data from TUV is the 205 main error source.
We also evaluated the accuracy of the ERA5 T m estimates in comparison to T m computed from the radiosonde data at the BCO. The mean difference and standard deviation of the resulting IWV estimates for 138 soundings were 0.03 kg m −2 (ERA5 -soundings) and 0.05 kg m −2 , respectively.
Compared to the differences in IWV estimates resulting from the use of different auxiliary data, the impact of the refractivity 210 constants (Thayer (1974) vs. Bock (2020a)) is rather a small and represents a bias of 0.045 kg m −2 .

ZTD and IWV intercomparisons at the BCO
In this Section we inter-compare both ZTD and IWV results from the five processing streams. The motivation for comparing the ZTD results is that it reflects the uncertainty due to the GNSS processing only while the IWV comparison includes the conversion errors discussed in the previous Section. The uncertainty in the ZTD data is of interest to the data assimilation 215 community since ZTD are currently assimilated in most NWP models (rather than IWV). On the other hand, the uncertainty in the IWV data may be of interest to users who wish to analyse the IWV directly (e.g. for process studies, intercomparison with other observational techniques, or verification of atmospheric model simulations).

RGP NRT results
As explained above, in the RGP NRT solution, each hourly time slot is processed six times while it is shifted from the most 220 recent position (referred to as NRT.6) to the oldest (NRT.1) in the 6-hourly time window. The highest precision is expected in the central time slots (NRT.2 to NRT.5) because these are farther from the edges and are best constrained by the observations.

ZTD intercomparison
The overall uncertainty in the NRT solution is mainly due to the ultra-rapid satellite products and small window-size (impacting the ambiguity resolution). It can be evaluated by comparison with the RGP daily or any of the GIPSY solutions. The upper part 225 of Table 2 shows the results for the ZTD estimates. The comparison to RGP daily reveals a trend in the mean difference from NRT.1 to NRT.6, with a large negative bias in NRT.1 of -0.04 m (RGP NRT ZTD < RGP daily ZTD). This bias was actually unexpected. It is an artifact due to the propagation of tropospheric gradient biases from one time slot to the next. This feature was corrected in July 2020 and the bias is no longer present in the current operational NRT product. The standard deviation of differences shows a minimum (9.4 mm) for the NRT.4 solution. The values are actually smaller in the central time slots as 230 expected. However, the formal errors (top row of Table 2) predict smaller uncertainty in the first time slot instead. This is a consequence of using the previous normal equations.
The comparison with GIPSY repro1 confirms these conclusions although there are small differences in the results: mainly, the bias is slightly offset by 4 mm. Since this comparison involves two different software and processing approaches, these small differences are not surprising. Note that the number of comparisons are also different by a factor of four due to the 235 different ZTD sampling of the three solutions. Compared to GIPSY repro1, the RGP NRT.4 has the smallest bias and standard deviation. We recommend thus to use this solution for near-real-time applications when timeliness is not too restrictive (e.g. in NWP assimilation schemes, the ZTD data from NRT.6 are used which are a bit less accurate). Figure 3 shows the ZTD and formal error time series from the RGP NRT.4 solution and the other two processing streams. It it seen that the temporal variations in ZTD are very consistent among all three data streams with small differences (< 0.01 m) 240 compared to the ZTD variations ( 0.1 m). On the other hand, the formal errors of the NRT solution are much larger than those of the RGP daily and the GIPSY repro1 solutions. GIPSY repro1 has also more stable formal errors than RGP daily.

IWV intercomparison
The lower part of Table 2 shows the IWV results. In the case of the comparison of RGP NRT to RGP daily, the IWV results are consistent with the ZTD results since both solutions use the same ZHD and K conversion data. The IWV differences are 245 here proportional to the ZTD differences with a factor of 162 kg m −3 which is close to the mean value of K=164 kg m −3 . On the other hand, the comparison with GIPSY repro1 includes the differences in the conversion data, and the ratio of IWV to ZTD results is not a constant factor. This point is further discussed in the next sub-section. Similar to the ZTD results, the IWV solution for NRT.4 is in good agreement with RGP daily and GIPSY repro1. Figure 4 compares the NRT.6 IWV estimates to GIPSY repro1 for the time period from 10 to 20 January 2020 marked by a 250 period of two days with high IWV around 45 kg m −2 . It is striking that the NRT.6 solution shows spurious oscillations with a period about 1 hour. This feature is due to the strong relative constraint (1 mm) in the NRT solution with an hourly update. The constraint is the same in the NRT.4 slot but the oscillation is strongly damped (not shown). This result also militates for using the NRT.4 solution rather than the NRT.6 solution when relevant.

255
Similar ZTD and IWV comparisons have been performed for the other processing streams which are of interest for offline applications.
The left part of Table 3 shows the ZTD comparison results. The comparison of RGP daily to GIPSY rapid and repro1 show very similar results. The agreement between the two processing software is quite good, with a small difference in the mean values of 2.1-2.4 mm for station BCON and 1.1-1.3 mm for station BCOS. The standard deviations of ZTD differences is about 260 5.7-5.8 mm for both stations. The differences reflect mainly the impact of using different satellite products, mapping functions, tropospheric model constraints, and elevation cutoff angles. The comparison of the three GIPSY solutions show comparatively much better agreement with almost no bias (the mean differences are around -0.1 to 0.3 mm) and very small standard deviations (0.8 to 2.8 mm). The comparison of the two operational solutions (GIPSY rapid and final) shows the smallest standard deviation because they use exactly the same processing options; the only difference is in the satellite products but they were produced by 265 the same analysis centre (JPL) and are highly consistent. The comparison of the GIPSY operational solutions to repro1 shows slightly larger differences which are due to small differences in the processing options.
IWV comparisons have been performed for the same five combinations. In addition to the ZTD differences due to the processing strategies these comparisons will include the effect of the IWV conversion methods: -RGP daily and GIPSY rapid differ by ZHD (PTU vs. TUV) and K (Bevis vs. TUV); strong impact is expected from 270 differences in both parameters; -RGP daily and GIPSY repro1 differ by ZHD (PTU vs. ERA5) and K (Bevis vs. ERA5); strong impact is expected from K; -GIPSY rapid and final use the same conversion parameters; -GIPSY rapid and final differ from repro1 by ZHD and K (TUV vs. ERA5); strong impact is expected from ZHD;

275
The right part of Table 3 shows the IWV results and confirms they include additional errors compared to the ZTD results (left part the Table). Most notably: (i) a change in the bias between RGP daily -GIPSY rapid and RGP daily -GIPSY repro1 because the two GIPSY solutions use different parameters. In the first comparison the bias is due to the ZHD data used in the GIPSY rapid solution (TUV) and the K data in the RGP daily solution (due to Bevis et al. (1994) formula). In the second comparison, only the K bias in the RGP rapid solution is involved. (ii) a noticeable bias and standard deviation between the 280 two operational GIPSY solutions and the repro1 due to the ZHD data used in the former (TUV). The two operational GIPSY solutions are highly consistent in IWV (similar to the ZTD solutions) but this result hides the fact that they include common ZHD errors. Table 4 quantifies more precisely the contributions from the difference error sources (ZTD, ZHD, and K) to the total IWV differences, where the mean and random components have been separated according to Eqs. 9 and 10. It can be seen that:

285
-In the RGP daily and GIPSY rapid comparison, all three error sources contribute in the same proportions to the bias, while the random errors are dominated by ZTD errors.
-In the RGP daily and GIPSY repro1 comparison, the small bias in RGP daily is a result of almost exact cancelling of a ZTD bias and a K bias in the RGP data. The random errors in the RGP data are again dominated by ZTD errors.
-In the GIPSY final and repro1 comparison, the bias in the GIPSY final is almost only due to a bias in ZHD, while the 290 random errors are due to ZTD and ZHD differences of the same magnitude.
In conclusion, compared to the GIPSY repro1 dataset which uses the most accurate data processing and conversion parameters, both operational streams show small systematic and random errors on the level of ± 0.5 kg m −2 and ± 1.0 kg m −2 , respectively.

GNSS compared to other IWV data sources 295
Two other instruments operated at the BCO facility are used here: the Vaisala RS41/MW41 radiosonde system (Stephan et al., 2020) and the RPG, HATPRO-G5, microwave radiometer (MWR) (Rose et al., 2005). The IWV retrievals from these systems were compared to IWV estimates from the BCON GNSS station (GIPSY repro1) and ERA5 IWV data for the period from 1 January to 29 February 2020. For the radiosondes, the level-1 pressure (P), temperature (T), and relative humidity (RH) data were used and IWV was computed as: where the integral extends from the GNSS station height to the top of the sounding ( 24 km on average), q(P ) is the specific humidity computed from RH and T using Tetens (1930) saturation pressure formula over water, g(P ) is the acceleration of gravity as a function of altitude. The 1-sec time sampling of the radiosonde data provide high vertical resolution (about 5 m in the lower troposphere). They were checked for consistency and thinned for including only increasing altitude points.

305
Very few data gaps were noticed in the vertical profiles which confirms that the sounding operations worked fine all along the campaign and did not need further correction or screening. The BCO soundings provided a nearly continuous sampling of the upper air conditions between 16 January and 19 February, every 4 hours (Stephan et al., 2020). The MWR worked for a more extended period from 1 January to 15 February. The brightness temperature measurements were used to retrieve IWV using a neural network algorithm provided by the manufacturer. In precipitating conditions, the measurements usually experience 310 strong biases due to wet radome emissions (see Fig. 5) and are screened out according the a rain detection index (including also sea spray detection). The IWV contents are nominally retrieved with a 1-sec sampling which we down-scaled to 5-min, by computing arithmetic means to be consistent with the GNSS sampling. Some outliers remained in the MWR IWV series which were subsequently removed from the comparisons by a simple outlier check of the IWV differences with limits mean ± 3 standard deviations).

315
The ERA5 IWV contents above the GNSS station were computed from the hourly pressure level data at the four surrounding grid points and interpolated bilinearly to the position of the GNSS antenna. The use of pressure level data instead of model level data induces a minor bias of 0.2 ± 0.1 kg m −2 but these data are more convenient to use and consistent with the Tm In addition, we also included the twice-daily (00 and 12 UTC) radiosonde measurements from the operational radiosonde 320 station at Grantley Adams International Airport (GAIA, WMO code 78954) located 11 km away from the BCO. This station used GRAW DFM-09 at the time of the campaign (Kathy-Ann Caesar, personal communication). The sonde data were retrieved from the University of Wyoming sounding archive and contained on average 85 ± 6 vertical levels up to an altitude of 29 ± 2 km. The vertical sampling of these data is coarser than the BCO data but fine enough to compute proper IWV contents. The latter were computed in a similar way as for the BCO soundings except near the surface where the ERA5 pressure level data 325 were used to complete the sounding profiles. This is because the GAIA station is located at an altitude of 57 m, slightly above the BCO which is at 25 m altitude. ). The agreement between the five data sets is good in general but ERA5 is seen to underestimate the IWV during the aforementioned peaks. Quite obvious is also the bias between the two sounding data, with the BCO IWV contents systematically higher than the GAIA measurements. Less easy to distinguish, but nevertheless significant, is the offset in the HATPRO MWR IWV retrievals before and after the 27 January. The HATPRO data are actually more consistent with the GAIA data before that date and instead more consistent with the BCO radiosonde data after. Analysis of the measured brightness temperatures 335 showed that during the former period, the measurements were less accurate due to an instrumental malfunctioning which was fixed on 27 January. Results from the GPS, the BCO and GAIA sondes, and the ERA5 comparisons are reported in Fig. 6. While IWV varies from 20 to 55 kg m −2 over the period, the bias between the two sondes is -2.89 kg m −2 (GAIA -BCO) with a slope of 0.94 and an offset of -0.69 kg m −2 . The GAIA data are drier than the BCO data over the full observation range. The comparison of profiles shows that the humidity measurements from the two sondes differ mainly in the lowermost 2.5 km where the mean difference in q is larger than 1 g kg −1 . One contribution to this difference may be the difference of trajectories of the sounding 350 balloons released from the two sites. Balloons released from the BCO site usually travel west and southwards over the Barbados island until they reach an altitude of 6-8 km when they enter the westerly jet. The balloons released from GAIA drift in similar directions but arrive earlier over the open sea as they are released from the Southern part of the island. A small "island effect" might show up here as the moisture profiles over the Barbados island seemed to be slightly moister than those over the nearby sea during the EUREC 4 A field campaign ( Fig. A4 in Stephan et al. (2020)). Another contribution might arise from moisture 355 transport associated with land and sea breezes as was previously evidenced in Colombo, Sri Lanka, where the land/see breezes contributed to a daytime boundary layer moistening and nighttime drying observed in radiosoundings (Ciesielski et al., 2014a).
The moisture observations at the BCO exhibit actually a day/night variation in IWV of 1.7 kg m −2 whereas the variation at GAIA is significantly smaller (1.1 kg m −2 ). Although both effects are not negligible, they are not large enough to explain the mean bias between the BCO and GAIA sondes. Instead we suspect the difference in sonde types to play a more central role.

360
The Vaisala RS41 sondes used at BCO are from the last generation of sondes and are considered to provide significantly improved temperature and humidity measurements compared to previous sonde types (e.g. Vaisala RS92) especially in cloudy conditions (Jensen et al., 2016). High confidence in the BCO soundings suggests that the GAIA soundings have a significant dry bias which may impact humidity analyses from NWP models that assimilate these observations. This seems to be the case for ERA5 and would explain the high consistency between ERA5 and GAIA IWV estimates (bias of 0.29 kg m −2 and near 365 unity slope) and the large difference of ERA5 compared to the BCO (bias of -2.33 kg m −2 and slope of 0.93). The GNSS IWV estimates are intermediate between the two sounding results. Using GNSS as the reference, we conclude on a wet bias in the BCO radiosonde data of 1.64 kg m −2 and a dry bias in the GAIA radiosonde data (-1.23 kg m −2 ) as well as in the ERA5 data (-0.77 kg m −2 ). On the other hand, using the BCO soundings data as a reference one would conclude on a dry bias in all other IWV estimates, except the HATPRO MWR after the instrumental fix on 27 January, where the MWR retrieval is slightly moist.
Since no reference water vapour measurements were made at the same time, it is difficult to establish the absolute sign of these biases. The GNSS -RS41 bias of -1.64 kg m −2 observed here represents a fractional bias of 5 % which is slightly larger the uncertainty of both systems and thus needs to be explained. From our previous experience comparing Vaisala RS92 sondes and GNSS IWV estimates (produced using the same approach as in this study) we observed slightly smaller biases of ± 0.5 to 2 % . Further investigation is needed in the case of the RS41 vs. GNSS comparisons.

375
The comparison of HATPRO IWV retrievals has been separated in two batches (before and after 27 January). Compared to GNSS, the bias before (after) is -1.28 (+2.06) kg m −2 with a slope significantly larger than one indicating that the bias increases with the amount of IWV. Similar slope values are obtained in the comparison to the other data sources, which might be attributed to the MWR training data set (Rose et al., 2005). Occurring biases could be further corrected by applying a clearsky brightness temperature offset correction based on sounding data. In terms of temporal variability, the MWR, GNSS, and 380 BCO sondes are in good agreement (standard deviation of differences ∼ 1 kg m −2 and correlation coefficient ≥ 0.98). The bias with respect to the BCO soundings is 0.48 kg m −2 (MWR too moist) during the second period, but this bias is within the known uncertainties.
ERA5 is biased low compared to both GNSS and BCO sondes with slopes lower than one, meaning that the dry bias increases at larger IWV values. The scatter plots in Fig. 6 exhibit a few outlying values which correspond to the situation when ERA5 385 underestimates the IWV peaks during the period 22-24 January (see Fig. 5). Quite surprisingly, the GAIA soundings during this period, although slightly drier, are much closer to GPS and BCO than ERA5. Inspection of the vertical humidity profiles (Fig.   7) shows that ERA5 is close to the GAIA observations in the lower 2-2.5 km but does not account for the vertical extension and the sharp drop in humidity at the top of the deep moist layer around 3-4 km on 22 January and around 5 km on 24 January. The mis-representation of the relative humidity profile is spectacular on 24 January 00 UTC. This might be due to the assimilation 390 of other data (e.g. satellite humidity sounders) that are biased low in the mid-troposphere during this event. Inspection of radar reflectivity measurements from the BCO reveals that on both dates heavy rain was occurring over an extended part of the lower troposphere (up to 3.5 km height on 22 January around 10-11 UTC and 5.0 km on 23 January around 22-23 UTC). This would explain the saturated air (RH=100 %) below 3 km for the former and below 5 km for the latter of the two BCO soundings. The GAIA soundings mimic the BCO soundings, but with a dry bias within these rainy layers. Above this layer, it seems that the 395 GAIA sondes have a moist bias, probably due to rain contamination of the humidity sensor, a problem that is corrected in the design of the newer RS41 sondes (Jensen et al., 2016). Above the moist layers, ERA5 fits well the BCO profiles.
Complementary statistics are reported in Table 5 where the data have been time-matched. The conclusions are essentially the same as discussed above, but the pair-wise biases can be more easily compared and combined.
4 Spatio-temporal distribution of IWV at the regional scale where ERA5 has a small dry bias (0.5 -1 kg m −2 ). On average over all stations, the mean difference (ERA5 -GNSS) is -0.34 ± 0.50 kg m −2 , pointing to a small dry bias of ERA5 in the region with some site-to-site variation. This result is consistent with the bias that is discussed in the previous Section based on BCO results, and also with the results of a previous study (Bosser and Bock, 2021). But further investigation is needed to explain the reason of this dry bias in ERA5 (e.g. its link with 410 the data assimilation). In terms of temporal variability, the agreement between ERA5 and GNSS is also quite good, except over Puerto Rico. There is significant spatial modulation of the magnitude of variability. Again the maximum is observed to the South over the tropical Atlantic Ocean. The average difference of standard deviation (ERA5 -GNSS) is -0.25 kg m −2 , meaning that the variability in ERA5 is slightly underestimated compared to GNSS. This may be partly due to a difference in representativeness between the GNSS point observations and the gridded reanalysis fields (Bock and Parracho, 2019) and 415 partly to some special situations where the reanalysis underestimates high IWV contents (see Fig. 5). Table S2 in the Supplement gives the comparison results for all stations. The mean and standard deviation of IWV differences (ERA5 -GNSS) have been cross-compared with the GNSS data quality diagnostics from Table S1 and no significant correlation was found. We therefore believe the main differences are not due to GNSS uncertainties but rather to differences in representativeness such as evidenced particularly over Puerto Rico. Closer inspection of the ERA5 orography for this is-   Table S1). For the islands associated with several stations (e.g. 15 for Guadeloupe and 6 for Martinique), the GNSS station with the most complete IWV series and the best data 435 quality was chosen.
The IWV times series associated with the five GNSS stations are shown in Fig. 9. They highlight substantial variability in the course of the January-February 2020 period, alternating moist (in excess of 50 kg m −2 ) and dry (below 20 kg m −2 ) episodes, sometimes within a few days as for instance in Barbados where GNSS-derived IWV was observed to decrease from 54.3 kg m −2 to 17.8 kg m −2 between 2300 UTC on 23 January and 2100 UTC on 27 January 2020 (Fig. 9d). Unsurprisingly,  (Fig. 10). On that day, St Croix was under the influence of a Fish while Barbados was surrounded by clear air. The clear-sky band in between two cloud clusters is actually part of the Fish pattern. Guadeloupe located off the southern edge of the Fish structure and was surrounded by Gravel, while Martinique and Grenade further south were bordered by Sugar. On other days such as 13 January 2020, Gravel was observed to be uniformly spread across the CA domain (Fig. 10). The classification performed for each of the five sites in January and February 2020 are provided in the Fig. 11.

465
For each of the five GNSS stations representative of the different parts of the CA, the time series of cloud classification is also shown, with Fish appearing in red, Gravel in green, Sugar in light blue, Flowers in dark blue and cloud-free conditions in white. As also observed by Stevens et al. (2020) between 2007 and 2016, Gravel-type was the dominant mode of cloud organisation over the CA during January and February 2020, with a number of cases ranging from 19 in Guadeloupe to 27 in Grenade (Table 6). Fish was the next most dominant pattern of cloud organisation with a number of cases ranging from 10 in 470 Martinique to 19 in St Croix and Guadeloupe. Sugar was the least observed cloud pattern (between 0 and 6 cases depending on the site) as also demonstrated by Stevens et al. (2020). Flowers were observed more often in the central and southern part of the CA (between 8 to 11 days), than in the northern part (on 3 days only). Finally, Cloud-free cases in between Fish and Flowers were also observed on more than 10 days during the period of interest, except in the southern part of the CA where only a few cloud-free days are observed.

475
From the time series of IWV and the cloud classification shown in Fig. 9, and if we consider the two most dominant modes, the picture emerges that Fish pattern (in red) is more systematically associated with higher IWV values than the Gravel pattern (in green). This visual impression is confirmed when computing the average IWV values associated with Fish and Gravel patterns over the five stations ( Table 6) Gravel-related IWV values range between 4.3 kg m −2 (Guadeloupe) and 7.4 kg m −2 (Grenade), and are found to be significant using a Student's t-test (see Table 7). Clear scenes over the islands of interest are also seen to be associated with rather dry In summary, based on a compilation of IWV values gathered from representative GNSS stations across the CA, we found that the environment of Fish cloud patterns to be moister than that of Flowers cloud patterns which, in turn, is moister than the environment of Gravel cloud patterns. This is consistent with the relative humidity profiles composited by Schulz et al. (2021).

495
Since Fish patterns are associated with weak winds (relative to Flowers or Gravel (Bony et al., 2020), it means that this high humidity is related to mass convergence within the column, associated with ascent. The differences in the IWV means between Fish and Gravel were assessed to be significant. Finally, the Gravel moisture environment was found to be similar to that of clear, cloud-free conditions. The moister environment associated with the Sugar cloud pattern has not been assessed because it was never prominent around the GNSS stations during the first two months of 2020 (but Sugar-like clouds occur very often 500 within mesoscale cloud patterns).
In the following, we focus on the specific period during which a large variation of IWV was observed at the GNSS station BCOS in Barbados from 54.3 kg m −2 to 17.8 kg m −2 between 2300 UTC on 23 January and 2100 UTC on 27 January 2020, associated with a transition from a Fish cloud pattern to clear, cloud-free conditions (see the cloud scene classification in Fig.   9d). For the period from 20 to 30 January, the extreme IWV values were indeed observed in Barbados (i.e. values given above), On 22 January, a large southeast-northwest oriented Fish feature is observed with MODIS to extend over a large portion of the CA, from Barbados to Guadeloupe, also covering Martinique (Fig. 13a), while Grenade is located in clearer air south of the Fish and St Croix is covered by another distinct Fish feature further north. Figure 14a shows the ERA5-derived IWV  (Fig. 12). In Barbados, a maximum IWV value of 48.6 kg m −2 was observed at 1200 UTC which is associated with a deep moist lower troposphere as observed from the radiosounding measurements made from the BCO (Figs. 7a and b). Using trajectory analyses from the LAGRANTO Lagrangian analysis tool (Sprenger and Wernli, 2015), Villiger et al. (2020) also evidenced that air parcels arriving at 1000-700 hPa above the BCO are transported from high latitudes towards the BCO by an extratropical surface cyclone/upper level 520 trough located off the US East coast. The initially dry air parcels descend from upper-levels into the boundary layer, where they experience a rapid moistening, before arriving at the BCO as anomalously humid.
On 23 January, a southeast-northwest oriented Fish feature is also observed with MODIS to extend between Barbados and Martinique. Guadeloupe is under the influence of another distinct Fish feature further north, while Sugar is observed to surround Grenade and cloud-free conditions are found over St Croix (Fig. 13b). Three distinct IWV plumes are seen in the ERA5 field, 525 one over Puerto Rico, one further east over the Netherlands Antilles almost reaching Guadeloupe and one extending from the southeast over Barbados and Martinique (Fig. 14b). Drier air masses are seen in St Croix (located between two IWV plumes to the north) and in Grenade along the southern edge of the southernmost plume. Comparison between ERA5 and GNSS IWV values suggest that the southernmost plume in ERA5 is located a bit too far south as ERA5 IWVs are underestimated in Barbados and Martinique, and slightly overestimated in Grenade. On that day also, the GNSS stations in Barbados, Martinique 530 and Guadeloupe all show IWV values in excess of 35 kg m −2 after 1200 UTC, while lower values are observed in Grenade to the south and St Croix to the north (Fig. 12). The maximum IWV value of 54.3 kg m −2 observed at the end of the day in Barbados is in good agreement with that derived from the high-resolution radiosounding performed at 2240 UTC at the BCO (Fig. 5). It exhibits a deep moist layer, with a nearly 100 % relative humidity in the first 5 km of the troposphere (Figs. 7c and d) which was also highlighted by Stephan et al. (2020), their Fig. 9. Like on 22 January, the moist air below 5 km above the 535 BCO is associated with transport from high latitudes towards the BCO by the extratropical surface cyclone/upper level trough off the US East Coast (Villiger et al., 2020).
On 24 January, the cloud scene over the CA is dominated by Fish features again, with the previously observed southern southeast-northwest oriented Fish pattern being shifted to the west, clearing most of the central CA and Barbados, while moving over Grenade (Fig. 13c). It merges over the Caribbean Sea with a southwest-northeast oriented Fish structure going 540 across the CA north of Guadeloupe. St Croix is underneath a distinct Fish at that time. As a result, Guadeloupe, Martinique and Barbados are in a mostly cloud-free air mass moving in westward behind the '<'-shaped Fish structure. As shown in Fig.   14c, the '<'-shaped structure of Fish pattern also reflects in the ERA5-related IWV field, with higher IWV values located to the south of Barbados and to the east of Grenade. The northern branch of the '<'-shaped plume appears to be located too far east in ERA5 as suggested by the slight overestimation of IWVs over Guadeloupe compared to the GNSS retrievals. Consistent with 545 the above described spatial distribution, the GNSS stations in Grenade shows the highest IWV values (reaching 50 kg m −2 ) after 1200 UTC, while lower similar values (and a similar decreasing trend) are observed in Barbados, Martinique, Guadeloupe and St Croix further north (Fig. 12).
On 25 January, Grenade and St Croix appear beneath the remains of two distinct Fish features, while cloud-free conditions are observed in the rest of the CA. The central part of the CA is bordered with Sugar to the east (Fig. 13d). The '<'-shaped 550 feature of higher IWVs is still visible in the ERA5 field further west with respect to the previous day, with most of the central part of the CA being located in a drier environment (Fig. 14d). On that day also, the northern branch of the '<'-shaped IWV plume appears to be located too far east in ERA5. The GNSS stations in Grenade shows the highest IWV values after 1200 UTC even though IWV is observed to decrease on that day as the '<'-shaped plume is moving west. IWV values in Barbados, Martinique and Guadeloupe are even smaller than the previous days as drier air masses continue to move in from the Tropical 555 north Atlantic (Fig. 12). A maximum of IWV is observed in St Croix in relationship with the presence of a Fish feature.
On 26 January, IWV values in Grenade, Barbados and Martinique continue to drop (Fig. 12)  50 kg m −2 , while larger IWV values than the previous day are also observed in Guadeloupe (Fig. 12), consistent with the spatial distribution of IWV obtained with ERA5 (Fig. 14e).
Finally, on 27 January, a well-defined Fish feature is observed in the north of the domain, covering the US Virgin Islands and St Croix (Fig. 13f), which is associated with a plume of IWV values between 40 and 50 kg m −2 (Fig. 14f). All the other stations to the south are located in the drier, mostly clear-free air mass to the east of the Fish pattern. Martinique and Guadeloupe are 565 surrounded by Sugar while the environment of Barbados and Guadeloupe is observed to be cloud-free (Figs. 11 and 13f). The GNSS station in St Croix highlights a IWV maximum of 50 kg m −2 for the second consecutive day, while in Barbados and Grenade the GNSS stations record the lowest IWV values of the period after 1200 UTC as the driest conditions are evidenced with ERA5 over the Tropical North Atlantic east of central CA (Fig. 14f). In Martinique and Guadeloupe, higher IWV values are observed (between 30 and 35 kg m −2 ) after 1200 UTC (Fig. 12) as both stations appear to be located in regions of IWV 570 gradients west of the drier air masses (Fig. 14f).
It is also worth noting that from 16 to 22 January the IWV features as modelled with ERA5 are moving westward across the CA domain (not shown). On 22 and 23 January, IWV features remain quasi-stationary, while from 24 January onward, there is a clear decoupling between the northern and southern part of the domain, with IWV features north of 14°N being advected eastward and IWV structures south of 14°N moving westward. This is the origin of the '<'-shaped IWV feature seen 575 on 24 and 25 January (Figs. 14c,d) as well as the the '<'-shaped Fish pattern observed in the MODIS images (Figs. 13c,d). The eastward motion north of 14°N on 24 and 25 January is likely related to the growing influence of an extratropical disturbance (sea surface pressure of 1005 hPa, not to be confused with the extratropical cyclone off the US East Coast) that formed north of Puerto Rico on 22 January, with its center located at 30°N, 65°W, and moved northeastward in the following days (the center of the disturbance is located at 35°N, 45°W on 25 January).

Discussion and conclusions
This paper describes the data processing streams and discussed the quality of GNSS ZTD and IWV retrievals from four operational streams run by IGN and ENSTA-B/IPGP and one research stream (GIPSY repro1), see Table 1. The two operational streams run by IGN (RGP NRT and RGP daily) provide data for stations BCON and BCOS, the two GNSS stations installed at the BCO for EUREC 4 A in October 2019, as a well as a few other sites in the French overseas territories of the Antilles.

585
The RGP NRT stream provides ZTD estimate within 45 min after the end of measurements which are generally available for assimilation into NWP models. We showed that this stream had some instabilities, especially in most recent 1-hour time slot, which are corrected in the t-2h time slot. This instability was found to be related with the processing scheme and was corrected since then by IGN. The other operational streams were shown to be in good agreement, with a small bias of 2 mm in ZTD between the IGN and ENSTA-B/IPGP solutions due to the use of two different software packages. IWV estimates are made 590 available for stations BCON and BCOS from all four operational streams since 31 October 2019. This offers the possibility to analyse the IWV time series both in near real time and on the long term. However, due to the use of various different ZTD to IWV conversion methods and auxiliary data, the uncertainty in the IWV estimates is variable and not optimal in these streams (see Table 4). It is planned in the near future to improve and homogenize the operational IWV conversion approach for all four operational streams and make these data available to the users.

595
The GIPSY repro1 research stream includes 49 GNSS stations covering the Caribbean Arc (see Fig. 1). The ZTD and IWV estimates from this stream, which is the main focus of this paper, have been analysed for the period 1 January-29 February 2020 and made available with 5 minute time sampling for scientific applications. This stream used a slightly improved processing approach and IWV conversion method and auxiliary data with reduced uncertainties. The uncertainty due to the IWV conversion is believed to be at the level of ± 0.2 kg m −2 (see the Appendix). The data have been thoroughly quality 600 checked and screened for outliers (see the Supplemental material).
The GNSS IWV estimates set have been compared to the Vaisala RS41 radiosonde measurements operated from the BCO and to the measurements from the operational radiosonde station at GAIA (see Table 5). A significant dry bias is found in the GAIA humidity observations, both with respect to the RS41 sondes (-2.9 kg m −2 ) and to the GNSS results (-1.2 kg m −2 ).
The RS41 sondes and GNSS measurements show also a systematic difference, with a mean of 1.64 kg m −2 or 5 %, where 605 the GNSS retrievals are drier than the sonde measurements. The IWV estimates from a colocated MWR agree with the RS41 results after an instrumental update on 27 January, while they exhibit a dry bias compared to GNSS and sondes before that date.
Understanding the origins of the various moisture biases requires further investigation that is beyond the scope of this paper.
However, the GNSS minus sonde bias resembles strongly the results reported by Ciesielski et al. (2014b) where GNSS IWV data were compared to corrected Vaisala RS92 radiosonde measurements at several sites in the Indian Ocean. These authors 610 found a similar dry bias in the GNSS data although good agreement was found with MWR measurements. However, the same authors also found systematic dry biases in satellite microwave data and in NWP model analyses/reanalyses compared to their corrected radiosonde data. Other studies found either dry or wet biases, or no bias at all, between GNSS, radiosondes, and other techniques (see e.g. Buehler et al. (2012) and references therein, and also Bock et al. (2007); Wang and Zhang (2008); Yoneyama et al. (2008); Ning et al. (2012)).

615
A global dry bias intrinsic in the GNSS technique is unlikely. Instead, site-specific error sources are thought to contribute to biases with variable signs and magnitudes. It has been observed that at sites with multipath, satellite visibility obstructions, and/or electromagnetic interferences, the ZTD estimates can be biased either dry or wet (Ning et al., 2011). Biases can also result from using wrong antenna models or inaccurate mapping functions. Such situations can be detected by performing cutoff tests and can be partly mitigated by using a higher cutoff angle. We checked the BCO GNSS data by reprocessing the 620 measurements using two different cutoff angles, one lower (5°) and one higher (10°) than the nominal value (7°) and found negligible impact on the mean ZTD estimates. This result excludes the hypothesis of any of those bias sources. However, GNSS is not an absolute remote sensing technique and unless the IWV estimates are compared to an adequate reference it is difficult to figure out which part of the observed bias is due to GNSS or to the RS41 sondes.
IWV estimates from the ERA5 reanalysis were also compared to GNSS data and to BCO and GAIA sonde data. It was found 625 that the IWV content from reanalysis over Barbados is overall close to the GAIA observations. Indeed, ERA5 assimilated the GAIA sondes but not the sondes from the BCO during EUREC 4 A (Irina Sandu, ECMWF, personal communication). At several occasions, ERA5 is also shown to significantly underestimate IWV peaks observed by all systems (sondes, GNSS, and MWR) by 5 to 8 kg m −2 (see Fig. 5). Two such events are documented (22 January and 23/24 January) during which a deep moist layer extended from the surface up to altitudes of 3.5 and 5 km (see Fig. 7). It was shown that ERA5 significantly underestimated the 630 moisture content in the upper part of these layers, possibly due to the assimilation of other data over the domain that were biased low. Overall, the reanalysis showed a small dry bias (0.34 ± 0.50 kg m −2 ) over the study area in comparison to the 49 GNSS stations. Although the assimilation of biased radiosonde data might be thought as a potential reason, recent experiments with the ECMWF Integrated Forecasting System show that removing all the radiosondes and dropsondes in the EUREC 4 A domain does not significantly impact the simulated humidity field (Savazzi and Sandu, personal communication). Understanding the 635 origin of the mean and occasional dry biases in ERA5 needs further investigation.
At the synoptic scale, ERA5 showed spatio-temporal variations in the IWV field over the domain which were in general in good agreement with the observations from the GNSS network. The link with the cloud organisation was studied using MODIS visible images inspired by the classification of Stevens et al. (2020). We found that the environment of Fish cloud patterns was moister than that of Flowers cloud patterns which, in turn, is moister than the environment of Gravel cloud 640 patterns. The differences in the IWV means between Fish and Gravel were assessed to be significant. Finally, the Gravel moisture environment was found to be similar to that of clear, cloud-free conditions. The moisture environment associated with the Sugar cloud pattern has not been assessed because it was hardly observed during the first two months of 2020.
These preliminary results prompt for a more systematic analysis of the cloud organisation and the lower and mid-tropospheric moisture field.

650
The ORPHEON GNSS RINEX data are provided for scientific use in the framework of the GEODATA-INSU-CNRS convention. The DOI references for the UNAVCO GNSS RINEX data are listed in Table A1  Special issue statement. This article is part of the special issue "Elucidating the role of clouds-circulation coupling in climate: datasets from the 2020 (EUREC 4 A) field campaign". It is not associated with a conference. We are grateful to AERIS, the French data and service centre for atmosphere, for providing the ERA5 reanalysis data and hosting the GNSS data. Thayer, G.: An improved equation for the radio refractive index of air, Radio Science, 9, 803-807, 1974.  from 1 January to 29 February 2020. Bottom: from 1 to 6 January 2020.    HATPRO data set has been separated in two parts, before (black symbols) and after (red) the 27 January when the system was fixed for an instrumental failure. Note that sample sizes are different between diagrams because each comparison is done with the highest temporal resolution for the period of available data between 1 January and 29 February 2020.            Table 4. Separation of bias (mean) and random errors (one standard deviation) of IWV differences into contributions from ZTD differences, ZHD differences, and K differences, for two operational data streams (RGP daily and GIPSY rapid) and the reprocessed stream (GIPSY repro1   Many authors published refractivity coefficients from experimental work performed between the 1950s and the 1970s. Smith and Weintraub (1953) compiled and averaged the early measurements, and Hasegawa and Stokesberry (1975) compiled and characterized a significantly larger number of experimental results. Thayer (1974) developed an alternative and hybrid approach which includes measurements extrapolated from optical frequencies. Bevis et al. (1994) revisited the data used by Hasegawa and Stokesberry (1975) and determined a new set of average values and associated uncertainties. Finally, Rueger (2002) pro-810 posed a new set of "best average" coefficients after reassessing the data set used by Bevis et al. (1994). While there has been a broad consensus on the value of k 1 among previous authors, Rueger's new k 1 coefficient is 0.115 % larger than previous values. The impact on ZHD would be an increase of about 2.6 mm at mean sea level (i.e. a bias in IWV of −0.4 kg m −2 ).
The impact is also significant on the determination of bending angles from GNSS radio-occultation measurements as discussed by Healy (2011). The latter author examined the origin of the increase in Rueger's k 1 and identified two obvious reasons: a 815 numerical inconsistency in the value of 0°C = 273 K instead of 273.15 K and neglecting CO 2 in the gas mixture composing the dry air in many previous studies. Healy (2011) highlights that although Rueger's estimate of k 1 appears to be more robust and defendable than the previous values, it also has one significant caveat as it does not account for non-ideal gas effects.
According to the significant work done by Rueger (2002) in re-assessing past measurements and re-evaluating the refractivity coefficients we believe that his results are the more accurate to date and will use them along with a correction for the non-ideal 820 gas effects, as suggested by Healy (2011), and an update for present CO 2 content.
Using a present day CO 2 mixing ratio of r c =408 ppm, k 1 and k 4 can be summed together to form k 1 , the refractivity constant for dry air including CO 2 , k 1 = k 1 ·(1−r c )+k 4 ·r c , to give k 1 = 77.6909 K hPa −1 . The last step is correcting for the non-ideal 825 gas effects using the compressibility factor given by Owens (1967), i.e. 1/Zd = 1.000588 for dry air at 273.15 K and 1013.25 hPa, and 1/Zw = 1.000698 for water vapour at 293.15 K and a partial pressure of 13.33 hPa (the conditions of measurements of refractivity used by Rueger (2002)). Finally, the updated refractivity coefficients become: k 1 = 77.6452 ± 0.0094 K hPa −1 k 2 = 71.2 ± 1.3 K hPa −1 830 k3 = (3.7520 ± 0.0076) · 10 5 K 2 hPa −1 where we included the uncertainties evaluated by Rueger (2002). If we assume that these uncertainties are fair estimates of the true absolute accuracy of the coefficients, the uncertainty in the IWV estimates due to k 1 , k 2 , and k 3 , would be: 0.04, 0.03, and 0.06 kg m −2 , respectively, under the average conditions of EUREC 4 A (IWV = 32 kg m −2 ). The systematic bias due to the uncertainty in the refractivity coefficients is thus extremely small. However, using the new estimates can make some systematic 835 difference compared to those used in past studies. Figure A1 compares the new refractivity coefficients to those from Thayer (1974) and Bevis et al. (1994), which are the two most widely used data sources in GNSS meteorology, and those from Rueger (2002) (without the compressibility factor correction and assuming a CO 2 mixing ratio of 375 ppm). It is seen that the k 1 from Rueger is about 0.09 K hPa −1 larger than Thayer (1974) and Bevis et al. (1994) and that the new value differs only by 0.04 K hPa −1 from the latter which represents 840 a ZHD difference of 1.22 mm and a IWV difference of -0.19 kg m −2 (the negative sign is because a larger ZHD correction decreases the IWV estimate). The new value for k 2 is in agreement with Bevis et al. (1994) and Rueger (2002), but this coefficient has anyway a small weight in the final IWV estimate. The new value for k 3 is in good agreement with Rueger (2002) and to a lesser extent with Thayer (1974) and Bevis et al. (1994) from which it differs by -0.2399 and +0.1301 ·10 5 K 2 hPa −1 , respectively, leading to a small fractional change in IWV of -0.64 % (-0.20 kg m −2 ) and +0.35 % (+0.11 kg m −2 ), 845 where the IWV biases in kg m −2 are computed assuming a mean value of IWV = 32 kg m −2 .
In conclusion, based on the differences between published refractivity coefficients and their uncertainties, we consider that the uncertainty in the absolute IWV values (i.e. the possible bias) retrieved from GNSS during EUREC 4 A due to these coefficients is at the level of ± 0.2 kg m −2 .