Tracing Water Sources and Fluxes in a Dynamic Tropical Environment: From Observations to Modeling

Tropical regions cover approximately 36% of the Earth’s landmass. These regions are home to 40% of the world’s population, which is projected to increase to over 50% by 2030 under a remarkable climate variability scenario often exacerbated by El Niño Southern Oscillation (ENSO) and other climate teleconnections. In the tropics, ecohydrological conditions are typically under the influence of complex land-ocean-atmosphere interactions that produce a dynamic cycling of mass and energy reflected in a clear partition of water fluxes. Here, we present a review of 7 years of a concerted and continuous water stable isotope monitoring across Costa Rica, including key insights learned, main methodological advances and limitations (both in experimental designs and data analysis), potential data gaps, and future research opportunities with a humid tropical perspective. The uniqueness of the geographic location of Costa Rica within the mountainous Central America Isthmus, receiving moisture inputs from the Caribbean Sea (windward) and the Pacific Ocean (complex leeward topography), and experiencing strong ENSO events, poses a clear advantage for the use of isotopic variations to underpin key drivers in ecohydrological responses. In a sequential approach, isotopic variations are analyzed from moisture transport, rainfall generation, and groundwater/surface connectivity to Bayesian and rainfall-runoff modeling. The overarching goal of this review is to provide a robust humid tropical example with a progressive escalation from common water isotope observations to more complex modeling outputs and applications to enhance water resource management in the tropics.


INTRODUCTION
Global warming scenarios and recent evidence suggest a large potential for spatio-temporal changes in the surface water balance and the consequent intensification of the global water cycle (Giorgi, 2006;Lau and Wu, 2007;Diffenbaugh and Giorgi, 2012;Huntington, 2016;Rifai et al., 2019;Stevens et al., 2019). Both observed and projected changes are characterized by abrupt shifts between aboveaverage rainfall and prolonged drought periods often exacerbated by El Niño-Southern Oscillation (ENSO) teleconnections (Khouakhi et al., 2017;Stephens et al., 2018;Boers et al., 2019;Duque-Villegas et al., 2019;Roderick et al., 2019). Moreover, predicted scenarios suggest warm ENSOlike conditions (Vecchi and Soden, 2007;Zheng et al., 2019) and evidence of the role of ENSO in the enhancement of the carbon cycle (Knutson and Manabe, 1994;Kim et al., 2017) rises concern over a potential acceleration of warming. The latter underscores an urgent need to better understand short and long-term regional and local ecohydrological responses González-Zeas et al., 2019) in both arid and humid tropical biomes (Kimble and Stewart, 2019).
Stable isotopes studies in the tropics have frequently been constrained by numerous limitations, such as site remoteness, accessibility disruptions during extreme weather, interrupted monitoring (i.e., lack of consistent long-term time series), insufficient local funding, and a variety of socio-economic and political challenges (Landholm et al., 2019;Lee et al., 2019). The latter, in large degree, limits the reliability and robustness of temporal and spatial isotopic assessments (Sánchez- Murillo and Durán-Quesada, 2019). Beyond the lack of scientific information, factors such as rapid population growth, abrupt land use changes (i.e., agricultural, urban development, and mining land expansion), and a large dependency on groundwater exploitation are severely affecting water security and sustainability (Barrow, 2016;Hund et al., 2018;Roberts et al., 2018) across the tropics. Contrasting to temperate environments, whereby hydrological responses are well-constrained by seasonal modes (i.e. sinusoidal patterns; Allen et al., 2018), across tropical regions a nearly-continuous combination of strong orographic effects, rainfall modes (i.e., ITCZ, MJO, monsoon, and tropical cyclones), local and regional moisture recycling (i.e., ocean and land), and intense evapotranspiration fluxes play an important role in influencing the isotopic ratios, which in turn poses an advantage for modeling applications (i.e., large isotopic differences between water reservoirs throughout the year).
The recognition of the essential value of long-term experimental data Riveros-Iregui et al., 2018) relies on the powerful combination of hydrometric and tracer data to address current and future water management challenges. This manuscript, focused on stable isotope hydrology, presents the case of Costa Rica as an example of humid tropical region and reviews seven years of continuous integrated isotopic and hydrometeorological monitoring, including key insights learned, main methodological advances and limitations (both in experimental designs and data analysis), potential data gaps, and future research opportunities with a humid tropical perspective.
The structure of this review includes five sections presented in a sequential approach. First, rainfall generation and stable isotope variability are presented from conventional high frequency observations (i.e., monitoring network), and then analyzed in the light of the current state of rainfall type convective/stratiform ratio rationale (Aggarwal et al., 2016;Munksgaard et al., 2019). Local and regional climatic features are discussed, including the effects of the Inter-tropical Convergence Zone (ITCZ), ENSO, and tropical storms on rainfall isotopic seasonality and interannual variability . Second, surface water and groundwater connectivity is examined from seasonal spatially-distributed samplings and second order variables. Similarly, new soil water isotope data from rainforests is presented and discussed in the context of the ecohydrological separation hypothesis (McDonnell, 2014;Evaristo et al., 2015) which proposes that soil water is separated into two distinct pools: a mobile pool of water that feeds groundwater and streams and an immobile pool that is tightly bound to the soil matrix and is used by plants. Recharge is analyzed from isotopic lapse rates and main percolation elevations, as well as from precipitation to groundwater isotopic ratios (hereafter P/GW; Jasechko and Taylor, 2015). Third, the influence of magmatic-derived fluids in volcanic mountainous aquifers is reviewed from an andesitic water contribution perspective (i.e., waters associated with more silicic andesitic magmatism typical of mature volcanic arc systems; Taran et al., 1989;Giggenbach, 1992) based on a Bayesian mixing model approach. Fourth, advances in tracer-aided spatiallydistributed models are discussed from a flow path, runoff generation, and water age standpoint (Birkel and Soulsby, 2015;Alo-Alo et al., 2017). Finally, a summary of concluding remarks based on current data and knowledge gaps, methods, experimental designs, data collection, and future research opportunities is presented.
Our work proposes a potential route for pilot tropical isotopic studies whereby conceptual models in atmospheric, hydrogeological and ecohydrological studies are often based on non-systematic monitoring efforts. In addition, this review introduces the ongoing assimilation of multiple isotopic endmembers into tracer-aided and Bayesian mixing models and applications to improve water resources management in a changing tropical climate.
RAINFALL GENERATION: LOCAL, REGIONAL, AND GLOBAL FEATURES

