GPS Precipitable Water Vapor Estimations over Costa Rica: A Comparison against Atmospheric Sounding and Moderate Resolution Imaging Spectrometer (MODIS)

Abstract

The quantification of water vapor in tropical regions like Central America is necessary to estimate the influence of climate change on its distribution and the formation of precipitation. This work reports daily estimations of precipitable water vapor (PWV) using Global Positioning System (GPS) delay data over the Pacific region of Costa Rica during 2017. The GPS PWV measurements were compared against atmospheric sounding and Moderate Resolution Imaging Spectrometer (MODIS) data. When GPS PWV was calculated, relatively small biases between the mean atmospheric temperatures (T_m) from atmospheric sounding and the Bevis equation were found. The seasonal PWV fluctuations were controlled by two of the main circulation processes in Central America: the northeast trade winds and the latitudinal migration of the Intertropical Convergence Zone (ITCZ). No significant statistical differences were found for MODIS Terra during the dry season with respect GPS-based calculations (p > 0.05). A multiple linear regression model constructed based on surface meteorological variables can predict the GPS-based measurements with an average relative bias of −0.02 ± 0.19 mm/day (R² = 0.597). These first results are promising for incorporating GPS-based meteorological applications in Central America where the prevailing climatic conditions offer a unique scenario to study the influence of maritime moisture inputs on the seasonal water vapor distribution.

1. Introduction

Although it constitutes only 0.001% of the planet’s water resources, water vapor plays an important role in atmospheric processes as it is one of the major radiative gases and a dynamic element in the atmosphere. Water vapor is a useful parameter to forecast severe weather conditions and precipitation formation and is also a key factor for studying the global water cycle, changing climatic conditions, and earth-atmosphere energy exchange [1,2,3]. Overall, water vapor is essential for the development of disturbed weather and influences the planetary radiative balance. In the lower atmosphere, it controls the heat exchange during the precipitation formation and the thermal structure of the troposphere, and it is the main source for precipitation in all weather systems [3,4]. Therefore, accurate estimates of atmospheric water vapor content are needed to improve the predictability of rainfall and the understanding of and feedback in climate related processes [5,6].

A quantifiable parameter useful for studying water vapor is the precipitable (or integrated) water vapor (PWV). Precipitable water vapor mainly comprises tropospheric water vapor and the less abundant stratospheric water and can be used to analyze water vapor variability and its contributions to climate change [6]. The classical approach to gather information about PWV is using atmospheric sounding based on radiosonde profiles [7]. However, due to high costs, radiosonde networks lack spatial and temporal resolutions and, thus, provide limited information to carry out detailed studies of weather and climate. For example, radiosondes are usually launched 1–2 times per day in monitoring stations spaced several hundred kilometers from each other. In recent years, the fast development of ground-based GPS networks allows a new source of water vapor information. As atmospheric water changes the atmospheric refractivity, satellite-receiver path delays provide a unique information on the total water vapor within the troposphere and stratosphere. Therefore, GPS has become a standard technique for measuring PWV with some noticeable advantages over radiosondes. For instance, GPS can be used in all weather conditions and has low operation costs, allowing for a high temporal resolution with numerous records throughout the daytime and nighttime [8,9,10]. In Costa Rica, there are 14 Global Navigation Satellite System (GNSS) stations in operation which are associated with the Sistema de Referencia Geocéntrico para las Américas (SIRGAS) network. Eight of these GNSS stations are officially administrated by the National Institute of Geography. Although there are other GPS stations operating in the country, access to these GPS data is rather limited.

Satellite remote sensing is also a feasible method to derive the PWV distribution. The Moderate Resolution Imaging Spectroradiometer (MODIS) installed at the Terra and Aqua satellites offers spatial and temporal PWV estimations [11,12]. Despite the high spatial coverage and resolution that these satellite-based PWV products offer, there are several sources of errors in water vapor column retrievals from these remote sensing platforms. These errors are mainly linked to an uncertainty in the spectral reflectance of the surface, an uncertainty in the sensor calibration, an uncertainty in the atmospheric temperature and moisture profile, and an uncertainty in the amount of haze [13,14]. Moreover, there are two other additional limitations related to a polar orbiting satellite like MODIS: i) most areas are sampled only once per day, depending on the latitude and the configuration of the instrument, and ii) the measurements are mainly restricted to cloud-free areas (especially during daytime) as clouds are opaque in the visible and NIR spectrum [15]. Unlike satellite-based water vapor estimations, the presence of clouds and precipitation does not affect GPS observations because the liquid water contribution to the refractivity is normally small, especially outside of clouds [16]. In order to assess the performance of satellite measurements, their PWV estimates have been evaluated against other conventional techniques (e.g., GPS PWV measurements) in several regions, for instance in China, in Spain, and in Tibet [6,9,12]. Nevertheless, limited knowledge exists for the Central American Isthmus regarding the application of remote sensing for PWV measurements and how well GPS delay data compare to classical water vapor measurements made by atmospheric sounding in complex tropical mountainous regions like those found in Costa Rica.

