**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–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–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 o fficially 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 o ffers spatial and temporal PWV estimations [11,12]. Despite the high spatial coverage and resolution that these satellite-based PWV products o ffer, 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 a ffect 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–19]. Strong orographic e ffects 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].

**Figure 1.** The location of the GPS stations (Liberia, LIBE and Central Valley, AACR, green circles) in Costa Rica and the atmospheric sounding site at the San José International Airport (International Civil Aviation Organization code: MROC, red triangle): The AACR and the MROC sounding site are situated in the central mountainous region of Costa Rica (Central Valley), whereas LIBE is located on the dry corridor of Central America (northern Pacific of Costa Rica).

#### *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 di fferences 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:

$$\text{PHID} = 1.013 \ast 2.27 e^{(-0.00016 \ast h)},\tag{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 cuto ff angle was set at 7.5◦.

The observations of AACR and LIBE stations were available from January 1st to December 31st, 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 di fferences between the GPS stations and the atmospheric sounding site are summarized in Table 1. In situ meteorological observations were measured with a Vantage Pro2 weather station (Davis Instruments, Hayward, CA, USA), with no significant height di fference between the GPS stations and the weather monitoring sites.


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

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

#### *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:

$$\text{ZTD} = \text{ZHD} + \text{ZWD} \tag{2}$$

*Climate* **2019**, *7*, 63

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–28]:

$$\text{ZHD} = \frac{0.002277 \,\text{P}\_{\text{surf}}}{1 - 0.00266 \cos(2.0) - 0.00028 \text{H}\_{\text{site}}},\tag{3}$$

where Psurf is the surface pressure (hPa), θ the geodetic latitude, and Hsite 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/m2, 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]:

$$\text{PWV} = \text{kZWD} \text{, k} = \frac{10^6}{\left(\frac{\text{k}\_3}{\text{T}\_m} + \text{k}\_2'\right) \text{R}\_\text{V} \rho} \tag{4}$$

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

$$\rm{T\_m} = 70.2 + 0.72T\_{sr} \tag{5}$$

where Ts 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 Tm values for the local atmosphere. As Tm depends both on the temperature profile and the vertical distribution of water vapor, Tm was calculated using the following equation [31,32]:

$$T\_{\rm m} = \frac{\int\_0^\infty \frac{\mathbf{P\_w}}{T} d\mathbf{z}}{\int\_0^\infty \frac{\mathbf{P\_w}}{T^2} d\mathbf{z}} \,\tag{6}$$

where Pw 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 (Ts) varied from 293–296 K during the dry season and from 289–296 K during the wet season. The corresponding Tm 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 Tm data to obtain a Tm–Ts relationship for MROC (black line in Figure 2B), we found a poor Spearman's correlation between Tm and Ts 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 Tm 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 Tm 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 Tm in the calculations of PWV at our study sites.

**Figure 2.** (**A**) A composite atmospheric temperature (K) profile constructed using radiosonde measurements at the MROC sounding site during the dry season (N = 70, red dots) and wet season (N = 140, blue dots). (**B**) The surface temperature (Ts, K) vs. mean temperature of the atmospheric column (Tm, K) used to calculate the Tm–Ts relationship for the Central Valley of Costa Rica (black line) and the Bevis equation [8] (red line).

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μ<sup>m</sup> 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.