Tropical Isotopic Lapse Rates
Despite the increasing computational power capabilities for data analysis and available satellite information, a significant lack of continuous and high frequency (i.e., weekly, daily or event-based) stable isotope monitoring of rainfall (Bailey et al., 2017) persists across the tropics. Water vapor isotopic measurements are still scarce due to the intrinsic challenges of upper air sampling and monitoring, while available water vapor/fog isotopic records are biased to the northern mid-latitudes (Wei et al., 2019;Zhang et al., 2020). These data are crucial to cross-validate Global Circulation Models (GCMs) and satellite measurements of current changes in atmospheric circulation, particularly during ENSO modes . Moreover, the implementation of isotope enabled GCMs, as well as Regional Circulation Models (RCMs), allows a significant improvement in the climate models with the incorporation of the physical basis of phase change processes at various scales (Yoshimura et al., 2016;Pfahl et al., 2012). Global initiatives (see Munksgaard et al., 2019) have increased spatial and temporal sampling resolution with the aim of improving our current understanding of the links between stable isotopes in rainfall and the atmospheric factors governing in tropical regions.
One such atmospheric factor where isotopes can provide additional understanding is the Rayleigh adiabatic cooling process. This process describes the preferential condensation and removal of isotopically heavier molecules during water vapor uplift (from a single source; assuming equilibrium under vapor and liquid water) across elevation gradients (Dansgaard, 1964;Jouzel and Merlivat, 1984). Figure 1 summarizes mainly tropical and subtropical isotopic lapse rates reported as δ 18 O in ‰ km −1 FIGURE 1 | Isotopic variations across the tropics. The upper panel shows isotopic lapse rates across the tropics based on published data (see Supplementary  Table S1 , 2020). In general, the δ-elevation relationship is mostly affected by the temperature vertical gradient (lapse rate; dT/dz) and the initial relative humidity (h) of the ascending air masses (Gonfiantini et al., 2001). After an air parcel has reached saturation and considering the simplification of an adiabatic atmosphere, the isotopic variation becomes a function of the adiabatic lapse and h. If a portion of the condensate remains in the cloud or moisture recycling occurs within the precipitable system, the δ-elevation relationship results in a weaker trend as a result of the heat exchange required to complete the phase transition. The latter is evident in the broad δ-elevation spectrum across complex mountainous settings in the tropics, particularly in high relief ranges (>3,500 m asl). Mountain shadow effects have also been reported to cancel isotopic variation in leeward slopes (Gonfiantini et al., 2001). On a seasonal basis, the δ-elevation relationship is apparently stronger during wet season rainfalls than dry season events, given precipitation type may exert an additional control over the isotopic variability. When analyzed from individual transects, tropical isotopic lapse rates are weaker (−1.87 ± 0.65‰ km −1 ; see Supplementary Table S1) than the global mean of −2.80‰ km −1 (Poage and Chamberlain, 2001) and the pantropical mean −2.20‰ km −1 (N 348 stations; Adj. r 2 0.48; p < 0.001). Stronger δ-elevation relationships are reported above 3,500 m asl (see Supplementary Figure S1A). Similarly, d-excess increases with elevation (h > 80%), ranging from near +10‰ (0-1,000 m asl) to >15‰ in high relief windward ranges, whereas d-excess < 10‰ is often reported in leeward regions due to strong sub-cloud evaporation (Bershaw, 2018) (see Supplementary Figure S1B). This effect is amplified in areas with rugged topography. Figure 2 shows the distribution of the monitoring stations of rainfall isotopes in Costa Rica, a good example of a complex ocean-land setting. Passive collectors were located across an elevation range from the coastal line up to ∼3,400 m asl, within a variety of maritime (Cocos Island; eastern tropical Pacific ocean) and continental ecosystems across the Caribbean and Pacific domains [see Esquivel-Hernández et al. (2017) for a detailed description of tropical biomes in Costa Rica]. Sampling frequencies were constrained mainly to daily/weekly intervals and event-based occurrences during a major tropical cyclone landfall (Hurricane Otto, 2016;Sánchez-Murillo et al., 2019b). The isotopic lapse rate (windward) is commonly stronger on the Caribbean slope, −1.90‰ in δ 18 O per km of orographic relief (Adj. r 2 0.77; p < 0.001), whereas the Pacific slope (leeward) accounts for only −1.10‰ in δ 18 O per km (Adj. r 2 0.59; p < 0.001). The stronger rain-out effect observed over the Caribbean domain is related to the direct influence of the trade winds (NE prevailing wind direction) (windward) and nearby moisture transport from the semi-closed Caribbean Sea basin, whereas in the Pacific slope the combination of i) the rain shadow effect (leeward), ii) more complex topography and recycled evapotranspiration (ET) moisture, and iii) seasonal deep convective activity (throughout August-October from Pacificoriginated storms) exhibits a weaker and more variable altitude trend. In Central America, Pacific lapse rates are commonly lower than −1.90‰ km −1 , while in the Caribbean and Gulf of Mexico domains windward lapse rates may reach up to −2.40‰ km −1 (Figure 1; Supplementary Table S1). Similar factors affecting leeward/windward isotopic lapse rates have been recognized in other tropical regions (Gonfiantini et al., 2001). The large sensibility of dδ/dz gradients to local circulation features, coupled with the strong seasonal and interannual bias, emphasizes the need for robust and long-term stations (>4 years; Putman et al., 2019), in particular when δ-elevation rates are used to derive critical water management information, such as mean recharge elevations (MREs), either from simple regression analysis or Bayesian methods (Yamanaka and Yamada, 2017;Sánchez-Murillo et al., 2020).

Moisture Transport Mechanisms
In Central America, moisture transport often acts as an "atmospheric bridge" that connects the semi-closed Caribbean Sea basin and the eastern Pacific warm pool. Inter-basin moisture transport is largely modulated by the ITCZ and sea surface temperature (SST) variations. Such interactions exert a large control on rainfall seasonality and distribution (Pan et al., 2018) with variability influenced led by ENSO at interannual scales. In this context, moisture transport was evaluated using the HYSPLIT Lagrangian model (Stein et al., 2015;Rolph et al., 2017). This model uses a three-dimensional Lagrangian air mass vertical velocity algorithm to determine the average position of the air mass, which is reported at an hourly time resolution over the trajectory (Soderberg et al., 2013). Air parcel trajectories were modeled 48 h backwards in time based on the proximity of the Caribbean Sea and the Pacific Ocean. In total, 593 air mass back Frontiers in Earth Science | www.frontiersin.org November 2020 | Volume 8 | Article 571477 trajectories were created and further divided into two main groups: dry (January-April) and wet (May-December) seasons, as a representative example of main moisture transport mechanisms across the Central American Isthmus and surrounding Mesoamerican regions ( Figures 3A,B).
The movement of air masses in the region is primarily dominated by the trade winds so that moisture transport from the Caribbean Sea is favored by the low level wind flow and defined by the Caribbean Low Level Jet (CLLJ, Durán-Quesada et al., 2010). The ITCZ also plays a role in the regional moisture transport as its position either disrupts or enables the moisture supply from the oceanic moisture sources (Durán-Quesada et al., 2017). During the dry season, the ITCZ is at its southernmost position and the low-level flow is mostly zonal from the Caribbean toward southern Central America so that air masses travel directly to the country ( Figure 3A). As the air masses enter the lowlands in the Caribbean domain, the moist air is lifted by the central mountain range; consequently, water vapor experiences a strong adiabatic expansion. This moisture transport pattern is associated with enriched rainfall events, Frontiers in Earth Science | www.frontiersin.org November 2020 | Volume 8 | Article 571477 5 whereby recycled moisture mixes with the air masses traveling along the Caribbean plains. In contrast, during the wet season, the rainfall regime is controlled by the activity of the ITCZ. As the ITCZ moves northward, the zonal component of the CLLJ is weakened and the Choco jet develops (i.e., a low-level westerly jet; Poveda and Mesa, 2000), cross-equatorial winds from the southern hemisphere recurve to transport Pacific parental moisture to Costa Rica ( Figure 3A) as shown by Durán-Quesada et al. (2017) and Morales et al. (2020). As a result, rainfall along the Pacific coast of Central America increases aided by the moisture transport enhancement, the intensification of convective development over the eastern tropical Pacific (Rapp et al., 2014), and the propagation of easterly waves that interact with the ITCZ (Raymond, 2017). Overall, the main moisture source is the southern Caribbean Sea (>90%; Figure 3B) (i.e., Lesser Antilles region). Along moisture transport paths over the isthmus, transpiration (T r ) and re-evaporation (RE) fluxes may reduce the land effect. In an attempt to model this process, a modified water balance approach can be expressed as follows (Leibundgut et al., 2011): where R and R o are isotopes ratios; f is the residual water fraction; α is a cloud-base temperature dependent fractionation factor; and Q a , P, and E correspond to moisture flux from the ocean, rainfall, and evaporation as a function of the trajectory length (τ). However, this model may not apply as the complexity increases (in space and time scales) when adding i) mixing of different vapor sources, ii) input of recycled moisture from T r (across dense vegetation patches), iii) relative contributions between stratiform and deep convective activity, iv) below cloud base evaporation, v) major tropical storm disturbances, and iv) local rural and urban ET fluxes. T r flux is of particular interest, since it accounts between 60-90% of global terrestrial ET (Jasechko et al., 2013;Wei et al., 2017), with the Leaf Area Index (LAI) and the growing vegetation stage collectively explaining 43% of the T r /ET variance (Wang et al., 2014) so that for tropical regions with dense vegetation, like most of Costa Rica, the T r /ET ratio is highly relevant and may require in situ water vapor measurements above and under the canopy cover. In Costa Rica, a recent isotope ecohydrological study in a pristine tropical rainforest, using the spatially distributed tracer-aided rainfallrunoff model for the Tropics (STARRtropics; Dehaspe et al., 2018), resulted in a water partioning of 65% of actual ET being driven by vegetation with greater transpiration rates over the dry season (Correa et al., 2019). Overall, this system was dominated by young water, ranging from 1-hr-transpiration to 3.3 years-old groundwater (Correa et al., 2019).