In this study, the objectives were i) to evaluate the GPS-based estimates of PWV against PWV based on radiosonde measurements and on the MODIS satellite radiometer, ii) to estimate the influence of the main circulation patterns in Costa Rica on the PWV variability using GPS-based estimations, and iii) to identify major meteorological variables controlling PWV seasonal variations. We selected two GPS stations located in the Pacific region of Costa Rica to calculate the mean daily PWV estimates during 2017. These GPS-based estimations were then compared to PWV measurements made using radiosondes at the only atmospheric sounding site in operation in Costa Rica, located in the Central Valley of the country. We further compared data from the MODIS satellite radiometer against the GPS and radiosonde estimations over the Central Valley of Costa Rica. GPS PWV estimates were also analyzed in combination with surface meteorological data and the Hybrid Single Particle Lagrangian Integrated Trajectory (HYSPLIT) model. We expect that this work will contribute to highlighting the opportunity of incorporating GPS-based meteorological applications in Central America, which can be useful to study the influence of maritime moisture inputs from the Caribbean Sea and Pacific Ocean on the seasonal water vapor distribution.

2. Materials and Methods

2.1. Climatic Characteristics of Costa Rica

Costa Rica is located in the tropics between 8°–11°N latitude and 82°–86°W longitude (Figure 1). The climate of Costa Rica is influenced by four regional air circulation types: NE trade winds, the latitudinal migration of the Intertropical Convergence Zone (ITCZ), cold continental outbreaks, and the sporadic Caribbean cyclones [17,18,19]. Strong orographic effects are caused by a NW to SE mountain range (or cordillera) with a maximum elevation of 3820 m above sea level (m a.s.l.), which divides the country into the Caribbean and Pacific regions, each region having a distinct precipitation regime. In the Pacific region of Costa Rica, the dry season ranges from December to April and the wet season ranges from May to November. There is a secondary humidity gradient along the Pacific coast where wetness increases from north to south [20,21]. The observed cyclic deviations in the ocean-atmosphere domain can be described as “wet” and “dry” years throughout Costa Rica and are mainly linked to changes in the sea surface temperature (SST), especially the warm/cold El Niño Southern Oscillation (ENSO) episodes [17,22].

2.2. GPS and Atmospheric Sounding Data

As stated above, there are 14 GNSS stations in operation in Costa Rica, which are associated with the SIRGAS network. We selected GPS data from two of these stations to estimate PWV: one located in the Central Valley (AACR) and one situated in the northern Pacific (Liberia or LIBE). As shown in Figure 1, AACR is located in the mountainous central region of the country known as the Central Valley and LIBE is situated on the northern Pacific region. We selected these two stations based on three criteria: i) at least one station must be as near as possible to an atmospheric sounding site (AACR), ii) at least one additional station must be included in the analysis and situated in the Pacific slope of the country (the climatic region where the Central Valley is located, LIBE), and iii) for each station, a weather station must be available to register the meteorological conditions, with no significant height differences between the GPS station and the weather station.

The GPS data were processed by the National Processing GNSS Data Center. The receiver and antenna types at AACR were Topcon TPS NET-G3A and Topcon TPSCR.G3 TPSH, respectively. At LIBE, the receiver and antenna model were Leica GRX 1200 + GNSS and Leica AT504GG LEIS respectively. GPS data were processed using the software GIPSY, version 6.4 from JPL [23], using the Precise Positioning Point (PPP) method based on the precise ephemerides computed by JPL. The parameters for the satellite and receiver antenna phase center calibration were set according to the JPL products. The tropospheric model incorporated a priori hydrostatic delay (PHD, m), computed as follows:

$$\mathrm{PHD}=1.013\ast 2.27e^{\left(-0.00016\ast h\right)},$$

(1)

where h is the station height (m) above the ellipsoid. The PHD value was estimated as 0.1 m. The tropospheric gradient was estimated based Bar-Sever et al. [24]. The global mapping function (GMF) troposphere mapping functions were implemented, and an elevation cutoff angle was set at 7.5°.

The observations of AACR and LIBE stations were available from January 1st to December 31^st, 2017. To obtain PWV radiosonde estimations, we used the only atmospheric sounding site in operation in Costa Rica, namely the International Airport of San José, Costa Rica (International Civil Aviation Organization code: MROC). The radiosonde launching was carried out by the National Meteorological Institute of Costa Rica using mainly Sprenger E085 (St. Andreasberg, Germany) sounding systems. The radiosonde data were obtained from the University of Wyoming [25]. The distances and elevation differences between the GPS stations and the atmospheric sounding site are summarized in