Temporal and Spatial Isotopic Variability
Water source domains were divided in three groups: Caribbean, Pacific, and Maritime (using Cocos Island data as a valid maritime example; see Corrales-Salazar et al., 2016). Out of 27 monitoring stations (2013-2019), one corresponds to a maritime setting in the eastern tropical Pacific ocean, nine sites were located within the Caribbean slope, four sites received significant inputs from both rainfall regimes (i.e., located near topographic depressions), and 13 stations were installed across the Pacific slope. This spatial distribution ( Figure 2) ensures that representative rainfall modes were captured in the isotopic compositions presented as part of this work. Figure 4 shows a dual δ 18 O-δ 2 H diagram of the main domains. Note that these domains are a common denominator across the Mesoamerican region and the northern portion of South America, and thus our information is pertinent and relevant for a large tropical region, sharing similar rainfall characteristics and facing the same challenges. Similarly, other regions in tropical Africa and southeast Asia are exposed to analogous topographic and ITCZ/ENSO dynamics. Caribbean domain rainfall exhibits a right-skewed density distribution toward more enriched isotope compositions due to the proximity of the Caribbean Sea, less intense Rayleigh-type distillation (<100 km from coast to the main mountain range), and strong convective activity over energy-limited forested watersheds (see Esquivel-Hernández et al., 2017 for a detailed separation of energy and water limited watersheds) ( Figure 4). Isotope ratios resulted in a highly significant meteoric water line (LMWL): δ 2 H 8.05·δ 18 O + 11.48 (Adj. r 2 0.98; p < 0.001). In contrast, Pacific domain rainfall results in a normal distribution with more depleted isotope ratios related to a) more intense Rayleigh-type distillation over the continental divide and b) the influence of a second moisture source. The LMWL for the Pacific slope: δ 2 H 7.92·δ 18 O + 9.75 (Adj. r 2 0.98; p < 0.001) ( Figure 4) suggests the signal to be influenced by sub-cloud evaporation in leeward domains as expected from the convective activity that characterizes rainfall over this area. Cocos Island's LMWL and related density distributions denote the obvious absence of a continental effect but reflect intense moisture re-evaporation from the surrounding eastern Pacific Ocean as expected (LMWL: δ 2 H 8.36·δ 18 O + 13.16, Adj. r 2 0.98, p < 0.001, Figure 4). The strong isotope separation between windward and leeward slopes in tropical landscapes offers unique and traceable isotopic pulses into soil water transport and plant water uptake, surface water networks, and shallow/deep aquifers, which in turn facilitate endmember hydrological modeling. Figure 5 shows a disaggregation of all monitoring sites in terms of δ 18 O and d-excess. In general, enriched right-skewed distributions correspond to the Caribbean domain, whereby d-excess values are commonly between 10 and 15‰. Pacific sites exhibit normal or bimodal distributions. The latter is exhibited at sites where the interaction between two rainfall regimes influences isotopic compositions on a seasonal basis. For example, sites located at topographic depressions and near the continental divide receive highly enriched rainfall from cold fronts between December-February and highly depleted storms during May-October, resulting in clear bimodal patterns. Influence of cold front rainfalls is negligible at Pacific low elevation sites. However, depleted rainfall (as low −22‰ in δ 18 O and up to 20‰ in d-excess) is presented at high elevation sites within the Pacific slope. These results highlight the high spatial variability even for a region considered small, stating the relevance of studying areas featured by a complex Frontiers in Earth Science | www.frontiersin.org November 2020 | Volume 8 | Article 571477  Frontiers in Earth Science | www.frontiersin.org November 2020 | Volume 8 | Article 571477 7 surface as is the case for Costa Rica. Even when the region is surrounded by large water bodies, the results make clear that the region cannot be categorized solely by maritime influences since the continental effect play a major role in driving the observed differences.
Climate seasonality from the dry (January-April) to the wet season (May-December) is characterized by two notable V-or W-shaped isotopic patterns in rainfall ( Figure 6A). These patterns are consistent with the intra-seasonal variability of Central America rainfall, that typically result in two or three depleted excursions during the wet season and two enriched pulses during the mid-summer drought (hereafter MSD; Magaña et al., 1999;Small et al., 2007). Spectral analysis of the δ 18 O time series exhibits a dominant rainfall cycle ranging from 60 to 15 days, with maximum peak activities between 41 and 22 days (Sánchez- Murillo et al., 2019a). Similarly, the local Maiden-Julian Oscillation (MJO) plays a predominant role in modulating rainfall during the MSD in Central America and Mexico (Barlow and Salstein, 2006;Zhao et al., 2019). The latter coincides with the ITCZ northward migration ∼8°N and the intensification of deep convection over the Caribbean domain between July and August, as a result of the activation of tropical easterly waves and the months of the strongest trade winds (January-February). For instance, by mid-May when the ITCZ passes over the region, the first sharp depletion in the isotope composition is commonly observed ( Figure 6A). Slight enrichment occurs during the transition to the MSD. Later in August and October, the isotopic composition shows a second strong depletion, often accompanied by the indirect influence of Caribbean/Atlantic tropical cyclones. As the ITCZ migrates southward by mid-November, a second enrichment is observed toward mid-December. These patterns are amplified within the Pacific slope, whereas in the Caribbean domain (i.e., juxtaposed coasts and inner Caribbean Sea islands) the isotopic composition is less variable throughout the year. The d-excess indicates moisture recycling during the dry season and the influence of cold fronts (December to April), both represented by a notable d-excess increase (up to 20‰) ( Figure 6B). After the onset of the rainy season in mid-May, a prominent decrease in d-excess is observed, primarily in the Pacific slope (leeward), whereby sub-cloud evaporation is the main driver. The latter reaches a minimum during the MSD period, and the second half Frontiers in Earth Science | www.frontiersin.org November 2020 | Volume 8 | Article 571477 of the rainy season features an increase influence of Pacific Ocean moisture transport as the Choco jet activates and the CLLJ is minimum (Figure 3), exhibited d-excess values near +10‰.

Stratiform and Convective Rainfall Fractions
To examine the influence of rainfall types on daily isotopic variability, stratiform rainfall area fractions (F st ) near the Caribbean and Pacific coasts of Costa Rica were analyzed. Point-by-point details on stratiform analysis are presented in Sánchez- Murillo et al. (2019a) and Munksgaard et al. (2020). Positive (enriched) and negative (depleted) inputs in the isotopic composition of tropical rainfall have been linked in paleoclimatic studies to the occurrence of dry and wet conditions, respectively (Frappier et al., 2007;Medina-Elizalde et al., 2016). These types of inferences are often not physically based, and interpretations are even more problematic as modern sampling frequency increases (i.e., monthly to daily). To highlight the last point, Supplementary Figure S2 shows the lack of correlation between daily rainfall, δ 18 O, and d-excess in the longest time series of Costa Rica (N 703; 2013-2019). Moreover, the absence of isotopic composition of water vapor to join rainfall samples allow further complications as the story is partially told, which often drives criticisms of the interpretation. As we are fully aware of the shortcomings and limitations of the stratiform/convective rainfall fraction analysis, we present what is known to propose the need for further research to link the rainfall isotopic composition to cloud microphysics properties. Daily rainfall δ 18 O observations for Costa Rica (N 703; 2013-2019), suggest that the "amount effect" can only explain, on average, 12% of the total variance. However, when aggregated to monthly means, the amount effect explains up to 88% of the variance observed in monthly isotopic composition in Costa Rica (Sánchez- Murillo et al., 2019a). In temperate and maritimeinfluenced environments such as the British Isles (Tyler et al., 2016), the overall correlation between daily δ 18 O and precipitation explained obly 17% of the variance, ranging from 5 up to 45% on spatial basis. Similarly, other continental temperate sites have reported weak precipitation amount correlations on daily basis (Bedaso and Wu, 2020;Adhikari et al., 2020;Bershaw et al., 2020). Overall, the amount effect is not a fully valid approximation over the tropics and is commonly less pronounced on a daily basis compared to monthly time scales as indicated by Araguás-Araguás et al. (2000) (Supplementary Figure S2B). This empirical negative correlation has been extensively used as a valid rationale for paleo-hydroclimate reconstruction in maritime and terrestrial isotope records (i.e., caves, corals, lake sediments, and even tree rings) (Nott, 2004), whereby longer water residence times reflect decadal or centennial precipitation inputs. On the other hand, in some coastal tropical regions, where rainfall modes are significantly affected by inter-annual variability of tropical storms, off-shore convection can result in large and enriched rainfall amounts during tropical cyclones passages and may potentially bias such reconstructions in favor of heavy water isotopes, which researchers may erroneously link to drier conditions via the amount effect.
This common linear regression analysis tends to mask other processes, such as moisture convergence and entrainment, resulting in weak correlations across the tropics with >80% of the variance unexplained when using daily data, whereas stronger correlations (30-70% variance explained) are reported when computing monthly means (Guy et al., 2019;Konecky et al., 2019). These complexities indicate that the isotopic variations in the "amount-effect" dominated region, including the tropics, are not directly controlled by rainfall amount, but rather are influenced by other processes that dominate the phase changes as convection occurs (e.g., cloud microphysics, in-cloud rainfall type differences, moisture transport, influence of up and downdrafts). In the tropics, rainfall can be considered a combination of convective (deep and shallow) and stratiform rainfall fractions (Houze, 1997). Strong convective cells (i.e., featured by strong updrafts >10 km) with potentially shorter moisture recycling times result in more enriched isotope compositions, whereas larger stratiform area fractions are represented by more depleted values (Bolot et al., 2013;Kurita, 2013). Figure 7 shows the relationships between the stratiform rainfall fraction (F st ) and δ 18 O near the Caribbean and Pacific coasts of Costa Rica. This mechanism is stronger in the Caribbean domain than the Pacific slope of Costa Rica (Figures 7A,B). Recent findings across the globe provide additional evidence of such a mechanism (Aggarwal et al., 2016;Zwart et al., 2018;Sánchez-Murillo et al., 2019a;Konecky et al., 2019;Munksgaard et al., 2019). Strong isotopic correlations with F st have recently been reported in Costa Rica, Ghana, Singapore, Malaysia, Bangladesh, Australia, and Vietnam (Munksgaard et al., 2019). Such findings highlight the importance of event-based sampling so that F st can be better interpreted in relation to the life cycle of convective activity, which is relatively short lived (1-6 h and 20-200 km; Orlanski, 1975).
Our isotopic observations demonstrate that enriched isotopic ratios occur during both large and small TC-rainfall events within the Central America region. Likewise, depleted isotopic ratios can also occur during both large and small rainfall events. However, during the study interval, tropical storms were the most common origin of significantly depleted isotopic ratios, whereas the fartherfield (i.e., long distance rainfall bands from the cyclone center) indirect effect of hurricanes was represented by more enriched compositions ( Figure 8A). As the structure of the cyclones becomes more organized, the variability of δ 18 O decreases to favor enriched values, reflecting the latent heat release that characterizes the organization of convection in a hurricane.
Additional evidence is provided in Welsh and Sánchez-Murillo (2020) and presented in Figure 8B, which depicts a dual isotope diagram that illustrates the isotopic evolution of Hurricanes Otto (2016) and The Bahamas (Nassau) (Lat: 25.0739, Long: 77.4138, elevation: 11 m asl) with event-based, 6-hourly and daily time intervals before, during, and after TCs, including both continental landfalls and maritime passages. These samples resulted in the first study of high-frequency isotope composition during extreme TCs within the basins of the Caribbean Sea and Atlantic Ocean (Welsh and Sánchez-Murillo, 2020). These cyclones were some of the strongest storms in recent history to affect this region, causing destruction in Abaco and Grand Bahama islands from Hurricane Dorian, and the Costa Rican northern lowlands from Hurricane Otto. This dataset represents the only isotopic fingerprint of these TC events. Maritime passages either in the Greater Antilles or near the Caribbean coast of Central America are represented mainly by enriched compositions, whereas continental landfalls exhibited more depleted ratios.
The interpretation of enriched values for TCs becomes relevant considering that most of the significant paleo-archives of the Mesoamerican region are located along the Caribbean coast. Across this region, off-shore convection can result in large and enriched rainfall amounts during TC passages and may potentially bias such reconstructions in favor of heavy water isotopes, typically linked to drier conditions. In this regard, it is important to point that enrichment can be associated to both extreme drought and the influence of a major hurricane, suggesting the need to re-evaluate the criteria used for the analysis of isotopic proxies in paleoclimate archives as extreme events could be marked by similar values.