Table 1. The location details for AACR and LIBE GPS receivers in the Central Valley and northern Pacific region of Costa Rica and for the MROC radiosonde site used in this study.

Station	AACR	LIBE	MROC
Latitude (decimal degrees)	9.9386	10.6305	9.9944
Longitude (decimal degrees)	−84.1179	−85.4380	−84.2079
Elevation (m a.s.l.)	1159	132	912
Δ distance (km) 1	11.5	152	--
Δ elevation (m a.s.l.) 1	+247	−780	--

¹ Δ the distance and Δ the elevation in relation to the MROC sounding site.

. In situ meteorological observations were measured with a Vantage Pro2 weather station (Davis Instruments, Hayward, CA, USA), with no significant height difference between the GPS stations and the weather monitoring sites.

2.3. GPS Data Processing

In general, GPS data processing is based on the physics of the atmospheric propagation delay. GPS radio waves are delayed by the ionosphere and troposphere when they travel through the atmosphere from the satellite to GPS ground-receivers. The so-called “total or zenith atmospheric delay” (or ZTD, in millimeters) of the signal emitted by a GPS satellite consists of two parts, “hydrostatic delay” or ZHD and “wet delay” or ZWD:

$$\mathrm{ZTD}=\mathrm{ZHD}+\mathrm{ZWD}$$

(2)

Overall, the ZHD is due to the effect of dry air, contributing to at least 90% of the total tropospheric delay, whereas the ZWD represents less than 10% of the signal. Therefore, the ZTD depends on the air mass between the receiver and satellite and can be expressed as a function of ground atmospheric pressure [26,27,28]:

$$\mathrm{ZHD}=\frac{0.002277{\mathrm{P}}_{\mathrm{surf}}}{1-0.00266\cos \left(2\mathsf{\theta }\right)-0.00028{\mathrm{H}}_{\mathrm{site}}},$$

(3)

where P_surf is the surface pressure (hPa), θ the geodetic latitude, and H_site represents the height (km) above the geoid [26]. Once the ZHD is calculated, ZWD is estimated by subtracting ZHD from ZTD.

Overall, the computation of ZWD using GPS delay is commonly related to the precipitable or integrated humidity along the altitudinal profile over the local atmosphere (known as GPS PWV). GPS PWV represents the total mass of water vapor in an atmospheric column with a unit area and is measured in kg/m², but it is usually reported as the height of an equivalent column of liquid water in millimeters [26,29]. We used the Bevis relationship to estimate GPS PWV using ZWD [8]:

$$\mathrm{PWV}=\mathrm{kZWD},\mathrm{k}=\frac{{10}^6}{\left(\frac{{\mathrm{k}}_3}{{\mathrm{T}}_{\mathrm{m}}}+{\mathrm{k}}_2^{\prime }\right){\mathrm{R}}_{\mathrm{v}}\mathsf{\rho }},$$

(4)

where k is a dimensionless water vapor conversion coefficient. In Equation (4), k₃ and k’₂ are empirical constants [30], Rv is the specific gas constant for water vapor, ρ is the liquid water density, and T_m is the mean temperature of the atmospheric column. To calculate T_m, we used the well-known Bevis equation [30]:

$${\mathrm{T}}_{\mathrm{m}}=70.2+0.72{\mathrm{T}}_{\mathrm{s}},$$

(5)

where T_s is the surface air temperature. In order to test the validity of this relationship for the tropical atmosphere of Costa Rica, we used the radiosonde profiles available at MROC during the study period (N = 210) to estimate the T_m values for the local atmosphere. As T_m depends both on the temperature profile and the vertical distribution of water vapor, T_m was calculated using the following equation [31,32]:

$${\mathrm{T}}_{\mathrm{m}}=\frac{{\int }_0^{\infty }\frac{{\mathrm{P}}_{\mathrm{w}}}{\mathrm{T}}\mathrm{dz}}{{\int }_0^{\infty }\frac{{\mathrm{P}}_{\mathrm{w}}}{{\mathrm{T}}^2}\mathrm{dz}},$$

(6)

where P_w is the water vapor pressure and T is the air temperature.