El Niño/Southern Oscillation Modes and Isotope Variability
In Central Costa Rica, our high-resolution rainfall isotope data are unique in the sense that they comprise ENSO's neutral, warm (very strong and weak), and cold (moderate) phases. In general, warm ENSO events result in a large rainfall deficit across the Pacific slope of Central America and a surplus over the Caribbean slope  Figure 8C).
Moreover, results exposed make clear that isotopic monitoring reveals the evolution of warm ENSO events with a better accuracy than for cold ENSO events. The relevance of this result is that continuous monitoring reveal information of the severity of the impact of ENSO at the local scale. Hence, systematic monitoring can leverage diagnosis of the local impacts of ENSO and potentially support seasonal forecasts based on SST to better inform the most affected sectors. In the case of Costa Rica, the latter is highly relevant as the evolution of warm ENSO is linked with severe drought development that causes negative impacts in the agricultural sectors and poses a threat to water resources availability across the country (Muñoz-Jiménez et al., 2019).

GROUNDWATER AND SURFACE WATER CONNECTIVITY
The distinct isotopic pulses or input functions provided by rainfall in the tropics are commonly damped or modified by different ecological and hydrogeological processes across the landscape. Connectivity between groundwater and surface water compartments have been widely studied in humid tropical watersheds (Gordon et al., 2015;Lagomasino et al., 2015;Batista et al., 2018;Oiro et al., 2018;Villegas et al., 2018;Welsh et al., 2018;Tsuchihara et al., 2019); however, large scale or regional approaches are uncommon across tropical biomes. Figure 9 shows locations of groundwater (N 1,156) and surface water (N 603) sampling efforts across Costa Rica (2016-2019). Samples were collected targeting baseflow periods (December to April) to better represent mean annual recharge and discharge isotopic values. In general, groundwater and surface water exhibited a strong meteoric origin, with secondary evaporation concentrated in Pacific rivers ( Figure 10A).
To better understand groundwater recharge and runoff generation processes, and through an exploratory approach, discrete soil samples from two different rainforests (central and northern Costa Rica) have recently been evaluated for water stable isotopes during two contrasting seasons. Two different methods were used to extract soil water: centrifugation (extracts the most available moisture in soils) and cryogenic distillation (extracts the complete moisture hold in soils) (see Orlowski et al., 2016 for a detailed description on methods). Soil water from a rainforest at central Costa Rica was extracted by centrifugation (topsoil sampling during wet season; red triangles Figure 10A; Sánchez- Murillo et al., 2019b), while soil water coming from a rainforest at northern Costa Rica was extracted through cryogenic distillation (sampled between 5 and 110 cm during dry season; yellow dots Figure 10A; Dehaspe et al., 2018). Indistinctively of soil depth, extraction method, season or region, water extracted from soil samples exhibited a clear meteoric origin in line with rainfall isotope seasonality and with minimal enrichment. Remarkably, the coincidence of ground and soil water along the CR MWL ( Figure 10A), suggests no strong evaporative effect and that there is a single soil water pool, rather than two soil water pools as suggested by the ecohydrological separation hypothesis (mobile vs. immobile water) reported by other studies worldwide (e.g., Evaristo et al., 2015;Luo et al., 2019). These results are coincident with those from a recent article by Jiménez-Rodríguez et al. (2019), who showed a lack of evaporation and coincidence with the LMWL of soil samples of a wet lowland forest of Costa Rica during the dry season. Nevertheless, the preliminary soil water isotope data presented here for rainforests heavily influenced by Caribbeantype rainfall regimes, invokes more in depth studies with higher sampling frequency, and across other tropical ecosystems whereby d-excess or lc-excess values may indicate stronger evaporation controls in the vadose zone. More efforts on studying how soil water patterns are traduced into plant water uptake and compartmentalization are also needed. In this sense, Jiménez-Rodríguez et al. (2019) already showed a compartmentalization of plant water use during the dry period between different plant types for the rainforest understudied. Further studies following a tracer-aided catchment perspective for tracking water fluxes with a high sampling frequency and across different tropical ecosystems are needed.
In the Caribbean slope, the strong connectivity between groundwater (−4.8‰) and surface water (−4.9‰) compartments is reflected as nearly equal median values ( Figure 10B). In contrast, in the Pacific slope strong secondary evaporation processes are reflected in more enriched δ 18 O median values in surface water (−7.2‰) when compared to groundwater (−7.9‰) ( Figure 10B). The latter is well represented by regional d-excess, in which the Caribbean slope exhibited values commonly >10‰ while values as low as −9‰ were reported in the Pacific domain ( Figure 10C). To emphasize regional differences between the Caribbean and Pacific domains, rainfall to groundwater spatial isotopic ratios (P/ GW) were constructed following Jasechko and Taylor (2015) to estimate the isotopic recharge bias across Costa Rica. Mean δ 18 O annual values in rainfall and groundwater were extracted from high-resolution (100 × 100 m grid) isoscapes   (Figure 11). Interpolations were based on a thin plate spline algorithm (Franke, 1982). The rainfall isoscape was calculated using amount-weighted isotope ratios, while the groundwater and surface water isoscapes were constructed using mean annual values. P/GW were computed using raster calculations based on rainfall and groundwater isoscapes. A P/GW > 1 indicates regions where infiltrated water is susceptible to near surface secondary evaporation, while a P/GW < 1 points toward recharge originating from more intense and more depleted rainfalls at high altitudes. A P/GW ∼ 1 indicates a rapid recharge via preferential flow paths. Caribbean and Pacific isotopic lapse rates were used to calculate potential mean recharge elevations (MRE in m asl) under the assumption that groundwater isotope ratios under baseflow conditions are representative of mean annual recharge conditions (Yamanaka and Yamada, 2017;Sánchez-Murillo et al., 2020). Isoscapes show a remarkable spatial variability of δ 18 O, with an enriched composition over coastal areas as a result of rainfall forming at lower altitudes and a low-level convergence favored by higher temperatures. Groundwater and surface water isoscapes reflect the strong orographic separation between the Caribbean and Pacific domains, as well as more enriched values in windward regions vs. depleted values in leeward slope (Figure 11). The interpolated P/GW isotope ratios ( Figure 12A) show areas where P/GW > 1, which indicates a significant soil matrix-water attenuation in the unsaturated zone and potentially implies shallow groundwater reservoirs, particularly across coastal aquifers within the Pacific domain. Areas relying on high altitude recharge with depleted waters are represented by P/GW < 1. These areas are mainly located across the highaltitude interior, such as the central mountainous range, which provides groundwater resources to over half of Costa Rica's population. The spatial pattern of P/GW isotope ratios that increase from the mountain range down to the coast are most prominent in the Central Pacific region and support the notion of groundwater that is mainly recharged at higher altitudes. P/GW isotope ratios close to 1 are an indication of rapid groundwater recharge often observed across the Caribbean lowlands. However, in high-elevation water systems located in the Talamanca range (∼3,400-3,600 m asl), and where the contribution of groundwater to the surface systems can be considered minimal, evaporation conditions are highly influenced by input from precipitation that is in isotopic equilibrium with local water vapor, yielding low evaporation to inflow ratios (<20%) and relatively short residence times of water in the glacial lakes (Esquivel-Hernández et al., 2018). FIGURE 10 | Dual isotope diagram and scattered box plots of rainfall, groundwater, and surface waters across the Caribbean and Pacific domains. (A) δ18O and δ2H relationship between groundwater (cyan squares), surface water (pink crosses) and soil water (yellow dots and red trinagles) compared to the Costa Rica meteoric water line (Sánchez- Murillo et al., 2013). LEL stands for local evaporation line. B) δ 18 O (‰) and C) d-excess (‰) box plots include 25th, 75th, median, and outliers for each group. Figure adapted using data from Sánchez- .
Frontiers in Earth Science | www.frontiersin.org November 2020 | Volume 8 | Article 571477 Mean recharge elevations (MREs) obtained from isotopic lapse rates (Yamanaka and Yamada, 2017) reinforce the P/GW previously inferred mechanisms, in which: a) high altitude recharge sites are located across the main mountain range (1,976-3,218 m asl), b) main urban areas exhibited MREs between 1,331-1,975 m asl, and c) coastal and distant lowland regions resulted in lower elevations (<1,330 m asl) (Figure 12). A large disproportionality exists between MREs and water use elevations, here denoted by the sampling site elevation (SE) ( Figure 12B). Median MRE exhibited a large spectrum with a median value > 1,000 m asl; however, MRE-SE differences can reach up to 2,500 m asl, indicating in some cases regional groundwater flowpaths.
In tropical regions elevation can impact the provision of water for domestic supply, which is complex and usually involves many sources (i.e., deep and shallow groundwater, surface water, and spring water). Overall, tap water at lowland urban settlements disproportionally depends on high elevation recharge from FIGURE 11 | Rainfall and groundwater annual spatial variations. Interpolated δ18O (A) rainfall, (B) groundwater, and (C) surface water isoscapes (Sánchez- . Figure adapted using data from Sánchez- . "mountain water towers" (Immerzeel et al., 2020). Isotopic metrics serve to better evaluate conservation and regulation needs in critical mountainous or plain recharge areas. However, in tropical regions significant temporal and spatial gaps require more systematic efforts that account for better: i) sampling gradients to improve isotopic lapse rates; ii) surface water and groundwater sampling during baseflow regimes to characterize mean annual conditions; iii) isoscapes in order to trigger a positive feedback between prediction and observation efforts; iv) forest conservation and land use practices; and v) protection laws in critical recharge elevations.