Based on the atmospheric sounding data available at MROC during 2017, the composite atmospheric temperature profiles depict lapse rate changes for the tropopause and troposphere for the dry and wet seasons (Figure 2A). During the dry season, the mean lapse rate was −5.0 °C/km from the ground to the tropopause level (approx. 15–20 km), whereas during the wet season, the mean lapse rate was −4.8 °C/km. Tropospheric temperature variations (up to 15 km) were similar during the wet season (range: 198–296 K, mean: 256 ± 27 K) and the dry season (range: 199–297 K, mean: 258 ± 28 K). At MROC, the mean surface temperature (T_s) varied from 293–296 K during the dry season and from 289–296 K during the wet season. The corresponding T_m values were in the range 279–289 K (dry season) and in the range 277–302 K (wet season). When we fitted a straight line to the T_m data to obtain a T_m–T_s relationship for MROC (black line in Figure 2B), we found a poor Spearman’s correlation between T_m and T_s for our data (r = 0.0257, p > 0.05). Therefore, we compared the relative bias of the Bevis equation (Equation 5, plotted as a red line in Figure 2B) to the T_m calculations from the radiosonde data. The mean relative bias using the Bevis equation was −0.009 ± 0.008 for the dry season estimations (N = 70) and 0.004 ± 0.009 for the wet season calculation (N = 140). We also calculated a RMSE of 3.50 K for the dry season, with a RMSE of 2.72 K for the wet season. The estimated relative biases are equivalent to the mean error values of −2.6 ± 2.4 K (dry season) and 1.1 ± 2.5 K (wet season). In terms of GPS PWV, the mean error associated to these T_m deviations were in the range of −0.2 and 0.4 mm. Therefore, as the relative biases for the Bevis equation are smaller than the estimated precision for the mean daily PWV calculations, we decided to apply the Bevis equation to estimate T_m in the calculations of PWV at our study sites.

In this work, we report mean daily PWV estimations based on hourly calculated ZTD and ZHD values at AACR and LIBE. We decided to carry out our analysis on a daily basis because there is limited atmospheric sounding data for Costa Rica, with only one radiosonde station in operation. Therefore, sounding data can be only considered representative of the average daily atmospheric conditions. The average precisions associated with the mean daily PWV calculations are 1.3 mm and 1.1 mm for AACR and LIBE, respectively.

2.4. MODIS Data

MODIS is a radiometer on board the Terra (launched in 1999) and Aqua (launched in 2002) satellite platforms. The MODIS instruments on the Terra and Aqua image the same area on Earth approximately three hours apart, observing the entire Earth’s surface every 1 to 2 days. Terra’s sun-synchronous, near-polar circular orbit passes the equator from north to south (descending node), whereas Aqua’s sun-synchronous, near-polar circular orbit crosses the equator from south to north (ascending node). The water vapor remote sensing method is based on detecting the absorption by the water vapor of the reflected solar radiation after it has transferred down to the surface and back up through the atmosphere. The total vertical amount of water vapor can be estimated from a comparison between the reflected solar radiation in the absorption channel and the reflected solar radiation in nearby non-absorption channels. The solar radiation between 0.86 and 1.24 µm on the sun-surface-sensor path is subjected to atmospheric water vapor absorption but also to atmospheric aerosol scattering and surface reflection. Therefore, in order to estimate column water vapor from measurements of the solar radiation reflected by the surface, the absorption and scattering properties of the atmosphere and the surface near 1 µm must be considered [33]. The PWV products are derived from infrared (IR) and near-infrared (NIR) measurements. NIR bands are used for daytime measurements (solar radiation reflected by Earth + atmosphere), and IR bands are used during nighttime conditions (radiation emitted by Earth + atmosphere). If clouds are present, other channels in the range of the 0.8−2.5μm region can be used in order to estimate the absorption due to water vapor above and within clouds [13,14].

Among the available MODIS products, the Level-3 MODIS Atmosphere Daily Global Product contains roughly 600 statistical datasets that are derived from approximate 80 scientific parameters from four Level-2 MODIS Atmosphere Products: Aerosol, Water Vapor, Cloud, and Atmosphere Profile. There are two MODIS Daily Global data product files: MOD08_D3, containing data collected from the Terra platform, and MYD08_D3, containing data collected from the Aqua platform [11]. In this study, the level-3 MODIS Terra and Aqua products of the daily mean (MOD08_D3 and MYD08_D3, respectively) global grid with a spatial resolution of (1° × 1°) were used to conduct the GPS PWV comparison during 2017. We selected the square region of 30 × 30 km dimensions centered on the MROC sounding site to calculate the satellite PWV estimations [34]. MODIS data estimates were calculated as area-averaged values and were processed using the Earth Observing System Data and Information System (EOSDIS) Giovanni website [35]. A total of 299 and 267 PWV Aqua and Terra satellite estimations were available from the MODIS data product, respectively, for the study period. The typical uncertainty of the MODIS PWV estimations is approximately 5–10% [13].

2.5. HYSPLIT Air Mass back Trajectory Analysis

Air mass back trajectory analyses were conducted using the HYSPLIT Lagrangian model developed by the Air Resources Laboratory (ARL) of the National Oceanic and Atmospheric Administration (NOAA, USA) [36,37]. Representative air parcel trajectories were estimated 72 h backwards in time due to the nearness of the Caribbean Sea and the Pacific Ocean. Each trajectory was calculated using NOAA’s meteorological data files (GDAS, global data assimilation system: 2006–present; 0.5° resolution) as input for the HYSPLIT model [38]. The ending altitude of air masses was set to the mean elevation of the Central Valley of Costa Rica (approx. 1,100 m a.s.l.). Trajectory analysis ending times at the Central Valley (AACR) were set to 12:00 UTC, which corresponds to a local time of 06:00 a.m. in Costa Rica. Given the estimated residence time of water in the atmosphere, ranging from around 4–10 days, weekly (N = 52, Figure 3), air mass back trajectories were calculated [38]. The ending dates for the trajectory analysis were set on Sunday of every week. These air masses were classified into two main groups, dry season (January–April) and wet season (May–December), to compare and identify the moisture transport pathways followed by the air masses that arrived at the Central Valley of Costa Rica.

2.6. Statistical Analysis

A Kruskal–Wallis non-parametric one-way analysis of variance on ranks was used to investigate if the GPS PWV stochastically dominates the other PWV estimations (i.e., atmospheric sounding and MODIS Aqua and Terra) during the dry and wet season [39]. A pairwise multiple comparison procedure was applied using Dunn’s method for those groups having a significant difference in PWV in order to isolate the stochastic dominance of the group or groups that differ from the others [40].

We also applied a multiple linear regression (MLR) model using surface meteorological data in order to identify the major variables controlling the PWV values in the Central Valley (AACR GPS station).

The cumulative annual precipitation for AACR during 2017 was 2586 mm with a mean daily precipitation during the dry season of 2 mm (range: 0 mm–54 mm) and 11 mm (range: 0 mm–85 mm) during the wet season. In LIBE, the corresponding cumulative annual precipitation was 2161 mm with a mean daily precipitation during the dry season of 1 mm (range: 0 mm–49 mm) and 9 mm (range: 0 mm–247 mm) during the wet season. At both sites, maximum daily precipitation values were recorded at the end of October (Figure 4A). Despite the differences in the cumulative annual precipitation, the relative humidity and air temperature variations were similar in the two regions. The mean relative humidity was 80.8 ± 5.6% (range: 39.0%–92.6%) and 82.4 ± 8.0% (range: 34.6%–93.6%) in AACR and LIBE, respectively (Figure 4B). The mean air temperatures were 21.0 ± 0.9 °C (range: 16.2 °C–23.2 °C) and 26.7 ± 1.3 °C (range: 22.5 °C–30.5 °C) in AACR and LIBE, respectively (Figure 4C). However, on a seasonal basis, the mean daily air temperatures and mean daily relative humidity were 6.3 °C and 2.1% greater in LIBE than in AACR during the dry season, respectively. During the wet season, the corresponding values were 5.2 °C and 1.2% greater in LIBE than in AACR, in that order.

3. Results

3.1. Seasonal Variations of GPS PWV in AACR and LIBE

During the study period, the HYSPLIT trajectory analyses identified that air masses reaching Costa Rica predominantly came from the southeastern Caribbean Sea with a less frequent contribution from the Pacific Ocean (Figure 3). Overall, the wind direction and speed in Costa Rica are mostly influenced by the seasonal migration of the ITCZ. Thus, during the dry season (December–April) when the ITCZ is located south of Costa Rica, air mass trajectories were associated with the influence of the NE trade winds. During the wet season (May–November), NE trade winds were weaker due to the passage of the ITCZ over Costa Rica, and cross-equatorial winds from the southern hemisphere transported moisture from the Pacific Ocean to the Central American Isthmus. This moisture transport pattern controlled the precipitation regimes observed at the Central Valley (AACR, Figure 4A) and the northern Pacific region of Costa Rica (LIBE, Figure 4A). During the study period, approx. 23% (N = 12) of the air masses arrived from the Pacific Ocean and the rest (approx. 77%, N = 40) came from the Caribbean Sea. Air masses arriving from the Pacific Ocean and the Caribbean Sea predominantly traveled over the eastern Pacific Ocean and the central and southern Caribbean Sea basins, respectively. No significant differences were found between the mean sea levels of the air masses that reached the Central Valley in the dry and wet seasons, with typical mean sea levels of 1500 m to 2000 m.