MAGMATIC FLUIDS IN VOLCANIC AQUIFERS
The geological/tectonic characteristics of Costa Rica offer the necessary heat and decompression cooling for storing fluids (i.e., meteoric water, sea water, magmatic water or mixed fluids including gas vents; Molina and Martí, 2016) within fractured volcanic aquifers across both Pacific and Caribbean slopes. Water stable isotopes studies have been used to investigate hydrothermal manifestations in Central America for almost four decades (Goff et al., 1987;Giggenbach and Soto, 1992;Rowe et al., 1995;Nieva et al., 1997;Tassi et al., 2005;López et al., 2006;Birkle and Bundschuh, 2007;Rouwet et al., 2009;Molina and Martí, 2016); however, these efforts have been primarily focused on the understating of volcanic and seismic activity.
Here, we present a Bayesian mixing isotope model (R package Simmr; Parnell and Inger, 2016) to determine the influence of andesitic water (i.e., waters associated with more silicic andesitic magmatism typical of mature volcanic arc systems; Taran et al., 1989;Giggenbach, 1992) in hydrothermal manifestations across Costa Rica. Hydrothermal samples (N 409) correspond to stable isotope archives and recent sampling efforts in Costa Rica (see Ramírez-Leiva et al., 2017 for sampling details). Sampling sites comprised mud-pools, hot springs, mountainous springs or seepages, and drilled and dug wells. End-members were classified as i) Pacific slope rainfall (δ 18 O -7.47 ± 3.84; δ 2 H -49.90 ± 30.75), ii) Caribbean slope rainfall (δ 18 O -4.26 ± 3.69; δ 2 H -22.82 ± 30.06), and iii) andesitic water (δ 18 O +10 ± 3; δ 2 H +20 ± 10) (Giggenbach, 1992). Simmr was implemented with 100,000 iterations discarding the first 10,000. No prior information was used resulting in an equal likelihood of all sources (Zhang et al., 2017). The estimated source water contributions to the volcanic aquifer were determined using a Markov Chain Monte Carlo function (Brooks, 1998) to repeatedly estimate the proportions of the various sources in the mixture and determine the values that best fit the mixture data (Parnell and Inger, 2016). The median (50% quantile) source water contribution from the posterior parameter distribution was used for practical comparisons.
In Costa Rica, δ 18 O in hydrothermal systems ranged from -11.5‰ to +6.2‰ with a mean of -4.5 ± 2.4‰, whereas δ 2 H composition varied from -82.3‰ to +7.6‰ with a mean of -34.2 ± 9.0‰ (Figure 13). The hydrothermal water line (HWL) (i.e., evaporation line) of Costa Rica can be described as: δ 2 H 4.70 · δ 18 O − 13.0 (r 2 0.79, N 409, p < 0.001). Similarly, d-excess varied from −63.1‰ to +28.4‰ with a mean of +1.8 ± 9.8‰ and lc-excess values ranged from −68.2‰ to +18.4‰ with a mean of −7.3 ± 9.0‰. In the tropics, hydrothermal manifestations often receive windward and leeward rainfall with distinct isotopic ratios. In the case of Costa Rica, the large spectrum of isotopic values in hydrothermal fluids are explained by the strong isotopic separation in rainfall between the Pacific and Caribbean slopes and further emphasize the importance of understanding recharge contributions from two distinct meteoric water sources under the influence of andesitic water contribution within the aquifer.
Bayesian mixing proportions denote a dominance of meteoric water over andesitic inputs. Meteoric waters are responsible for approximately 78% of the total contribution to volcanic aquifers and can be divided by origin as follows: Pacific rainfall (42.3 ± 2.4%) and Caribbean rainfall (36.0 ± 3.2%), whereas andesitic water contributes 21.7 ± 1.6%. Most of the hydrothermal manifestations are in the front volcanic arch (Pacific slope). Previous studies reported mean andesitic water contributions between 16.4% and 19.7%, based on simple two end-member models (Giggenbach and Soto, 1992;Molina and Martí, 2016;Ramírez-Leiva et al., 2017). The source of the andesitic water is most likely recycled seawater, which enters the subduction Central America zone as water-bound in hydrous mineral sub-ducting slab (Ryan and Chauvel, 2013;Zamboni et al., 2016). The large spectrum of andesitic water contribution (10-30%) ( Figure 13B) within the hydrothermal systems of Costa Rica largely reflects the degree of interaction and exchange with isotopically "pre-shifted" recycled fluids (i.e., fluids generally enriched in δ 18 O than meteoric waters). Further validation of this stable isotope approach should involve weekly to monthly sampling in hydrothermal fields, including surface water discharge and deep/ shallow wells to better constrain the temporal evolution of the isotopic ratios (both in rainfall and hydrothermal inputs) during the rainy and dry seasons, but most importantly, to separate the effect of surface kinetic fractionation vs. underground evaporation and mixing with "pre-shifted" recycled subduction fluids. This information may be used to better understand the hydrothermal/ volcanic activity, as well as providing new spatially-oriented insights of recycled fluids across tropical subduction zones.