Seasonal GPS PWV variations were clearly defined at AACR and LIBE (Figure 5A). During the dry season, GPS PWV values varied from 14.8 mm to 40.9 mm in AACR (mean value: 27.6 ± 6.3 mm), whereas in LIBE, the variation was in the range 20.2 mm–55.5 mm (mean value: 36.9 ± 7.6 mm). We observed an increment in the GPS PWV values at both sites at the end of April and at the beginning of May that coincides with the onset of the wet season in Costa Rica, namely the passage of the ITCZ. During the wet season, at AACR, the GPS PWV estimates ranged from 24.3 mm to 46.2 (mean value: 39.7 ± 3.7 mm), and at LIBE, these GPS PWV values fluctuated from 31.5 mm to 62.6 mm (mean value: 54.1 ± 5.2 mm). At the end of November, when the transition wet-to-dry season began, we registered a decrease in the GPS PWV estimations related to the beginning of the dry season and the influence of the NE trade winds. Overall, the GPS PWV values were greater at LIBE than at AACR due to the elevation difference between the GPS stations (Δ1,027 m a.s.l.). For example, the mean differences between the estimations for AACR and LIBE were −9.5 ± 4.5 mm in the dry season and −14.4 ± 3.0 mm in the wet season. As shown in Figure 5A, these observed differences in the GPS PWV measurements at the GPS stations were more evident during the wet season when the ITCZ predominantly influenced the air circulation over Costa Rica. During the dry season, on the other hand, the differences between the GPS PWV values for AACR and LIBE were relatively more difficult to separate. However, although the GPS stations are situated approx. 160 km from each other (one in the Central Valley, AACR, and the other one in the northern Pacific region of Costa Rica, LIBE), we found a good Spearman’s correlation (r = 0.929, p > 0.001) between the GPS PWV values estimated for AACR and the corresponding estimations calculated for LIBE (Figure 5B). The best performing linear regression model shown in Figure 5B overall explained 86.9% of the variance for the GPS PWV estimates calculated for LIBE using the GPS PWV values at AACR. Overall, this finding confirms that the PWV variations at both sites are controlled by the climatic conditions of the Pacific slope which is also reflected in the precipitation patterns and air temperature/relative humidity variations shown in Figure 4A–C. This is an important result that demonstrates the applicability of PWV to monitor changes in the hydrometeorological conditions at regions that share similar climatic conditions. Additionally, our HYSPLIT analysis is able to identify the seasonal PWV variations at AACR. For instance, air masses arriving from the Pacific Ocean between May and October are associated with high PWV estimations with values between 39 and 44 mm/day. These values are practically equal to or greater than the 75th percentile of our data set (41 mm/day). In turn, air masses coming from the Caribbean Sea were associated with greater variations in the PWV estimations registered between November and April but also to smaller PWV estimations (up to 14 mm/day, Figure 5A).

3.2. GPS PWV Comparison to Other Estimations Methods and MRL Analysis

GPS PWV observations at the Central Valley of Costa Rica (AACR) compared well to the atmospheric sounding measurements during the dry and wet season but only to the MODIS Terra estimations during the dry season. As shown in Figure 6A, GPS PWV observations followed the seasonal variations registered using the radiosonde data. The best performing satellite-based estimations were those retrieved from the MODIS Terra, which also followed the seasonal variations in the GPS and radiosonde PWV observations. Unlike MODIS Terra, MODIS Aqua PWV estimations showed a systematic positive bias with respect the GPS PWV values and the radiosonde data. To better identify the seasonal differences found after applying these PWV estimation methods, we split our data set into two groups: dry season and wet season estimations (Figure 6B,C). For the dry season, the GPS PWV median value (26.5 mm) was not significantly different from the radiosonde PWV median value and the MODIS Terra PWV median estimation (27.0 mm and 25.8 mm, respectively; p > 0.05). However, it was significantly different from the median value estimated using the MODIS Aqua PWV values (29.7 mm, p > 0.001). In turn, for the wet season, the GPS PWV median value (40.3 mm) was significantly different from the MODIS Terra and MODIS Aqua PWV median estimations (36.0 mm and 51.4 mm, respectively; p < 0.001) but not significantly different from the radiosonde PWV median value (41.4 mm, p > 0.05). The mean relative biases for MODIS Aqua PWV and MODIS Terra PWV were also calculated using the GPS PWV as a reference. During the dry season, these values corresponded to 0.16 ± 0.24 mm and 0.02 ± 0.30 mm, respectively, and were equivalent to RMSE values of 7.43 mm and 7.21 mm, in that order. During the wet season, the mean relative biases were 0.30 ± 0.24 mm and −0.06 ± 0.19 mm, respectively, corresponding to 15.2 mm and 8.05 mm, respectively.

Using the available surface meteorological data at AACR, we conducted a multiple linear regression (MLR) analysis to identify the major drivers controlling the seasonal variability of GPS PWV measurements made at the Central Valley of Costa Rica (AACR, Figure 7A). The best performing MLR model was calculated as

$$\text{GPS PWV}=4.257\left(\mathrm{T}\right)+0.355\left(\mathrm{RH}\right)-0.0486\left(\mathrm{FLUX}\right)-0.257\left(\mathrm{P}\right)-999.125,{\mathrm{R}}^2=0.597$$

(7)

where T is the mean daily air temperature (K), RH is the mean daily relative humidity (%), FLUX is the mean daily downward solar radiation flux (W/m²), and P is the mean daily atmospheric pressure (hPa). The mean relative bias associated with the estimations of GPS PWV values using this model during the dry season was 0.10 ± 0.25 mm, whereas for the wet season the mean relative bias was −0.04 ± 0.09 mm. The corresponding RMSE estimated for the dry and wet seasons were 6.09 mm and 4.02 mm, respectively. When this model was applied to estimate the sounding PWV measurements, the mean relative bias during the dry season and wet season were 0.12 ± 0.27 mm and −0.06 ± 0.10 mm, respectively. The RMSE values calculated for the dry and wet season were 6.83 mm and 4.72 mm, respectively. As shown in Figure 7B, the correlation between the GPS PWV data at MROC and the corresponding PWV values estimated from the MLR model was better for the values between 30 and 45 mm. As these values were mostly registered during the wet season, it seems that our MLR model performs better when the atmospheric conditions in the Central Valley are controlled by the seasonal migration of the ITCZ and worse during the less stable atmospheric conditions linked to the influence of NE trade winds.

4. Discussion

Because of the isthmian geographical environment of Costa Rica (with only a Pacific-to-Caribbean coast distance of approx. 400 km), the good correlation between the GPS PWV estimates at AACR and LIBE was an expected result of our analysis, as both sites are located on the Pacific slope and share similar climatological features (e.g., analogous precipitation patterns, Figure 4A). Moreover, we also confirm the good agreement between the GPS PWV calculations and the radiosonde-based measurements reported by others [6,9,26,27,41]. For example, Spearman’s correlation coefficients for the GPS PWV and the radiosonde-based calculations were 0.913 (p < 0.001) and 0.902 (p < 0.001) during the dry and wet season, respectively. With respect to the MODIS satellite estimations of PWV, our analysis yielded significant biases depending on the season of the year, which are related to the annual cycle of water vapor, the NE trade winds influence, and the ITCZ activity over the Central American Isthmus. For instance, only the MODIS Terra PWV estimations recorded during the dry season were not significant biased with respect the GPS PWV calculations. However, we also found good correlations between the MODIS Aqua and the GPS-based calculations, with Spearman’s correlation coefficients of 0.735 (p < 0.001) and 0.621 (p < 0.001) for the dry and wet season, respectively. The dry and wet season MODIS Terra Spearman’s correlation coefficients were 0.591 (p < 0.001) and 0.368 (p < 0.001), respectively. These correlation values were similar to those reported over different regions of the Iberian Peninsula, including island environments like Mallorca and several coastal sites [9]. Such relationships also allowed a further adjustment of the data to fit the observations by adopting a spatial bias (error) correction method like the one applied to precipitation data [21,42,43]. As mentioned above, due to the location of Costa Rica on the narrow land-bridge of Central America, the MODIS near-infrared water vapor retrieval algorithm could be greatly affected and the derived column water vapor values over coastal or water areas may vary significantly due to the lower signal-to-noise ratios of the measured spectra [13]. This effect on the MODIS retrieval algorithm was particularly evident in the MODIS Aqua PWV estimation which was somewhat high during both the dry and wet season. MODIS Terra also showed deviations but were related to an underestimation during the wet season which can be related to the so-called shielding effect (i.e., clouds are probably occulting water vapor underneath them) [13,44]. The differences between the MODIS Aqua and MODIS Terra estimations could be attributed to their different passing times over the Central American Isthmus (MODIS Aqua crosses the equator in the afternoon, whereas MODIS Terra does it in the morning) and to the use of different radiations to estimate the water vapor during the day and night. The MODIS Aqua estimations could be higher than the corresponding MODIS Terra values because the algorithm uses IR radiation during nighttime, which could be affected by the presence of clouds with water vapor, leading to overestimations. Overall, our HYSPLIT air mass trajectory analysis is consistent with the prevailing regional moisture transport mechanism during the dry season, the Caribbean Low Level Jet (CLLJ). During the wet season, in turn, there is an intensification in the genesis and development of deep convection systems on the Pacific coast of Costa Rica which is generally is associated with the presence of the “Chorro del Occidente Colombiano” or CHOCO jet [45]. These circulation patterns produced the two rainfall maxima observed on the Pacific slope, one in June and one in September, which were interrupted by a relative minimum between July–August, known as the Midsummer Drought, due to the intensification of trade winds over the Caribbean Sea [46]. The radiosonde data were also useful to validate the atmospheric conditions controlling the GPS PWV estimations. First, the composite temperature profiles calculated using the radiosonde data are in agreement with the previously reported structure of the upper troposphere and lower stratosphere over Costa Rica [47]. As shown in Figure 2A, the temperatures in both the dry and wet season are roughly the same at 25 km, but below this level (e.g., 15–20 km), the boreal winter (December to April) temperature profile is colder than in boreal summer (from May to October). This finding was previously attributed to the influence of wave-induced vertical motions across strong vertical gradients, the source variability in the air masses arriving at Costa Rica (e.g., tropical western Pacific or midlatitudes) resulting from horizontal transport and changes induced along parcel paths due to physical and/or chemical processes [47,48,49]. Secondly, despite these differences in the thermal structure of the tropical atmosphere of Costa Rica, the T_m calculations using the Bevis equation showed small differences with respect to the corresponding calculations based on radiosonde data. This finding also agrees with the calculations made in Algeria and Argentina where Namaoui et al. and Fernández et al. estimated the uncertainty of the T_m values and found that variations up to 15K produced small differences in the final estimation of GPS PWV, which did not exceed 2 mm [27,40]. Thirdly, the poor correlation observed between T_s and T_m at MROC deserves further discussion. It is generally considered that the most accurate method to obtain T_m is by using both temperature and humidity profiles from radiosonde data [42,50]. Therefore, we have confidence that our T_m estimations are good approximations of the temperature profiles over the Central Valley of Costa Rica. A possible explanation for this finding is the mountainous and isthmian characteristics of the Costa Rica territory. The atmospheric sounding site is located on the southwestern area of the Central Valley. From this site, the distance to the Pacific coast is only 55–60 km. Radiosondes typically head in that direction after they are launched. Therefore, it seems that the atmospheric profiles estimated from MROC are representative not only of atmospheric conditions over the Central Valley but also of the Pacific coast of Costa Rica. There is also a limitation regarding the time of day when the sounding is performed. At MROC, atmospheric sounding is only done once a day, typically at 12Z or 7:00 a.m. Central American time. Therefore, T_m estimations with respect to T_s represent only the atmospheric conditions prevailing during the morning when constant surface temperatures are observed (approx. 294K ± 1K; Figure 2B). In consonance with these results, it was decided to rely on the Bevis equation to estimate the hourly T_m values for the Pacific slope of Costa Rica as this model has been extensively applied to estimate weighted atmospheric temperatures in several regions.