TRACER-AIDED RAINFALL-RUNOFF MODELING
Here, we present a short overview of models that used water stable isotope data to estimate flow pathways, runoff generation processes and water ages in Costa Rica and beyond as examples from the humid tropics. The simplest approaches to assess flow pathways and runoff generation in the tropics used mixing models with work as early as 1978 by Bonell and Gilmour (1978) from the famous Babinda experimental catchment in northern Queensland, Australia. Elsenbeer et al. (1995) followed on from Bonell and Gilmour (1978) and expanded their work to the Amazon with similar methods (Elsenbeer and Lack, 1996). In the humid tropics of Costa Rica at the La Selva biological station, Genereux and Pringle (1997) were the first to assess flow pathways using in their case Frontiers in Earth Science | www.frontiersin.org geochemical tracers and mixing models. The latter work was expanded using a suite of environmental tracers for age dating of waters at La Selva by Solomon et al. (2010). Among the first attempts to estimate transit times (TT) of a tropical catchment in southern Costa Rica was to the best of our knowledge , who used a simple lumped convolution integral model (Małoszewski and Zuber, 1982) with a Gamma distribution as transfer function (Kirchner et al., 2000). Elsewhere, Goller et al. (2005) used stables isotope modeling to infer flow pathways in a rainforest catchment in Ecuador. Similarly,  found a relatively simple exponential distribution to be an effective model to simulate the baseflow TT of a high-elevation Andean Páramo catchment in Ecuador. The short-term event TTs of a tropical catchment in Colombia was assessed by Roa- García and Weiler (2010). Soil water TTs of a montane tropical cloud forest in Mexico were studied by Munõz-Villers and McDonnell (2012). Also, in Mexico but from a much more seasonally-dry tropical site, Farrick and Branfireun (2015) used TT models to assess flow pathways. Duvert et al. (2016) used stable isotope TT models in combination with tritium and chloride to assess groundwater contributions to tropical streams in Australia. Most of the TT modeling is based on the assumption of stationarity and few attempts were made to assess non-stationary TT in the tropics.  applied the gamma model with a moving time window assessing non-stationary TTs driven by the climatic seasonality of a pronounced tropical dry and wet season.
Following this research and to specifically account for timevarying water ages that do not necessarily follow an a priori distribution such as used with the lumped convolution integral approach, Birkel and Soulsby (2016) incorporated the observed stable isotope time series into lumped conceptual tracer-aided models of varying complexity assessing dominant flow pathways and runoff generation processes in the same southern Pacific catchment in Costa Rica. Increasing in complexity but still a conceptual model approach, Dehaspe et al. (2018) modified the spatially-distributed model of van Huijgevoort et al. (2016) to incorporate tropical rainforest interception (the INT upper blue reservoir in Figure 14) for the small, experimental San Lorencito catchment in northwestern Costa Rica. The same model was further adapted to incorporate eco-hydrological water partitioning accounting for transpiration fluxes (the green TR reservoir in Figure 14) based on stable isotope mass balances. The latter model allowed to explicitly quantify transpiration fluxes and additionally calculates water flux ages and the ages of stored water using a flux tracking approach similar to Hrachowitz et al. (2013).
The latter approach tags every incoming rainfall with a time stamp as water moves through the two-reservoir model cascade toward the stream conceptualizing soil and groundwater stores ( Figure 14). The models described here, ranging from simple two-parameter approaches to fully spatially distributed with 12 calibrated parameters, are presented and discussed in more detail in the review paper by Birkel and Soulsby (2015). Recent approaches used stable isotope tracers of carbon in combination with mixing models, TT models and conceptual hydrological models to assess fluvial tropical carbon export. Examples from Australia include among others, Duvert et al. (2020) and Birkel et al. (2020), and from Costa Rica, Genereux et al. (2013) and Sanchez-Murillo et al. (2020). Albeit other tracer-aided hydrological modeling approaches exist, such as particle tracking (Davies et al., 2011) and the recent Storage Age Selection (SAS) analytical models (e.g., Benettin et al., 2017), none to the best of our knowledge was yet applied in the tropics.

Future Challenges
Tracer hydrology studies are substantially increasing in the tropics (Sánchez-Murillo and Durán-Quesada, 2019), although FIGURE 13 | Andesitic water contribution. (A) Dual isotope diagrams for hydrothermal mixtures (gray empty circles) and mean values (error bars) of three endmembers: Pacific rainfall (purple), Caribbean rainfall (cyan), and andesitic water (yellow) as proposed by Giggenbach (1992). (B) Source contribution (as proportion) box plots. Figure  Frontiers in Earth Science | www.frontiersin.org November 2020 | Volume 8 | Article 571477 significant challenges remain and several aspects of ecohydrological processes are still poorly understood, such as: a) physically-based meaning of 17 O-excess variations, b) rainfall generation (i.e., deep and shallow convective development, stratiform fractions, and in-cloud fractionation kinematics) and its correlation with isotope variations in the context of high resolution water vapor data, c) coordination of continuous long-term stable isotope monitoring stations (precipitation, soil, plant, surface, and groundwater), d) standardization of isoscaping modeling and spatial coverage of sampling efforts, e) appropriate incorporation of stable isotopes measurements into GCMs, paleo-climatic, and ecohydrological models, f) high frequency measurements during extreme events such as tropical cyclone passages and landfalls, g) tracer-aided and Bayesian mixing models, h) incorporation of tracer hydrology courses in a wide range of related university majors, and i) foster a more direct connectivity with the ecology community.
Data Gaps: Water Vapor, Fog, Soil, and Plants In terms of data gaps, water vapor and fog isotope data remain poorly constrained across the tropics, as well as studies tracking water flowing through the soil-plant-atmosphere continuum (SPAC) with isotopes. However, despite the inherent sampling complexity of vapor, recent studies go from in-situ sampling (Lawrence et al., 2004), on-line measurements (Thurnherr et al., 2020), aircraft (Iannone et al., 2009), and ship (Benetti et al., 2017a) based to satellite (Worden et al., 2006). In the tropics, such measurements and estimations have been along with other meteorological information used to understand moistening in the mid-troposphere (Noone, 2012) and moisture anomalies (Bailey et al., 2017), atmospheric circulations (Dee et al., 2018) to mention some applications. In tropical montane cloud forests, fog featuring the cloudy nature of these ecosystems, has been shown to be highly sensitive to temperature variations (Brujinzeel et al., 2011;Rapp and Silman (2012)). Furthermore, such forests are often located in the higher elevation of the tropics where significant warming trends have been identified (Hu and Riveros-Iregui, 2016;Los et al., 2019). In this context, understanding the connectivity between temperature changes for the water vapor saturation and their effect in montane cloud forest fog is deemed as highly relevant to better understand the potential impacts of climate change related warming for this unique ecosystem. Isotopic composition of water vapor in these environments of high saturation comes in as a handy tool to understand this connectivity. The combination of the isotopic composition of fog  (Dehaspe et al., 2018) and further modified to include an eco-hydrological water partitioning. The coupled model simulates streamflow (observations in red and simulations in black and blue) and isotopes (here: δ 2 H composition) as catchment integrated outputs (upper right panels). The lower two panels on the right exhibit stream water age estimates in days as the model can be used to track the age of water in fluxes and storages. The STARRtropics simulated cumulative distribution function (CDF) of stream water ages is compared to a best-fit lumped convolution integral model using a Gamma distribution as transfer function. Both models result in comparable mean transit times (MTT) of around 1 year with a relatively similar form of the transit time distribution. Figure adapted using data from Dehaspe et al. (2018).
Frontiers in Earth Science | www.frontiersin.org November 2020 | Volume 8 | Article 571477 along with other endmembers has been proved to segregate fog contribution to the balance of saturation required in montane cloud forests to preserve their moisture profiles. Studies reveal that the isotopic composition of fog is featured by an enriched signal compared to rainfall because of differences in condensation and that the deviation between fog and rainfall isotopes are enlarged in presence of rainfall linked with synoptic systems, in contrast with more locally generated fog (Scholl et al., 2010). Studies on isotope ecohydrology across the tropics are emerging, and so far they have addressed: plant water uptake and partitioning (Meinzer et al., 1999;Evaristo et al., 2016;De Wispelaere et al., 2017;De Deurwaerder et al., 2018;Bodé et al., 2020;Jiménez-Rodríguez et al., 2019), soil water dynamics and mixing (Mosquera et al., 2020), and soil evaporative losses (Hasselquist et al., 2018). Most of these studies come from the dry tropics; however, research coming from mountain and wet rainforest environments, as well as highland Páramos, are particularly scarce. Across the tropics, and in particular in the mentioned environments, fog is a major driver of water inputs and outputs (Aparecido et al., 2018). Given the expected difference in fog's isotopic composition compared to rainfall (Scholl et al., 2010), the combined use of isotopic data with other ecophysiological measurements (e.g., sap flow) can enlighten our understanding of fog's relevance on tropical plant water partitioning and mechanisms underlying, like foliar water uptake and hydraulic redistribution (see Oliveira et al., 2014 for an overview). Moreover, field studies measuring together fog's presence and tracking it through the SPAC with stable isotopes could provide insightful evidence of tropical ecosystems sensitivity to climate change, as shown through glasshouse experiments by Eller et al. (2016). Estimating the residence times of water in plants and soil in the tropics is still poorly constrained. Overall, by combining samplings at different spatial scales and while integrating ecology and hydrology, a dual isotope approach for tracking water through the SPAC is promising to improve our understanding of the ecophysiology and ecohydrology of the changing tropics .
The Sampling Dilemma Figure 15 shows a conceptual diagram of the spatiotemporal sampling dilemma, including a detailed overview between time-step selection and study area (i.e., lakes, groundwater reservoirs, streams, springs, soil-plants, and urban environments) across multiple precipitation scale processes. From an experimental design and tropical perspective, sampling frequencies (i.e., real-time, event-based, subhourly, daily, weekly, or monthly composites) and time series length constrain the hydro-meteorological process to be analyzed. Frequently, research questions are not fully coupled with the appropriated endmember sampling design. Likewise, different sampling frequencies may represent similar isotope variations and may involve distinct sampling efforts (see Supplementary Figure S3 for a daily vs. weekly rainfall sampling in a tropical setting).
FIGURE 15 | Conceptual diagram of the spatiotemporal sampling dilemma. X and Y axes show the potential scales of water cycle components and artificial features and water residence times, respectively. Right panel includes different features of precipitation regimes from decadal climate variability to cloud thermodynamics. Dashed ovals denote a potential time/area boundary for each water compartment or feature. Gray parentheses represent suggested stable isotope sampling frequencies.
Frontiers in Earth Science | www.frontiersin.org November 2020 | Volume 8 | Article 571477 For instance, cloud thermodynamics and dynamical processes may require continuous monitoring at a resolution of at least 1 Hz to capture the primary cloud processes and vapor mixing , while monthly composited samples isotopic composition may only represent large scale regional processes or the signature impinged by an isolated extreme weather event such as tropical cyclones across tropical semiarid basins. A similar situation has recently been observed by  for ecohydrological studies, as new isotopic data sets with higher resolution than monthly sampling or seasonal snapshots are needed to properly cover different soil/ vegetation stages along drying-rewetting cycles around the year.
Weekly to hourly sampling frequencies are suitable in studies comprising urban runoff, wastewater cycles, and aqueduct dynamics and are usually linked to short-term strong convective variations. Weekly to monthly sampling is more appropriate for springs, canals, lakes/reservoirs and groundwater systems, whereby longer residence times (from few months to years) may also be correlated with seasonal to inter-annual precipitation modes. Plants and soil sampling can be more appropriate between hours to weeks or months, depending on the objective of the study, if it aims at detecting diel fluctuations, event response or seasonal patterns. Long-term isotope data sets (>10 years) are needed to capture warm and cold ENSO phases and potential oceanic decadal oscillations; however, these types of data sets are rather limited (Sánchez-Murillo et al., 2019a;Ellis et al., 2020;Oza et al., 2020).

Modeling Advances and Limitations
Tracer-aided and Bayesian mixing models are surging as the most prominent tool for the assimilation of tropical isotope observations to scale up (from small to meso-and large scale watersheds) the understanding of hydrological processes including sediment and solute transport (i.e., trace elements, dissolved organic matter). Recent advances in multi-source (N > 3) partitioning uncertainty estimations using large isotopic data sets (Correa et al., 2019) also offer an alternative method to the simple sum of analytical errors or the Bayesian uncertainty approach. However, future modeling success will crucially depend on more and novel measurements in tropical regions particularly of soil water and plant water isotopes to test our models against (Berry et al., 2018;Kuppel et al., 2018).
Furthermore, the ability of water stable isotopes to provide a multi-scale connectivity have been largely constrained to rainfall, groundwater and streamflow sampling with very little evidence of vapor isotopic composition so that the closure of the water budget has been constrained. Even when combined with more complex techniques to measure surface fluxes, the mass balance at shorter time scales is missing the local effect of diffusion in absence of continuous water vapor measurements. In this sense, field campaigns can provide a good opportunity to assess this problem and enable the inclusion of the vapor phase into more complex models which would greatly benefit the improvement of the present representation of the surfaceatmosphere coupling both over land and ocean.

AUTHOR CONTRIBUTIONS
All authors listed have made a substantial, direct, and intellectual contribution to this review, and approved it for publication.

ACKNOWLEDGMENTS
The authors thank various anonymous helping hands that contributed to rainfall, surface water, hydrothermal, and groundwater sampling in Costa Rica in different environmental challenging conditions. Support from the Isotope Network for Tropical Ecosystem Studies (ISONet) funded by the Research Council of Universidad de Costa Rica is also recognized. Monitoring at various National Parks was conducted in agreement with the National System of Conservation Areas (SINAC) of Costa Rica.