The MLR model estimated for GPS PWV data at AACR clearly matched the seasonal changes correctly, simulating smaller GPS PWV values during the dry season (from December to April) and much greater values during the wet season (from May to November). The best-performing and most parsimonious model included, as expected, near-surface (T and RH, Equation 7) and vertical atmospheric predictor variables (FLUX and P, Equation 6). The GPS PWV values were positively correlated with air temperature (T) and relative humidity (RH), with Spearman’s correlation coefficients of 0.210 and 0.426 (p < 0.001), respectively, and were negatively correlated with solar radiation (FLUX) and air pressure (P), with Spearman’s correlation coefficients of −0.360 and −0.175 (p < 0.001), respectively. These correlation results can be considered physically meaningful and can explain the overall model performance, although it is worth mentioning that, like the MODIS satellite estimations, it suffers from seasonal biases, specially during the dry season when the small PWV measured by the ground GPS receivers were not reproduced. This worse performance of the model during the dry season compared to the wet season was also evident after biases and RMSE values were additionally estimated using the sounding PWV measurements.

5. Conclusions

The combined analysis of PWV using GPS-based estimations, MODIS satellite products, and atmospheric sounding in the Pacific region of Costa Rica provides the first comparison between different water vapor calculation techniques for the Central American region. The evaluation of GPS-based estimates of PWV confirms the good performance of these estimations in comparison to the traditional and standard technique based on radiosondes, with no significant differences during the dry and wet seasons. These first results demonstrate the feasibility of incorporating GPS-based meteorological applications in order to improve the study of moisture inputs on the seasonal water vapor distribution in Central America. However, the performed evaluation identified significant biases between the GPS PWV estimates and the MODIS Aqua PWV estimations under both dry and wet season conditions and only the MODIS Terra PWV estimations recorded during the dry season were not significantly biased relating to the GPS PWV calculations. These results open the opportunity to evaluate other satellite products that provide higher spatial and temporal resolutions in order to give better insights into the causes of disagreements. Our analysis was also able to identify the influence of the main circulation patterns in Costa Rica, namely the trade wind regime and the ITCZ passage on PWV variability, which resulted in the relatively greater variability of the smaller PWV values during the dry season in comparison to the relatively smaller variability of the greater PWV values observed during the wet season. The influence of these moisture transport patterns was identified using the HYSPLIT analysis done for the Central Valley of Costa Rica. The multiple linear regression model successfully applied to this region can simulate the seasonal PWV variations using major meteorological variables, namely the mean daily air temperature, the mean daily relative humidity, the mean daily downward solar radiation flux, and the mean daily atmospheric pressure. We consider that a further analysis based on hourly GPS data could better analyze these relations between water vapor and HYSPLIT calculations and could refine the mathematical modeling presented in this work.