Articles | Volume 23, issue 1
Research article
24 Jan 2019
Research article |  | 24 Jan 2019

Multi-scale temporal variability in meltwater contributions in a tropical glacierized watershed

Leila Saberi, Rachel T. McLaughlin, G.-H. Crystal Ng, Jeff La Frenierre, Andrew D. Wickert, Michel Baraer, Wei Zhi, Li Li, and Bryan G. Mark

Climate models predict amplified warming at high elevations in low latitudes, making tropical glacierized regions some of the most vulnerable hydrological systems in the world. Observations reveal decreasing streamflow due to retreating glaciers in the Andes, which hold 99 % of all tropical glaciers. However, the timescales over which meltwater contributes to streamflow and the pathways it takes – surface and subsurface – remain uncertain, hindering our ability to predict how shrinking glaciers will impact water resources. Two major contributors to this uncertainty are the sparsity of hydrologic measurements in tropical glacierized watersheds and the complication of hydrograph separation where there is year-round glacier melt. We address these challenges using a multi-method approach that employs repeat hydrochemical mixing model analysis, hydroclimatic time series analysis, and integrated watershed modeling. Each of these approaches interrogates distinct timescale relationships among meltwater, groundwater, and stream discharge. Our results challenge the commonly held conceptual model that glaciers buffer discharge variability. Instead, in a subhumid watershed on Volcán Chimborazo, Ecuador, glacier melt drives nearly all the variability in discharge (Pearson correlation coefficient of 0.89 in simulations), with glaciers contributing a broad range of 20 %–60 % or wider of discharge, mostly (86 %) through surface runoff on hourly timescales, but also through infiltration that increases annual groundwater contributions by nearly 20 %. We further found that rainfall may enhance glacier melt contributions to discharge at timescales that complement glacier melt production, possibly explaining why minimum discharge occurred at the study site during warm but dry El Niño conditions, which typically heighten melt in the Andes. Our findings caution against extrapolations from isolated measurements: stream discharge and glacier melt contributions in tropical glacierized systems can change substantially at hourly to interannual timescales, due to climatic variability and surface to subsurface flow processes.

Please read the corrigendum first before continuing.

1 Introduction

Glaciers supply water resources to over 600 million people worldwide (Messerli et al.2004). By melting during dry seasons and drought years, they supplement streamflow (Chen et al.2017; Escher-Vetter et al.1994; Fountain and Tangborn1985; Jansson et al.2003; Juen et al.2007; Lang1986; Soruco et al.2015) and ensure reliable water supplies (Bury et al.2013; Mark and Mckenzie2007; Mark and Seltzer2003). This has led to the commonly held conceptual model, called the “glacier compensation effect” (Lang1986), in which meltwater buffers discharge variability.

Climate change can disrupt the glacier compensation effect, and tropical glacierized watersheds that already experience year-round melt (Kaser and Osmaston2002) may be the most vulnerable. Climate models predict amplified temperature increases at high altitudes in low latitudes (Bradley2006; Pepin et al.2015). The retreat of these glaciers temporarily results in increased runoff (Braun et al.2000; Carey et al.2017; Mark2008; Polk et al.2017), but gradually depletes the storage of these mountain “water towers”. Over time, this reduction in storage capacity can render these glaciers unable to supply sufficient dry-season meltwater discharge for the communities that depend on it (Barnett et al.2005; Bradley2006; Luce2018; Mackay2008; Ostheimer et al.2005). Indeed, observations already reveal reduced and fluctuating flows in glacierized watersheds (Baraer et al.2012, 2015; Huss et al.2008; Mark and Seltzer2003; Rabatel et al.2013; Soruco et al.2015), threatening the water security of millions of people (Carey et al.2017; Immerzeel et al.2010; Vuille et al.2018).

Of all glaciers in the tropics, 99 % are located in the Andes (Kaser1999), often in remote regions, where resource-limited populations often rely on their meltwater (Bury et al.2011; La Frenierre and Mark2017). Despite over a decade of research in Peru's heavily glacierized Cordillera Blanca (Baraer et al.2012, 2015; Juen et al.2007; Mark and Mckenzie2007; Mark and Seltzer2003; Mark et al.2005), many of the processes linking variability in climate, glacier melt, and stream discharge remain uncertain. For example, groundwater is also a major contributor to discharge in many glacierized mountainous watersheds around the world (Andermann et al.2012; Baraer et al.2009, 2015; Clow et al.2003; Engel et al.2016; Harrington et al.2018; Hood et al.2006; Huth et al.2004; Liu et al.2004; Pohl et al.2015; Schmieder et al.2018; Somers et al.2016; Tague and Grant2009; Tague et al.2008). This can further modulate discharge through baseflow, but its capacity to do so as glaciers respond to climate change is complicated by largely unconstrained relationships between glacial meltwater and groundwater recharge (Baraer et al.2015; Favier et al.2008; Gordon et al.2015; Harrington et al.2018; Minaya2016).

Understanding how different surface and subsurface pathways influence the timing of meltwater and groundwater contributions to streamflow is critical for predicting how climate change will impact the reliability of watershed discharge. A major challenge in evaluating these spatiotemporal effects in tropical glacierized watersheds is the relative data sparsity and resource limitations in these regions compared to better instrumented mountainous systems in North America and Europe. Many studies in tropical and other remote glacierized settings rely on focused field campaigns using methods such as synoptic water chemistry tracer sampling (Baraer et al.2009, 2015; Mark and Mckenzie2007; Wilson et al.2016), but these provide only snapshots of the hydrologic state. Even though physically based hydrologic models can provide greater spatiotemporal coverage in mountainous settings (He et al.2018; Liu et al.2004; Lowry et al.2010, 2011; Markovich et al.2016; Omani et al.2017; Pribulick et al.2016; Suecker et al.2000; Tague and Grant2009; Tague et al.2008), their application is relatively limited in Andean watersheds (Buytaert and Beven2011; Minaya2016; Ng et al.2018; Omani et al.2017) due to the lack of extensive monitoring infrastructure. With these obstacles, there remains limited understanding of how stream discharge in tropical glacierized watersheds varies over timescales ranging from hours to years, and how this variability is driven by dynamic inputs of glacial meltwater and precipitation through a combination of surficial and subsurface pathways.

In this study, we probe the multiple timescales of hydrological processes in a subhumid glacierized watershed on Volcán Chimborazo in the tropical Ecuadorian Andes. Prior to this work, there have been no comprehensive efforts on Chimborazo to quantify glacier melt as a component of watershed discharge. In contrast to the well-studied crystalline-cored Cordillera Blanca in the outer tropics, Chimborazo is a stratovolcano located in the inner tropics and therefore experiences less-pronounced seasonality in precipitation (Kaser and Osmaston2002) and more persistent ablation due to higher humidity (Favier2004; Harpold and Brooks2018; Vuille et al.2003). Higher humidity can enhance ablation rates by increasing net longwave radiation and condensation (Harpold and Brooks2018). Most mixing model analyses of melt contributions in the outer tropics have been limited to the dry season, leaving wet season effects less well understood. In the inner tropics, coincident glacier melt and precipitation inputs throughout the year could lead to multiple processes simultaneously driving discharge variability that are difficult to disentangle. Furthermore, Andean volcanoes may feature fractured bedrock aquifers, which in other parts of the world have been found to support greater groundwater storage and baseflow than bedrock in crystalline-cored mountainous watersheds (Markovich et al.2016; Tague and Grant2009), adding another factor to be reconciled. A growing body of work at Volcán Antisana, also located in the inner tropics, has begun to shed light on its hydrogeologic (Caceres et al.2006; Cauvy-Fraunié et al.2013; Favier2004; Favier et al.2008) and ecohydrologic (Minaya2016) conditions, but comprehensive understanding of mountain hydrology in the inner tropics still greatly lags that in the outer tropics.

Here, we implement field and computational methods to answer the following two questions. (1) What is the temporal variability of relative glacier melt contributions to discharge, from hourly to multi-year timescales, in a subhumid glacierized watershed on Volcán Chimborazo? (2) What hydroclimatic factors control this variability? Our approach comprises three methods: mixing model analysis applied to repeat synoptic sampling, time series analysis of hydroclimatic data, and numerical watershed modeling. Each method interrogates a distinct temporal relationship, and synthesizing their results illuminates how the dominant surface and subsurface processes driving the hydrological response of a tropical glacierized watershed vary as a function of timescale.

2 Study area

Volcán Chimborazo is a glacierized stratovolcano in Ecuador (Fig. 1a) whose glaciers serve as the headwaters for four major river systems – the Río Mocha (NE flank), Río Colorado (NW flank), Río Guano (SE flank), and Río Chimborazo (SW flank) – that supply water to a population of over 200 000 (INEC2010). Located in the inner tropics, Chimborazo's climate is characterized by minimal intra-annual temperature variation ( 2 C) and moderately seasonal precipitation, with two wetter seasons of unequal length (February–May and October–November) (Clapperton1990) and two intervening drier seasons that have less but not negligible amounts of precipitation. Moisture mostly originates from the Amazon basin to the east (Smith et al.2008; Vuille and Keimig2004), which produces a steep northeast (up to 2000 mm yr−1) to southwest (<500 mm yr−1) precipitation gradient across the mountain (Clapperton1990). Driving interannual climatic variability at the regional scale, El Niño generally brings drier and hotter conditions throughout the Andes (Bradley et al.2003; Francou2003; Smith et al.2008; Vuille and Bradley2000; Vuille and Keimig2004; Wagnon et al.2001), which enhances glacier ablation (Favier2004; Veettil et al.2014b; Wagnon et al.2001).

Figure 1(a) Satellite image of Volcán Chimborazo, with the study watershed Gavilan Machay outlined in red, and its location in Ecuador shown in the inset map. The glacierized Gavilan Machay watershed is a relatively humid watershed compared to the western flank of Chimborazo. (b) Land cover and locations of monitoring stations and water sampling within the Gavilan Machay watershed.


Records since 1980 indicate that, consistent with the rest of the tropical Andes, temperatures have warmed 0.11 C decade−1 around Volcán Chimborazo (La Frenierre and Mark2017; Vuille et al.2008). This likely caused the 21 % reduction in ice surface area and 180 m increase in mean minimum elevation of clean ice observed between 1986 and 2013 (La Frenierre and Mark2017). Although regional precipitation gauges show no notable change over time, local residents report a reduction in precipitation, which could further drive glacier mass balance changes (La Frenierre and Mark2017). Historical records of glacier melt are not available. Under current conditions, only 4 of Chimborazo's 17 glaciers, including the 2 largest, Reschreiter (2.55 km2) and Hans Meyer (1.33 km2), generate perennial surface discharge, nearly all of which flows northeast into the Río Mocha watershed. The lowest 16 % of Reschreiter Glacier is debris covered, providing insulation that stabilizes ice at lower elevations (4480 m a.s.l.) than would be expected for clean ice, given current climatic conditions. Our study focuses on the 7.5 km2 Gavilan Machay sub-catchment on the subhumid northeast flank of Chimborazo (Fig. 1b), which is 34 % glacierized by Reschreiter and is of concern because it discharges into the main Río Mocha channel just upstream of the Boca Toma diversion point (3895 m a.s.l. elevation) for an irrigation system.

In addition to glacier melt, groundwater and ecological conditions also control the hydrology of the Gavilan Machay watershed. Springs are prevalent below 4400 m a.s.l. Geologic maps and stratigraphic interpretations (Barba et al.2005; Samaniego et al.2012) support field evidence for aquifers within unsorted glacial deposits and fractured bedrock (McLaughlin2017). Extensive areas of páramo, the biologically rich grasslands endemic to the tropical Andes above  3500 m a.s.l., are common across the watershed. Wet páramos commonly contain homogeneous Andosol soils of volcanic origin that can accumulate elevated organic carbon content; this typically gives rise to high porosity, infiltration capacity, hydraulic conductivity, and water retention (Buytaert and Beven2011; Buytaert et al.2006). Absorbent páramo soils are considered to very efficiently regulate watershed discharge throughout Andean Ecuador (Buytaert and Beven2011; Buytaert et al.2006; Minaya2016).

3 Methods

3.1 Hydroclimatic data

Precipitation, temperature, and relative humidity data were collected from October 2011 to February 2017 from weather stations installed at 4515 m a.s.l. on the debris-covered portion of Reschreiter Glacier and at 3895 m a.s.l. at Boca Toma (Fig. S1 and data files in the Supplement). Logistical obstacles prevented instrumentation of the upper three-quarters of the watershed above Reschreiter. The Boca Toma weather station was deployed with an Onset Hobo Pendant® Event Data Logger starting on 16 June 2015. The Reschreiter weather station was deployed with an Onset Hobo Micro Station starting October 2011, but temperature and precipitation data recovery was discontinuous. The short data records from Reschreiter were primarily used together with Boca Toma data to determine a lapse rate for precipitation. Precipitation in mountainous watersheds typically exhibits piecewise linear relationships with elevation, with positive (negative) lapse rates below (above) the elevation of maximum precipitation (Wang et al.2016). We calculated a negative lapse between the two weather stations, which we applied over the entire watershed with the assumption that the elevation of maximum precipitation lies below the watershed. Uncertainty in this approach arises from the actual unconstrained elevation of maximum precipitation (which requires more than two weather stations) and from unquantified precipitation measurement errors that may be caused by wind and freezing temperatures at high elevations. Glacier melt was separately estimated through model calibration to stream discharge data (Sect. 3.3), which may compensate for errors in the precipitation inputs. Temperature lapse rates were calculated using data collected at glacier ablation stakes (described further below) over June 2016 to November 2016. Relative humidity measured at the Boca Toma station was applied over the entire watershed due to the lack of measurements elsewhere. Discharge simulations should be less sensitive to errors in relative humidity compared to precipitation and temperature, which directly control water inputs to the watershed. We obtained unmeasured meteorological variables (wind speed, solar radiation, longwave radiation, and air pressure) from the Global Land Data Assimilation Systems (GLDAS) (Rodell et al.2004). A discharge gauging station equipped with a Solinst Levelogger Junior pressure transducer was established at Gavilan Machay, 1.1 km upstream of the Río Mocha confluence. Solinst Barologger measurements at Boca Toma were applied to correct for atmospheric pressure, and standard USGS rating curve techniques (Andrews1981) were used to convert water depth to discharge over the period of record (Fig. S1).

In June 2016, glacier ablation stakes were installed at two elevations on the Reschreiter Glacier tongue (4792 and 4820 m a.s.l.) and one on the Hans Meyer Glacier tongue (4925 m a.s.l.), all in clean ice. Each stake included a temperature sensor, along with a look-down ultrasonic sensor for measuring changes in distance to the ice surface to estimate glacier mass loss in clean ice. The stakes were deployed using the open-source Arduino-compatible ALog data logger (BottleLogger v1.4.0, an intermediate model between v1.0.0, Wickert2014, and v2.2.0, Wickert et al.2018). All sensors were mounted at the top of 3 m long PVC tubes, which were inserted into holes drilled to about 2.5 m depth. In addition to the clean-ice mass loss determined with ablation stakes, glacier volume change in debris-covered ice was estimated by differencing a GPS-validated photogrammetric digital elevation model (DEM) in 1997 (Jordan et al.2010) and terrestrial laser scanner (Riegl LMS-Z620) surveys in 2012 and 2013 (La Frenierre2016). Because of the sparse spatiotemporal coverage of these glacier melt measurements, these were used only as comparisons for calibrated glacier melt and were not directly applied in our analysis.

Hydroclimatic data collected in the watershed were directly assessed using statistical analyses and implemented as inputs to the integrated hydrologic model. For the statistical analysis, we calculated cross-correlations to probe how the discharge time series may be driven by different climatic factors. Spectral analysis provided further insight on the timescales, represented by time frequency, over which these interactions occur. Specifically, we examined the magnitude squared coherence (Cxy) over time frequency (f):

(1) C x y ( f ) = | S x y ( f ) | 2 S x x ( f ) S y y ( f ) ,

where Sxx(f) and Syy(f) are auto-spectral densities of variables x and y, respectively, and Sxy(f) is the cross spectral density of x and y. Like the square of a correlation, the magnitude squared coherence varies between 0 and 1, with 0 indicating the weakest relationship between the two variables at frequency f and 1 indicating the strongest relationship.

3.2 Hydrochemical and isotopic tracers

3.2.1 Field sampling

Water samples were collected for use in the Hydrochemical Basin Characterization Model (HBCM) (Sect. 3.2.3), a hydrochemical mixing model that spans the stream network and requires synoptic water sampling over a sufficiently short time period such that data reflect spatial and not temporal variability. We carried out five synoptic sampling campaigns during 1–8 January 2012, 7–9 July 2012, 12–15 June 2015, 25–30 June 2016, and 4–7 February 2017. The June and July (January and February) samples represent the longer (shorter) dry season. Dry seasons were targeted because of water resource interests during these periods; integrated hydrologic model simulations served to extend the analysis to wet seasons. In addition to limiting the number of days spanned during a sampling campaign, synoptic sampling should avoid hourly timescale hydrochemical fluctuations. All samples were collected between mid-morning and mid-afternoon. In February 2017, we confirmed that 1 min resolution-specific conductivity changes over a 24 h time period at the Reschreiter Glacier tongue were 1 order of magnitude smaller than the spatial variability across the Gavilan Machay subcatchment (McLaughlin2017). Logistical difficulties prevented similar measurements farther downstream in the watershed, where dynamic melt versus groundwater contributions likely caused greater hydrochemical variability (Sect. 4.1).

During each of the five campaigns, we collected water samples from meltwater (which may contain both glacier melt and snowmelt), springs, and precipitation, as well as at stream confluence mixing points (locations shown in Fig. 1b). Spring samples from concrete capture boxes or natural valley wall seeps represent groundwater, which consists of an unconstrained mix of shallow saturated soil water from páramo areas, morainic debris aquifer water, and deeper fractured bedrock aquifer water. Precipitation samples were collected using evaporation-proof totalizing rain gauges deployed for 3–6 days at Boca Toma, near the Reschreiter weather station, and near Hans Meyer Glacier (at 4780 m a.s.l.). Each field campaign covered most of the same sampling locations between the Reschreiter Glacier tongue and the Gavilan Machay confluence. The 2012 and 2017 sampling periods included additional stream samples between some confluences to estimate groundwater contributions along shorter stream reaches (Fig. 2b). For each sampling site, 30 mL of water was collected, filtered in the field using either 0.45 µm (before 2017) or 0.2 µm (2017) filters, and stored in Nalgene bottles that were capped and sealed with electrical tape, and then stored near 4 C as soon as possible.

Figure 2Different computational cell configurations for the HBCM mixing model based on the available samples (site codes in squares; see Fig. 1) for each of the five periods. The upper solid line for each period represents the main channel. The lower solid lines depict tributary links for confluence cells, and the lower dotted lines show groundwater inputs to the main channel for reach cells.


3.2.2 Laboratory analysis

In 2012, major dissolved ions (Li+, Na+, K+, Mg2+, Ca2+, F, Cl, NO3-, PO43-, and SO42-) were measured using a Dionex DX-500 ion chromatography system at the Water Isotope and Nutrient Laboratory at the Ohio State University, and stable isotopes of water (δ18O and δ2H) were measured using Picarro L2130-i cavity ring-down spectroscopy (CRDS) isotope analyzers at the Water Isotope and Nutrient Laboratory and at the Byrd Polar and Climate Research Center. In 2015–2017, cations (Na+, K+, Mg2+, Ca2+) were measured using an Agilent 7700X inductively coupled plasma mass spectrometer (ICP-MS) at Gustavus Adolphus College, anions (F, Cl, SO42-) were measured using a Dionex ICS1000 ion chromatographer also at Gustavus Adolphus College, and stable isotopes of water (δ18O and δ2H) were measured at the University of Minnesota using an LGR DLT-100 liquid water analyzer (a laser spectroscopy system). We calculated the bicarbonate (HCO3-) concentration as the charge balance residual. Reported isotope ratios are relative to the Vienna Standard Mean Ocean Water (VSMOW) and typical precisions are ±1.0 ‰ for deuterium/hydrogen values and ±0.25 ‰ for 18O∕16O values. Checking for consistency across instruments, we confirmed that bulk concentrations at each location and spatial trends for each analyte were similar across sampling periods (Fig. S3). Certain analytes did exhibit a discernible systematic bias for a particular sampling period, which may be related to laboratory instrument, but this should not pose a problem for the mixing model, because it is implemented only using data within the same sampling period (measured on the same instrument). All hydrochemical data are available in the Supplement.

3.2.3 Mixing model: Hydrochemical Basin Characterization Model (HBCM)

Naturally occurring dissolved ions and stable isotopes of water (δ18O and δ2H) are commonly used to track the relative contributions of different surface source waters to total watershed discharge (Baraer et al.2009; Hooper and Shoemaker1986; Mark and Mckenzie2007; Mark and Seltzer2003; Ryu et al.2007), as well as to identify groundwater flow paths (Baraer et al.2009, 2015; Clow et al.2003; Crossman et al.2011; Kendall et al.2003). Here, the proportion of glacier and snowmelt versus groundwater in discharge at the Gavilan Machay watershed is quantified using HBCM, a multicomponent hydrochemical mixing model developed for use in data-sparse, glacierized tropical watersheds (Baraer et al.2009). Given source (or end-member) and outflow chemistries at different mixing points throughout the watershed, HBCM solves an over-constrained set of mass balance equations for multiple tracers to determine the relative flow contributions of each source. Details are provided in the Supplement (HBCM section).

As with all hydrochemical mixing models, HBCM's calculation of relative source contributions depends on three fundamental assumptions: (1) end-member chemistry is unique and spatially homogeneous within the analysis area, (2) tracers are chemically conservative within the analysis, and (3) end-member mixing is instantaneous and complete (Christophersen et al.1990; Soulsby et al.2003). A unique feature of HBCM is that it represents spatial information within a watershed through a series of cells that are interconnected by having outflow from one become inflow to a subsequent downstream cell. There are two types of cells, both of which have streamflow at the upgradient end of the cell as a source and streamflow at the downgradient end of the cell as the mixed water output. “Reach cells” have groundwater as their other source, and “confluence cells” have tributary water as their other source (see illustration in Fig. S2). Note that this makes end-members for a particular HBCM cell different than for the full watershed, which has only meltwater and groundwater as sources contributing to discharge at the outlet. Figure 2b shows the conceptual schematics of the Gavilan Machay cell configuration for the five sampling periods.

Although HBCM only requires the three assumptions to be met on a cell scale, we carried out a preliminary watershed-level analysis considering groundwater and meltwater as sources to all mixed stream samples, in order to identify potential conservative tracers that are reasonable candidates for all cells. Appropriate tracers should show end-members appearing on opposite ends of a line formed by mixed samples in bivariate plots, and samples for different end-members should group separately from each other in hierarchical cluster analysis diagrams (Christophersen et al.1990; Hooper2003; James and Roulet2006). Stable isotopes were excluded as tracers for reach cells, because the groundwater likely has a range of isotopic values due to different recharge elevations.

To bracket some of the uncertainty in the method, HBCM generates estimates of fractional contributions from each source using different combinations of potential tracers. The final result consists of a range of estimates that produce similar (within about 3 times the minimum) cumulative residual errors between the measured tracer concentrations in the mixed outflow water and that predicted by the over-constrained mixing model. This quantifies uncertainty due to the model's inability to distinguish among equally good optimization results but represents only a lower limit of error, because it does not account for the mismatch between the observed and predicted mixed concentration outflows. There are no straightforward methods to convert the sum of residual concentration flux errors to estimated source contribution errors.

3.3 Integrated hydrologic modeling

Spatially distributed watershed models can integrate surface hydrology and groundwater flow through time to evaluate their joint impacts on water resources. Over the 1-year period of June 2015–June 2016 when continuous air temperature and precipitation measurements are available in the watershed, we implemented Flux-PIHM version 0.5.0 (Shi et al.2013), an intermediate complexity watershed model that balances mechanistic parameterizations with computational efficiency. Full details about Flux-PIHM can be found in Qu and Duffy (2007) and Shi et al. (2013); here, we summarize the major features. Flux-PIHM couples physically based equations for canopy interception, infiltration, surface and subsurface water flow, and snowmelt with the energy balance scheme of the Noah land-surface model (LSM) (Ek et al.2003) for more accurate simulation of evapotranspiration. Flux-PIHM employs a semi-discrete finite volume approach on an unstructured grid that performs efficiently on steep topographies. Channel and overland flow are represented by diffusion wave approximations to St. Venant equations, shallow groundwater flow follows a 2-D Dupuit approximation, and unsaturated zone flow is based on a 1-D form of the Richards equation. The model simulates water storage in one vertically integrated unsaturated zone layer and one vertically integrated saturated zone layer, providing a “2.5-D” distributed model. Other PIHM family codes include the watershed reactive transport module RT-Flux-PIHM (Bao et al.2017; Li et al.2017), as well as the hydrologic and landscape evolution module LE-PIHM (Zhang et al.2016), among a suite of other functional modules (Duffy et al.2014).

Flux-PIHM determines snowmelt based on energy balance. Use of coarse-scale GLDAS radiation inputs introduces errors, but as will be discussed in Sect. 4.3.1, precipitation limitations may make snowmelt calculations less sensitive to radiation input uncertainties compared to glacier melt calculations. Due to the unavailability of high-resolution radiation input measurements as well as intensive source-code modifications required to couple energy balance calculations for ice melt into the Flux-PIHM, we added a separate module to simulate glacier melt using a temperature-index scheme (NRCS2009). Although the accuracy of a temperature-index glacier melt model for tropical glaciers can be uncertain due to uncaptured effects of solar radiation, cloud cover, humidity, topography, and aspect (Fernández and Mark2016; Gabbi et al.2014; Hock1999, 2005; Huss et al.2009; Pellicciotti et al.2005; Sicart et al.2008), it remains the most feasible approach in poorly instrumented watersheds given its simplicity and limited field data requirement compared to an energy balance approach (Fernández and Mark2016; Hock2005; Reveillet et al.2017). The temperature-index glacier melt model includes


where FI is the ice melt rate (m h−1), Ta is air temperature in the grid cell containing ice (C), MI is the melt factor parameter (m h−1C−1) to be calibrated, and TM,I (C) is set to 0 C as the air temperature threshold for ice to melt. Over the simulated time period, we assume that there is an unexhausted supply of ice that can melt in the glacier-covered grid cells below the equilibrium line altitude (ELA), located at  5050 m a.s.l. (Frenierre and Mark2014), which is a reasonable approximation over the 1-year simulation period. The melt simulated with the temperature-index model was added to the precipitation amount for the Flux-PIHM forcing inputs.

We used the PIHMgis software (Bhatt et al.2014) to construct an unstructured domain of 188 cells over the Gavilan Machay subcatchment using a 30 m resolution Shuttle Radar Topography Mission (SRTM) digital elevation model (Farr et al.2007). Although a major feature of PIHMgis is its tight integration with spatial and temporal datasets for model inputs such as soil properties and meteorological forcing, these datasets only cover densely monitored regions, mostly within North America and Europe. For meteorological forcing, we used the spatially distributed inputs described in Sect. 3.1. Vegetation mapping by McLaughlin (2017), based on 30 cm resolution aerial photo surveys conducted by the Sistema Nacional de Información de Tierras Rurales e Infrastructura Tecnológica (SIGTIERRAS; last access: 2 January 2019), provided land-cover types and boundaries. Built-in land-cover parameters from the Noah LSM were used for the “grassland/herbaceous” type at the lowest elevations corresponding to páramo, the “barren/sparsely vegetated” type for intermediate elevations with rock/dirt/gravel, and the “perennial ice/snow” type for the ice-covered areas. This approach simplifies the mix of tussock grasses, acaulescent rosettes, and cushion plants that make up the páramo into a single representative grassland/herbaceous type in order to reduce the calibration burden. For the grassland/herbaceous land-cover type, the default monthly leaf area index (LAI) values were replaced with measurements from MODIS (Vermote2015) to avoid using incorrect seasonal changes from the original model settings for this tropical region. Hydraulic parameters were manually calibrated to match observed discharge at Gavilan Machay.

4 Results and discussion

This section presents the respective insights gained from each of the three methods on the temporal relationship among meltwater, groundwater, and discharge: the mixing model analysis offers discrete multi-year estimates over 5 years, the time series analysis shows fine-scale hourly resolution correlations over nearly a year with continuous hydroclimatic observations, and the integrated hydrological model explores intermediary weekly to seasonal processes within a 1-year simulation period containing a strong El Niño event. A complete interpretation of the multi-scale temporal variabilities and their hydroclimatic controls emerges in Sect. 4.3.2 when evaluating the model simulations in relation to the mixing model and time series analysis results.

4.1 Mixing-model analysis of meltwater contributions to discharge

Total cation concentrations provide a summary representation of hydrochemistry results from the five dry-season synoptic sampling periods in Fig. 3. These plots show that, even though hydrochemical conditions vary over the different periods, groundwater samples, which geochemically interact with soil and rocks, consistently contain much higher ion concentrations than meltwater samples. The distinctive chemistries of groundwater and meltwater make it possible to use the mixing model approach to estimate their relative contributions to streamwater, which shows an increase in ion concentration while moving downgradient due to the cumulative addition of groundwater (see Fig. 1b for sample locations). We chose as tracers those analytes that most consistently showed the mixed sample visually falling close to the line between its two source samples in the bivariate plots in Figs. S4–S8: sum of monovalent cations, Mg2+, Ca2+, Cl, and HCO3-. Hierarchical cluster analysis lends confidence that these five potential tracers can be used in the mixing model analysis to distinguish between groundwater and melt samples as different watershed-level end-members (Fig. S9). Precipitation was not included as an end-member, because precipitation samples tended to plot outside the range of stream samples bracketed tightly by groundwater and meltwater samples in bivariate plots. Rather than directly add to streamflow, most of the precipitation that fell close to the sampling time likely evapotranspired or infiltrated and contributed to streamflow through groundwater.

Figure 3(a) Hydrochemistry at different sampling locations along the Gavilan Machay main channel shows variability in space and over the five sampling periods; site codes are ordered from highest (0 for meltwater) to lowest (16 for Boca Toma) elevation (see Fig. 1b). (b) Hydrochemistry for different tributary (sites 3, 6, 9, and 15) and groundwater (spring) samples (sites 5 and 13), ordered from highest to lowest elevation, shows variability in space and time.


HBCM results in Fig. 4 illustrate the importance of both meltwater and groundwater in the watershed. Surficial meltwater comprises between 23 % and 66 % of discharge during the five dry-season sampling periods, with groundwater constituting the remaining 34 %–77 % at any given time. Notable differences were observed across the sampling periods. The higher relative melt contribution during February 2017 compared to January 2012 could reflect the accelerating melt rates observed on Chimborazo (La Frenierre and Mark2017). However, the absolute melt contribution, determined by applying estimates of relative melt contributions to average observed weekly discharge measurements around the sampling time, was in fact lowest in February 2017, because of significantly lower total discharge compared to the other sampling periods (Fig. 4b). The lower total discharge was likely due to lower precipitation and temperature during the weeks around the sampling period compared to during the other sampling periods (Fig. S1). Our findings across the five sampling periods demonstrate that one single synoptic tracer test should not be directly generalized or interpreted without considering temporal dynamics and groundwater conditions.

Figure 4The HBCM mixing model predicts a range in relative surficial meltwater contributions to discharge over five discrete sampling times, which may reflect both temporal changes and uncertainties. Error bars bracket HBCM estimates that produced similar best matches to observed tracer concentrations; however, actual uncertainties are higher because of residual errors. Absolute meltwater discharge contributions can vary in time very differently than relative inputs, in part due to varying groundwater contributions.


Our results show Volcán Chimborazo to deviate from trends found at the well-studied Cordillera Blanca, likely due to its distinct climatic and geologic conditions. When compared to an exponential fit between relative groundwater contribution and glacierized fraction for four watersheds in the Cordillera Blanca (Baraer et al.2015), our estimates for groundwater contributions in Gavilan Machay are approximately twice as large. Also, the glacierized Gavilan Machay sub-catchment of the upper Río Mocha watershed has a specific discharge that is less than half of that in the non-glacierized portion, in contrast to the greater specific discharge generally found with greater glacierized areas in the Cordillera Blanca (Baraer et al.2009; Mark and Seltzer2003).

HBCM implementation with five different field campaigns enabled us to evaluate uncertainties due to distinct sampling plans (details in the Supplement – Table S1). We found that changing HBCM cell configurations could generate up to a 23 % melt fraction difference in estimates. Also, having fewer HBCM analysis cells (e.g., longer stream reaches) and fewer groundwater samples consistently led to greater HBCM residual errors. Too few groundwater samples becomes problematic when groundwater is not a homogeneous end-member throughout the watershed, which is the case in Gavilan Machay, which contains springs with somewhat higher solute concentrations at lower elevations (Fig. 3b). Further, errors in the estimated groundwater contribution grow when using fewer and longer reach cells, because with additional and shorter reach cells, observations can reset the stream channel chemistry to correct concentrations. These results demonstrate the importance of adequately measuring the spatial variability of the surface and subsurface flow network, and they prompt the use of alternative methods to help constrain uncertainties in HBCM analysis results.

4.2 Time series and spectral analysis of hydroclimatic controls

Because of the uncertainties and long time gaps in the HBCM analyses, we applied statistical analyses to the continuous data available from July 2015 to March 2016 at the Boca Toma weather station and Gavilan Machay gauging station to further infer characteristic trends between meltwater and discharge and their climatic controls. Considering air temperature as a proxy indicator of meltwater, the hourly cross-correlation of air temperature leading discharge at Gavilan Machay in Fig. 5a shows a strong diurnal signal, with peak discharge occurring 4 h after the warmest part of the day at an average rate of 0.1 m3 s−1 that is about twice the magnitude of average morning discharge.

Figure 5Cross-correlations (with 95 % confidence interval shown) of observed discharge at Gavilan Machay with (a) hourly temperature and (b) weekly precipitation show that discharge has a clear diurnal link with temperature at about a 4 h lag and a strong relationship with weekly precipitation, respectively. Magnitude squared coherence (with various confidence intervals shown) between discharge and (c) temperature exhibits a high peak at 24 h corresponding to the cross-correlation result, as well as a strong peak at 12 h, and a more moderate peak at a multi-day scale. The coherence between discharge and (d) precipitation peaks between 100 and 200 h (about 1 week) and may also be significant at scales approaching 1 year.


To determine if melt could be driving discharge variability beyond diurnal timescales, we examined the magnitude squared coherence between temperature and discharge in Fig. 5c, which quantifies their correlation at a certain period (inverse time frequency). As expected from the time series results, the most prominent feature is the peak at a period of 24 h. Interestingly, the next highest coherence is at a period of 12 h. While this could be a harmonic artifact of the dominant 24 h trend, it is supported by a slight temperature increase commonly observed around 23:00. The resulting increase in discharge, while detectable, is inconsequential to the total daily water balance. Other very narrow coherence peaks at periods less than 12 h are likely spurious, because the power spectral densities of both temperature and discharge are low over that range (Fig. S10). However, the smaller but broader peak around 50 to 60 h suggests that multi-day warming may also drive multi-day discharge events, though this link is much weaker than the diurnal response.

The substantial groundwater contributions to discharge inferred from the HBCM analysis prompts a look at not only melt but also precipitation controls on multi-day discharge. Hourly precipitation and discharge are very weakly correlated (Fig. S11a); however, a significant correlation of 0.5 appears for weekly averages with zero time lag (Fig. 5b). Correspondingly, a high (above 95 % confidence interval) and broad coherence peak between precipitation and discharge can be seen over periods of about 1 week (168 h) (Fig. 5d). Together, these results suggest that sustained rain events influence discharge over 1 week and therefore that rainwater tends to infiltrate instead of flow quickly overland. Other statistically significant (at the 95 % confidence interval) coherences between discharge and both temperature and precipitation across multi-week to multi-month periods further support the role of even slower subsurface flow pathways.

Combining the time series analysis with the HBCM results suggests that streamflow at Gavilan Machay is heavily influenced by both surficial meltwater and groundwater, and that the latter is mostly driven by precipitation. Furthermore, time series and spectral analyses highlight temporal links not easily found through a fieldwork-intensive tracer-based approach: meltwater feeds discharge at Gavilan Machay on an hourly timescale, while weekly discharge responds most strongly to precipitation events. There are, however, limitations to this statistical assessment that result from the short 9-month dataset, which precludes robust examination of any seasonal to multi-year responses to bimodal wet and dry seasons and El Niño effects. Hydrologic modeling in the following section can address this, as well as questions about the quantitative role of melt contributions or groundwater buffering periods of low rainfall.

4.3 Integrated hydrologic model simulations

4.3.1 Calibration results

Matching the observed discharge dynamics with the hydrologic model required calibration of two different melt factors based on time period: a lower value of 7.10 mm w.e. (millimeters water equivalent) C−1 d−1 over December 2015–February 2016 and a higher value of 8.64 mm w.e. C−1 d−1 over the rest of the simulation. The estimated melt factors fall within the range of melt factors calculated for other glaciers in the tropics (3.5–9.9 mm w.e. C−1 d−1) reported in Fernández and Mark (2016). Our simulation period coincides with a strong El Niño event that generated the warmest and driest conditions from late November 2015 to the start of February 2016 (Fig. 6a–b). The lower melt factor at this time dampens the intensity of glacier melt, possibly because of the absence of heat transfer from rain (Francou2004), but it nonetheless simulates among the highest glacier melt volumes of the simulation period (Fig. 6a), consistent with other studies showing increased glacier melt in response to El Niño events in the Andes (Francou2004; Manciati et al.2014; Maussion et al.2015; Veettil et al.2014a, 2016). Huss et al. (2009) similarly estimated a lower melt factor in the Swiss Alps under warmer conditions, which they attributed to warmer conditions being driven by longwave rather than shortwave radiation. Overall, in Gavilan Machay, the average specific glacier melt rate (in water equivalence) simulated over glacierized areas below the ELA was 1.5 m yr−1. This falls within the range measured at the Reschreiter Glacier tongue, bracketed by mass balance estimates on slower-melting debris-covered ice of 0.87 m yr−1 (1997–2013) and 0.54 m yr−1 (June 2012–January 2013), as well as average ablation stake observations on faster-melting clean ice of 3.4 m yr−1 (June–November 2016). Although useful for comparing against the calibrated melt model, these measurements do not cover sufficient areas and periods of time to constrain separate melt factors for debris-covered and clean ice.

Over the entire watershed, the resulting calibrated glacier melt production is equivalent to 68 % of the precipitation input and 567 % of the simulated snowmelt amount during the simulation period. Based on temperature, Flux-PIHM partitioned the precipitation input into 12 % snowmelt and 88 % rainfall. The much smaller amount of simulated snowmelt compared to glacier melt supports the earlier suggestion (Sect. 3.3) that snowmelt could be precipitation limited rather than energy limited. This helps justify our separate approach of simulating snowmelt through Flux-PIHM's energy balance module while simulating glacier melt through a calibrated temperature-index model, because energy balance calculations of glacier melt would be much more sensitive to uncertainties in coarse-scale radiation inputs than snowmelt. We do acknowledge, however, that snowmelt simulations depend on lapse-rate-determined precipitation inputs that have their own uncertainties (Sect. 3), and so our calibrated glacier melt may include some amount of snowmelt that is not represented in the model.

Figure 6Time series of weekly moving average of (a) average air temperature over the ablation zone (glacier-covered areas below the ELA (5050 m a.s.l.) and simulated glacier melt production; (b) precipitation (solid line) and precipitation+glacier melt production (dashed lines) (c) discharge at Gavilan Machay from observations (gray), calibrated simulations (blue solid), and simulations with no ice melt (black solid); groundwater contribution to discharge for calibrated (blue dashed) and for no-ice simulations (black dashed); and (d) simulated percent glacier melt contribution to discharge at Gavilan Machay. The shaded block shows the range of percent melt contributions estimated for five sampling periods over 2012–2017 from the mixing model; these estimates may include snowmelt in addition to glacier melt, but omit meltwater contributions through groundwater.


The calibration procedure also involved soil parameter adjustments. Hydraulic parameter estimates in páramo environments are scarce, and their characterization can be uncertain (Buytaert et al.2006). For an initial estimate, we applied pedotransfer functions used in Flux-PIHM (Wosten et al.1999) to a range of páramo soil measurements from a study area 20 km northwest of Chimborazo (3800 to 4200 m a.s.l.) (Podwojewski et al.2002) and from a study watershed on the glacierized Volcán Antisana also in the Ecuadorian Andes (4000–4600 m a.s.l.) (Minaya2016) (see Table S2 for full details). We then calibrated the model for three mapped land-cover zones corresponding to páramo, rock/dirt/gravel, and ice (Fig. 1). In the páramo zone (Table 1), matching observed discharge required lower hydraulic conductivity and greater water retention (expressed in van Genuchten hydraulic parameters) than initially estimated. This is likely due to high organic matter content supported by the study area's subhumid conditions and the well-recognized retentive hydraulic properties of páramo soils (Buytaert et al.2006; Podwojewski et al.2002). In the sparsely vegetated and ice-covered zones, the calibration yielded higher hydraulic conductivities and lower water retentions than in the páramo zone, corresponding to reduced organic matter fraction and fractured bedrock, though the hydraulic conductivities were still lower than the initial estimates from the pedotransfer functions.

Table 1Soil hydraulic parameters calibrated for the three soil-type areas compared against the range predicted using the pedotransfer function with measured páramo soil textures in Ecuador (Minaya2016; Podwojewski et al.2002). Parameters include hydraulic conductivities for vertical infiltration (KINFV), vertical saturated zone flow (KSATV), and horizontal saturated zone flow (KSATH); porosity; residual soil moisture (θres); and shape parameters (α and β) for the van Genuchten moisture retention curve: θ=θres+porosity×11+|αψ|β(1-1β), with water content θ and pressure head ψ.

a Predicted using pedotransfer functions with nine páramo soil texture measurements from Ecuador. Samples are from three locations in a glacierized watershed in Volcán Antisana (Minaya2016) and four locations about 20 km northwest of Chimborazo (Podwojewski et al.2002) (see Supplement for more details). b Calibrated to match observed discharge at Gavilan Machay.

Download Print Version | Download XLSX

Simulation results in Fig. 6c show that the calibrated model parameters closely produced the observed weekly discharge, including lower discharge under the drier and warmer El Niño conditions in December 2015–January 2016. The single major model mismatch occurred at a precipitation-driven discharge peak in February 2016. This could reflect uncertainties from the soil parameter calibration, as well as from our use of a lapse rate-based precipitation field over complex terrain, in which high-altitude precipitation events may not all be recorded at the low-altitude rain gauge. On shorter timescales, hour-of-day simulation results in Fig. 7b demonstrate that the model does produce a diurnal trend, but with slightly less than half the average range and at a 6 h later peak compared to observations (Fig. 7a). These hourly discrepancies can be attributed to weaknesses in the simple melt model. Hock (2005) argued that temperature-index models can successfully capture seasonal glacier melt trends but struggle with diurnal fluctuations, which are strongly driven by solar radiation dynamics. Although our simulations cannot reliably produce the timing of hourly discharge, they can provide informative lower bounds on the size of the average diurnal range.

Figure 7Box plots over local hour of day showing the 25th to 75th percentiles with boxes and maximum and minimum with whiskers for (a) observed discharge over the 9-month observation period, (b) simulated discharge over the 9-month observation period, (c) simulated discharge over the 12-month simulation period, and (d) simulated percent glacier melt contribution to discharge over the 12-month period. Simulations do capture diurnal patterns, but with an underestimated magnitude and shifted peak. Simulated percent glacier melt contributions to discharge closely mirror the simulated diurnal fluctuations in discharge.


4.3.2 Simulations of relative glacier melt contribution to discharge

To quantify relative glacial meltwater contribution to stream discharge using Flux-PIHM, we compared the calibrated discharge simulations at Gavilan Machay (QCalib) with simulations that omitted glacier meltwater in the forcing inputs (QNoGlacierMelt). We then calculated relative glacial meltwater contribution to discharge via

(3) % Glacier Melt = Q Calib - Q NoGlacierMelt Q Calib × 100 % .

Apart from the addition of glacier melt in the calibrated case, the two model scenarios include the same model inputs, including air temperature, precipitation (including identical snowmelt), land cover, and hydrologic properties. Thus, our calculation of change between the two scenarios isolates the effect of having glacier melt versus not having glacier melt. Overall, over the 1-year simulation period, an average of 52 % of stream discharge in Gavilan Machay can be attributed to glacier melt, compared to an 8 % contribution by snowmelt. Because the simulated glacier melt amount was calibrated in addition to precipitation inputs, we consider this glacier melt contribution to originate from the preexisting ice reservoir at the start of the model period. However, as noted above (Sect. 4.3.1), the calibrated glacier melt contribution could include some amount of snowmelt not fully accounted for in the model due to uncertainties in precipitation inputs.

The estimate provided by Eq. (3) may not exactly correspond to the HBCM result for meltwater contribution for a number of reasons. First, the water balance impact of glacier melt conveyed in Eq. (3) may not equal the proportion of meltwater in a sample of discharge water if, for example, melt inputs facilitate more runoff of precipitation-sourced water. Also, while Eq. (3) aims to isolate glacier melt of preexisting glacier ice, HBCM estimates could include snowmelt and melt of freshly accumulated ice, depending on the composition of the meltwater sample taken just below the glacier tongue. Lastly, any melt that infiltrates is considered as part of the groundwater rather than melt contribution in HBCM estimates, while Eq. (3) includes the effect of both surficial and groundwater contributions of meltwater to discharge.

These conceptual discrepancies complicate comparisons between the two methods, but it is possible to assess the 20–25 June 2016 HBCM result against the simulation (other HBCM periods fall outside the 1-year model period). During 20–25 June 2016, the 45 %–64 % estimate with HBCM is higher than the average simulated relative melt contribution of 29 %, but it falls within the simulated hourly range of 13 %–70 % over that time. Considering that samples were collected during the daytime when melt contributions were generally high, the results from these two approaches reasonably agree.

While process-driven temporal patterns were difficult to glean from the sparsely spaced and uncertain HBCM results, Flux-PIHM simulations in Fig. 6d indicate considerable variability in weekly glacier melt contributions of 25 %–61 % over June 2015–June 2016. This range compares very closely with the 23 %–66 % range bounded by the five HBCM estimates spanning January 2012 to February 2017 (indicated by the shading in Fig. 6d), lending further confidence in the consistency between the watershed and mixing model results, despite the differences in temporal and other types of representation noted above. Hourly distributions of simulated glacier melt contribution in Fig. 7d show an average diurnal range of 25 %–68 %; given that the actual range is likely even broader due to underestimations by the temperature-index model, diurnal fluctuations in relative melt contributions may be of similar or even greater magnitude than changes on weekly to monthly timescales.

Comparing Fig. 6c and d reveals a remarkably strong correlation (0.89) between simulated weekly discharge and relative glacier melt contributions, indicating that major discharge peaks over multi-day timescales are melt-driven, a timescale link that was masked in the time series analysis in Sect. 4.2 by the overwhelming diurnal signal between temperature and discharge. Figure 8b shows strong coherence over 30- to 80-day periods between simulated glacier melt production and relative glacier melt contributions to discharge, highlighting those timescales as the most prominent for direct glacier melt inputs beyond diurnal timescales. Although occurring at a longer multi-month timescale, these melt contributions to discharge nonetheless appear to happen mostly through fast surface runoff; the model showed most (86 %) of the glacier melt contribution to discharge occurring through surface runoff processes, consistent with the hourly correlations found between observed discharge and temperature in Sect. 4.2. With the model, the influence of glacier melt production on discharge is clear during times such as early January 2016, when an uptick in discharge occurs despite low precipitation, due to warm temperatures and greater melt production and contribution (Fig. 6).

Figure 8Magnitude squared coherences between simulated percent glacier melt contribution to discharge with (a) precipitation and (b) simulated glacier melt production. Calculations used daily average data to avoid known uncertainties in modeling diurnal fluctuations of glacier melt production. Precipitation and glacier melt production appear to be related to glacier melt contributions at complementary timescales.


Unexpectedly, the peak relative glacier melt contributions (Fig. 6d) do not always align with glacier melt production patterns (Fig. 6a) (e.g., in early February 2017, percent glacier melt contribution peaks just when temperature and glacier melt production dip), and their weekly correlation is not strong at a 95 % confidence level (Fig. S11b). Squared coherences in Fig. 8a indicate that this is because relative glacier melt contribution also responds to precipitation, and this happens at weekly and multi-month timescales when correlations with glacier melt production are low (Fig. 8b). In fact, the strong coherence between precipitation and glacier melt contribution around 6- to 14-day periods reveals that the high correlation between observed weekly precipitation and discharge in Sect. 4.3 may relate in part to glacier melt contributions. Inspection of time periods such as mid-February and the start of April 2016 shows that rainy periods can augment glacier melt contributions to discharge (Fig. 6b, d). During those times, increases in relative glacier melt contributions to discharge coincided with precipitation events during local drops in glacier melt production, possibly because week-long precipitation creates antecedent moisture conditions that enhance the fraction of glacier meltwater that runs off over the surface. Surface runoff can be expected to contain mostly glacier meltwater or high-elevation precipitation, because runoff generally does not occur on low-elevation páramo soils except under intense precipitation (Harden2004; Sarmiento2000).

At greater timescales, significant coherences between precipitation and simulated relative glacier melt contributions around 120+ days (Fig. 8a) correspond to the overall lowest relative glacier melt contributions during the dry El Niño period (December 2015–February 2016) and highest relative contributions during the two wet seasons (February–May and October–November), possibly due to a transfer of heat from rain to ice. These glacier melt contribution responses to bi-seasonal to interannual climatic patterns likely occur through slower groundwater pathways; the simulation with glacier melt produces 18 % greater groundwater contributions to discharge than that without glacier melt (Fig. 6c). This supports the hypothesis that some amount of glacier meltwater infiltrates and recharges to groundwater that then discharges farther downgradient to the stream channel. It should be noted, however, that the subsurface component of Flux-PIHM contains only a single shallow saturated zone, which mostly closely corresponds to a perched water table in the soil zone. This model implementation may not accurately represent additional deeper fractured bedrock aquifer systems, which could discharge groundwater to streams over shorter, multi-month timescales (Andermann et al.2012), in contrast to the relatively constant groundwater dynamics simulated in Fig. 6c.

Our results may be particular to the moist climate in the Gavilan Machay watershed, which supports constant groundwater discharge to the stream throughout the year and reflects compounded effects of melt inputs and rainfall conditions. However, although glacier melt intensifies discharge variability in this watershed, a baseline 25 % glacier melt contribution throughout the simulation period indicates that a constant minimum level of glacier melt helps prevent episodes of even more extreme low flows during drought times, such as during El Niño conditions.

5 Summary and conclusions

Although meltwater is typically credited with modulating stream discharge by buffering periods of low precipitation, we demonstrate using a combination of methods that relative meltwater contributions may drive nearly all the variability in discharge (correlation of 0.89) over a range of hourly to multi-year timescales while buffering only against extreme low-discharge periods in a subhumid glacierized watershed in the Ecuadorian Andes. Hydrochemical mixing model results for five sampling periods spanning 2012–2017 showed the meltwater fraction in discharge may have varied over approximately 20 %–65 %. Hydrologic model simulations over June 2015–June 2016 produced a nearly identical range for weekly contributions. The model also predicted a very similar average diurnal range, which likely provided a lower bound of actual variability based on hydroclimatic data.

This multi-scale variability in melt contributions can be attributed to dynamic climate forcings that also contain a range of temporal patterns (Fig. 9). We found a strong correlation between diurnal temperature and discharge changes that likely reflects melt production and supports the use of a temperature-index melt model (Fig. 9a). Although such a simple melt approach somewhat underestimated hourly fluctuation extremes with a lag, it led to reasonable weekly discharge predictions when implemented with a seasonally variable melt factor that possibly accounts for additional heat transfer from rainfall to ice during wet seasons. Spectral coherence analysis of the model results showed that not only were diurnal discharge patterns responding to melt production, but relative melt contribution and discharge variations over 30–80-day periods were controlled by extended glacier melt production periods that also contribute to discharge through surface runoff (Fig. 9c).

Figure 9Relative melt contributions drive nearly all the variability in discharge in Gavilan Machay, mostly through surface runoff of glacial meltwater. What drives the variability in relative melt contributions to discharge? Our results show that this depends on the timescale. (a) Hourly timescale variability is controlled by radiation-driven (red arrows) melt production (light blue slab at upper left), which readily runs off overland (thick blue arrow) and eventually reaches the watershed's discharge point (circle with cross). (b) Weekly timescale variability is controlled by weekly precipitation events, which likely generate antecedent moisture conditions (light blue shading) that promote greater meltwater runoff. (c) Monthly timescale variability is driven by monthly trends in melt generation, which contributes to discharge mostly through surface runoff. (d) Seasonal to interannual variability is driven by long-term precipitation, which can enhance melt by transferring heat from rain (blue arrows to the right of the glacier tongue) and augment subsurface melt contributions through increased groundwater flow (thick blue arrow below water table).


Unexpectedly, model results showed that precipitation also boosted melt contributions to discharge, but on weekly and semiannual timescales that complement the hourly and monthly timescales controlled by temperature and melt production. Weekly precipitation events likely generate antecedent wet conditions that facilitate greater amounts of meltwater runoff (Fig. 9b), while longer-term precipitation patterns appear to drive slow changes in melt additions to groundwater (Fig. 9d). Most (86 %) melt contributions to discharge occurred through surface runoff in the model, but some meltwater recharged to groundwater, helping to support a relatively steady groundwater discharge to the stream that is about 18 % greater with glacier melt than without glacier melt. As expected, strong El Niño conditions corresponded to some of the highest simulated melt inputs, but less easy to predict was that Gavilan Machay exhibited its lowest discharge during this time. Melt prevented streamflow from dropping below a baseline level during the warm and dry El Niño event, but discharge was much higher during wetter periods, in part because of the rainfall-enhanced melt contributions described above.

In glacierized watersheds, slower glacier melt contributions through groundwater are poorly constrained. Past studies have identified a component of meltwater in the groundwater (Baraer et al.2015; Favier et al.2008; Harrington et al.2018; Lowry et al.2010; Minaya2016), but, to our knowledge, our work is the first to quantify this component. Generally, fractured young volcanic bedrock can support extensive groundwater (Frisbee et al.2011; Markovich et al.2016; Tague et al.2008), and, thus, other volcanic systems may contain similar meltwater fractions as found here. Meltwater–groundwater interactions may be even more ubiquitous. Prominent groundwater pathways have also been identified in fractured crystalline bedrock (Andermann et al.2012; Pohl et al.2015; Tague and Grant2009), morainic deposits (Favier et al.2008; Minaya2016; Somers et al.2016), and alpine meadow soils (Gordon et al.2015; Loheide et al.2009; Lowry et al.2010). Even in settings with limited groundwater networks, talus slopes and rock glaciers can serve as localized areas of meltwater recharge (Baraer et al.2015; Clow et al.2003; Harrington et al.2018). More arid settings than the subhumid Gavilan Machay could have a higher proportion of glacier melt in groundwater due to less precipitation inputs, but our results also indicate that precipitation may serve to enhance meltwater contributions.

The multi-scale temporal variability of relative melt contributions to discharge has important implications for how to determine the hydrologic role of glaciers in watersheds, as well as for water resource management in fast-changing glacierized systems. Care must be taken in the implementation and interpretation of commonly employed tracer analyses. Potentially large diurnal fluctuations make it imperative to collect samples over consistent times of day, and weekly to interannual variability complicates extrapolations from single synoptic sampling estimates. Recharge by glacier melt further confounds the interpretation of groundwater as a source entirely distinct from surficial meltwater. These uncertainties, along with additional errors caused by heterogeneous groundwater chemistry and the choice of sample locations, limit the ability of tracer-based analyses to constrain dynamic melt contributions to discharge. Model simulations provide ideal temporal resolution, but they suffer from their own disadvantages, such as intensive data input requirements and uncertainties in parameter calibration.

For water resources, weekly to multi-year melt contributions to discharge are of greater interest than hourly fluctuations. At those timescales of concern, rain events and wet periods can accentuate relative melt contributions to streamflow in humid glacierized systems. This signifies a bonus in water yield, but it also intensifies discharge variability over weekly and seasonal timescales, which can pose challenges if water storage infrastructure is unavailable. On an interannual basis, melt can augment discharge during warm and dry El Niño events, but glacierized watersheds will likely experience the greatest melt contributions during wetter times, when they are under the combined effects of enhanced ablation (due to transfer of heat from rain to ice) and antecedent moist conditions (due to precipitation). Looking at downstream implications, the small (7.5 km2) Gavilan Machay headwater catchment contributes a range of 9 %–26 % of the discharge to the Boca Toma diversion point, which corresponds to surficial meltwater making up a range of 4 %–15 % of the water to the irrigation system, based on mixing model estimates. Although this amount may seem small, La Frenierre (2014) showed that farming communities cannot afford to lose any of the water; already, the irrigation system consistently fails to deliver its current allocations. Furthermore, if groundwater at Gavilan Machay also contains meltwater, as our simulations suggest, the actual total amount of meltwater contribution could be even higher than that estimated for surficial runoff of meltwater by the mixing model. Additional downstream monitoring would enable further extensions of the watershed model to investigate how meltwater contributions and temporal discharge variabilities found in Gavilan Machay propagate downgradient to successively larger watersheds, as non-glacier-sourced groundwater contributions increase. Spatial patterns of groundwater contributions within Gavilan Machay reveal sharp increases where geologic features likely create localized discharge points (Fig. S12). This indicates that extrapolations downstream will likely depend on geological conditions that control groundwater, in addition to watershed size and climate inputs.

In the future, should Reschreiter Glacier disappear completely, overall discharge in Gavilan Machay could decrease by up to about 50 %, even if precipitation and temperature remain the same and relatively constant groundwater flow continues. The exact decrease will depend on how much of our melt estimate may include contributions by snowmelt and melt from freshly accumulated ice that will persist as discharge sources. Overall, our findings suggest that, in response to glacier loss in a warming climate, glacierized watersheds in the humid inner tropics may eventually experience steadier discharge, but potentially at significantly decreased rates.

Code and data availability

The Flux-PIHM version 0.5.0 (Shi et al.2013) used for this paper is available on GitHub at (last access: 2 January 2018). The PIHMgis software (Bhatt et al.2014) used to create input files is available at (last access: 2 January 2018). Hydrochemical and hydroclimatic data are available in the Supplement. The corresponding author maybe contacted to obtain model outputs.


The supplement related to this article is available online at:

Author contributions

LS implemented the time series analysis and hydrologic modeling, assisted with the mixing model analysis, and contributed to writing the manuscript. RTM executed the field work in 2016–2017, implemented the mixing model analysis, and contributed to writing the manuscript. GHCN conceived of the project together with JLF, oversaw the analysis, and contributed to writing the manuscript. JLF also established the study site and led the fieldwork. ADW provided field instrumentation for 2015–2017 and contributed to fieldwork. MB provided the HBCM mixing model software and guided its implementation. WZ and LL assisted with the Flux-PIHM hydrologic model implementation. BGM supported establishment of the field site and supervised the initial work.

Competing interests

The authors declare that they have no conflict of interest.


Funding from NSF (EAR-1759071) helped support this work. Rachel T. McLaughlin received travel support from the University of Minnesota's Walter H. Judd Fellowship. G.-H. Crystal Ng and Andrew D. Wickert received financial support from the University of Minnesota. Jeff La Frenierre received financial support from the National Science Foundation's Doctoral Dissertation Research Improvement Grant (1103235), the Fulbright Commission of Ecuador, the Geological Society of America, and Gustavus Adolphus College. The authors would like to acknowledge field assistance from Chad Sandell, Casey Decker, Helen Thompson, and Abigail Michels; lab support from Jeff Jeremiason, Chris Harmes, and Scott Alexander; and land-cover mapping assistance from Josh Zoellmer.

Edited by: Wouter Buytaert
Reviewed by: two anonymous referees


Andermann, C., Longuevergne, L., Bonnet, S., Crave, A., Davy, P., and Gloaguen, R.: Impact of transient groundwater storage on the discharge of Himalayan rivers, Nat. Geosci., 5, 127–132,, 2012. a, b, c

Andrews, E. D.: Measurement and computation of bed-material discharge in a shallow sand-bed stream, Muddy Creek, Wyoming, Water Resour. Res., 17, 131–141,, 1981. a

Bao, C., Li, L., Shi, Y., and Duffy, C.: Understanding watershed hydrogeochemistry: 1. Development of RT-Flux-PIHM, Water Resour. Res., 53, 2328–2345, 2017. a

Baraer, M., McKenzie, J. M., Mark, B. G., Bury, J., and Knox, S.: Characterizing contributions of glacier melt and groundwater during the dry season in a poorly gauged catchment of the Cordillera Blanca (Peru), Adv. Geosci., 22, 41–49,, 2009. a, b, c, d, e, f

Baraer, M., Mark, B. G., McKenzie, J. M., Condom, T., Bury, J., Huh, K.-I., Portocarrero, C., Gómez, J., and Rathay, S.: Glacier recession and water resources in Peru's Cordillera Blanca, J. Glaciol., 58, 134–150,, 2012. a, b

Baraer, M., McKenzie, J., Mark, B. G., Gordon, R., Bury, J., Condom, T., Gomez, J., Knox, S., and Fortner, S. K.: Contribution of groundwater to the outflow from ungauged glacierized catchments: a multi-site study in the tropical Cordillera Blanca, Peru, Hydrol. Proc., 29, 2561–2581,, 2015. a, b, c, d, e, f, g, h, i

Barba, D., Samaniego, P., Eissen, J.-P., Robin, C., Fornari, M., Cotten, J., and Beate, B.: Geology and structure of the late Pleistocene to Holocene Chimborazo stratovolcano (Ecuador), 6th International Symposium on Andean Geodynamics (ISAG 2005), 12–14 September 2005, Barcelona, Spain, 90–93, 2005. a

Barnett, T. P., Adam, J. C., and Lettenmaier, D. P.: Potential impacts of a warming climate on water availability in snow-dominated regions, Nature, 438, 303–309,, 2005. a

Bhatt, G., Kumar, M., and Duffy, C. J.: A tightly coupled GIS and distributed hydrologic modeling framework, Environ. Modell. Softw., 62, 70–84,, 2014. a, b

Bradley, R. S.: CLIMATE CHANGE: Threats to Water Supplies in the Tropical Andes, Science, 312, 1755–1756,, 2006. a, b

Bradley, R. S., Vuille, M., Hardy, D., and Thompson, L. G.: Low latitude ice cores record Pacific sea surface temperatures, Geophys. Res. Lett., 30, 2–5,, 2003. a

Braun, L. N., Weber, M., and Schulz, M.: Consequences of climate change for runoff from Alpine regions, Ann. Glaciol., 31, 19–25,, 2000. a

Bury, J., Mark, B. G., Carey, M., Young, K. R., McKenzie, J. M., Baraer, M., French, A., and Polk, M. H.: New Geographies of Water and Climate Change in Peru: Coupled Natural and Social Transformations in the Santa River Watershed, Ann. Assoc. Am. Geogr., 103, 363–374,, 2013. a

Bury, J. T., Mark, B. G., McKenzie, J. M., French, A., Baraer, M., Huh, K. I., Zapata Luyo, M. A., and Gómez López, R. J.: Glacier recession and human vulnerability in the Yanamarey watershed of the Cordillera Blanca, Peru, Clim. Change, 105, 179–206,, 2011. a

Buytaert, W. and Beven, K.: Models as multiple working hypotheses: hydrological simulation of tropical alpine wetlands, Hydrol. Proc., 25, 1784–1799,, 2011. a, b, c

Buytaert, W., Célleri, R., De Bièvre, B., Cisneros, F., Wyseure, G., Deckers, J., and Hofstede, R.: Human impact on the hydrology of the Andean páramos, Earth-Sci. Rev., 79, 53–72,, 2006. a, b, c, d

Caceres, B., Francou, B., Favier, V., Bontron, G., Tachker, P., Bucher, R., Taupin, J.-D., Vuille, M., Maisincho, L., and Delachaux, F: Glacier 15, Antisana, Ecuador: its glaciology and relations to water resources, in: Climate Variability and Change – Hydrological Impacts, 308, IAHS, edited by: Siegfried Demuth, Alan Gustard, Eduardo Planos, Fred Scatena & Eric Servat., IAHS Publication, 479–482, 2006. a

Carey, M., Molden, O. C., Rasmussen, M. B., Jackson, M., Nolin, A. W., and Mark, B. G.: Impacts of Glacier Recession and Declining Meltwater on Mountain Societies, Ann. Am. Assoc. Geogr., 107, 350–359,, 2017. a, b

Cauvy-Fraunié, S., Condom, T., Rabatel, A., Villacis, M., Jacobsen, D., and Dangles, O.: Technical Note: Glacial influence in tropical mountain hydrosystems evidenced by the diurnal cycle in water levels, Hydrol. Earth Syst. Sci., 17, 4803–4816,, 2013. a

Chen, J., Knight, R., and Zebker, H. A.: The Temporal and Spatial Variability of the Confined Aquifer Head and Storage Properties in the San Luis Valley, Colorado Inferred From Multiple InSAR Missions, Water Resour. Res., 53, 9708–9720,, 2017. a

Christophersen, N., Neal, C., Hooper, R. P., Vogt, R. D., and Andersen, S.: Modelling streamwater chemistry as a mixture of soilwater end-members-a step towards second-generation acidification models, J. Hydrol., 116, 307–320, 1990. a, b

Clapperton, C. M.: Glacial and volcanic geomorphology of the Chimborazo-Carihuairazo Massif, Ecuadorian Andes, T. Roy. Soc. London, 81, 91–116,, 1990. a, b

Clow, D., Schrott, L., Webb, R., Campbell, D., Torizzo, A., and Dornblaser, M.: Ground Water Occurrence and Contributions to Streamflow in an Alpine Catchment, Colorado Front Range, Ground Water, 41, 937–950,, 2003. a, b, c

Crossman, J., Bradley, C., Boomer, I., and Milner, A. M.: Water Flow Dynamics of Groundwater-Fed Streams and Their Ecological Significance in a Glacierized Catchment, Arct. Antarct. Alp. Res., 43, 364–379,, 2011. a

Duffy, C., Shi, Y., Davis, K., Slingerland, R., Li, L., Sullivan, P. L., Goddéris, Y., and Brantley, S. L.: Designing a Suite of Models to Explore Critical Zone Function, Proced. Earth Plan. Sci., 10, 7–15,, 2014. a

Ek, M. B., Mitchell, K. E., Lin, Y., Rogers, E., Grunmann, P., Koren, V., Gayno, G., and Tarpley, J. D.: Implementation of Noah land surface model advances in the National Centers for Environmental Prediction operational mesoscale Eta model, J. Geophys. Res.-Atmos, 108,, 2003. a

Engel, M., Penna, D., Bertoldi, G., Dell'Agnese, A., Soulsby, C., and Comiti, F.: Identifying run-off contributions during melt-induced run-off events in a glacierized alpine catchment, Hydrol. Proc., 30, 343–364,, 2016. a

Escher-Vetter, H., Reinwarth, O., and Rentsch, H.: Two decades of runoff measurements (1974 to 1993) at the pegelstation Vernagtbach/Oetztal alps, Zeitschrift für Gletscherkunde und Glazialgeologie, 30, 99–107, 1994. a

Farr, T. G., Rosen, P. A., Caro, E., Crippen, R., Duren, R., Hensley, S., Kobrick, M., Paller, M., Rodriguez, E., and Roth, L.: The shuttle radar topography mission, Rev. Geophys., 45, 8755–1209,, 2007. a

Favier, V.: Glaciers of the outer and inner tropics: A different behaviour but a common response to climatic forcing, Geophys. Res. Lett., 31, L16403,, 2004. a, b, c

Favier, V., Coudrain, A., Cadier, E., Francou, B., Ayabaca, E., Maisincho, L., Praderio, E., Villacis, M., and Wagnon, P.: Evidence of groundwater flow on Antizana ice-covered volcano, Ecuador/Mise en évidence d'écoulements souterrains sur le volcan englacé Antizana, Equateur, Hydrol. Sci. J., 53, 278–291,, 2008. a, b, c, d

Fernández, A. and Mark, B. G.: Modeling modern glacier response to climate changes along the Andes Cordillera: A multiscale review, J. Adv. Model. Earth Syst., 8, 467–495, 2016. a, b, c

Fountain, A. G. and Tangborn, W. V.: The Effect of Glaciers on Streamflow Variations, Water Resour. Res., 21, 579–586,, 1985. a

Francou, B.: Tropical climate change recorded by a glacier in the central Andes during the last decades of the twentieth century: Chacaltaya, Bolivia, 16S, J. Geophys. Res., 108, 4154,, 2003. a

Francou, B.: New evidence for an ENSO impact on low-latitude glaciers: Antizana 15, Andes of Ecuador, J. Geophys. Res., 109, D18106,, 2004. a, b

Frenierre, J. L. and Mark, B. G.: A review of methods for estimating the contribution of glacial meltwater to total watershed discharge, Prog. Phys. Geogr., 38, 173–200,, 2014. a

Frisbee, M. D., Phillips, F. M., Campbell, A. R., Liu, F., and Sanchez, S. A.: Streamflow generation in a large, alpine watershed in the southern Rocky Mountains of Colorado: Is streamflow generation simply the aggregation of hillslope runoff responses?, Water Resour. Res., 47, 1–18,, 2011. a

Gabbi, J., Carenzo, M., Pellicciotti, F., Bauder, A., and Funk, M.: A comparison of empirical and physically based glacier surface melt models for long-term simulations of glacier response, J. Glaciol., 60, 1140–1154,, 2014. a

Gordon, R. P., Lautz, L. K., McKenzie, J. M., Mark, B. G., Chavez, D., and Baraer, M.: Sources and pathways of stream generation in tropical proglacial valleys of the Cordillera Blanca, Peru, J. Hydrol., 522, 628–644,, 2015. a, b

Harden, D. R.: California geology, Prentice Hall, 2004. a

Harpold, A. A. and Brooks, P. D.: Humidity determines snowpack ablation under a warming climate, P. Natl. Acad. Sci. USA, 115, 1215–1220,, 2018. a, b

Harrington, J. S., Mozil, A., Hayashi, M., and Bentley, L. R.: Groundwater flow and storage processes in an inactive rock glacier, Hydrol. Proc., 32, 3070–3088,, 2018. a, b, c, d

He, Z., Vorogushyn, S., Unger-Shayesteh, K., Gafurov, A., Kalashnikova, O., Omorova, E., and Merz, B.: The value of hydrograph partitioning curves for calibrating hydrological models in glacierized basins, Water Resour. Res., 54, 1–52,, 2018. a

Hock, R.: A distributed temperature-index ice- and snowrnelt model including potential direct solar radiation, J. Glaciol., 45, 101–111,, 1999. a

Hock, R.: Glacier melt: a review of processes and their modelling, Prog. Phys. Geogr., 29, 362–391, 2005. a, b, c

Hood, J. L., Roy, J. W., and Hayashi, M.: Importance of groundwater in the water balance of an alpine headwater lake, Geophys. Res. Lett., 33, L13405,, 2006. a

Hooper, R. P.: Diagnostic tools for mixing models of stream water chemistry, Water Resour. Res., 39, 1055,, 2003. a

Hooper, R. P. and Shoemaker, C. A.: A Comparison of Chemical and Isotopic Hydrograph Separation, Water Resour. Res., 22, 1444–1454,, 1986. a

Huss, M., Bauder, A., Funk, M., and Hock, R.: Determination of the seasonal mass balance of four Alpine glaciers since 1865, J. Geophys. Res., 113, F01015,, 2008. a

Huss, M., Funk, M., and Ohmura, A.: Strong Alpine glacier melt in the 1940s due to enhanced solar radiation, Geophys. Res. Lett., 36, L23501,, 2009. a, b

Huth, A. K., Leydecker, A., Sickman, J. O., and Bales, R. C.: A two-component hydrograph separation for three high-elevation catchments in the Sierra Nevada, California, Hydrol. Proc., 18, 1721–1733,, 2004. a

Immerzeel, W. W., van Beek, L. P. H., and Bierkens, M. F. P.: Climate Change Will Affect the Asian Water Towers, Science, 328, 1382–1385,, 2010. a

INEC: Cencos de Poblacion y Vivienda 2010, Quito, Ecuador, 2010. a

James, A. L. and Roulet, N. T.: Investigating the applicability of end-member mixing analysis (EMMA) across scale: A study of eight small, nested catchments in a temperate forested watershed, Water Resour. Res., 42, W08434,, 2006. a

Jansson, P., Hock, R., and Schneider, T.: The concept of glacier storage: a review, J. Hydrol., 282, 116–129,, 2003. a

Jordan, E., González, J., Castillo, K., Torres, J., Ungerechts, L., Vélez, F., Blanco, D., and Cruz, M.: Ortofotomapa del Chimborazo y su valor como diagnóstico para cambios climaticos en relacion con otros glaciares tropicales [Orthophotomap of Chimborazo and its value as a diagnosis for climatic changes in relation to other tropical glaciers], Glaciares, Nieves y Hielos de America Latina: Cambio Climatico y Amenazas, 239–260, 2010. a

Juen, I., Kaser, G., and Georges, C.: Modelling observed and future runoff from a glacierized tropical catchment (Cordillera Blanca, Perú), Glob. Planet. Change, 59, 37–48,, 2007. a, b

Kaser, G.: A review of the modern fluctuations of tropical glaciers, Glob. Planet. Change, 22, 93–103,, 1999. a

Kaser, G. and Osmaston, H.: Tropical glaciers, Cambridge University Press, 2002. a, b

Kendall, C., Doctor, D., Drever, J., Holland, H., and Turekian, K.: Stable isotope applications in hydrological studies, Treatise on geochemistry, 5, p. 605,, 2003. a

La Frenierre, J.: Rapid downward deflation of a tropical-debris covered glacier: an analysis from Volcán Chimborazo, Ecuador, in: AGU Fall Meeting Abstracts, 2016. a

La Frenierre, J. and Mark, B. G.: Detecting Patterns of Climate Change at Volcán Chimborazo, Ecuador, by Integrating Instrumental Data, Public Observations, and Glacier Change Analysis, Ann. Am. Assoc. Geogr., 107, 979–997,, 2017. a, b, c, d, e

La Frenierre, J. D.: Assessing the Hydrologic Implications of Glacier Recession and the Potential for Water Resources Vulnerability at Volcán Chimborazo, Ecuador, PhD diss., The Ohio State University, Prog. Phys. Geog., 255, 2014. a

Lang, H.: Forecasting Meltwater Runoff from Snow-Covered Areas and from Glacier Basins, in: River flow modelling and forecasting, edited by: Kraijenhoff, D. A. and Moll, J. R., Springer Netherlands, Dordrecht, 99–127,, 1986. a, b

Li, L., Bao, C., Sullivan, P. L., Brantley, S., Shi, Y., and Duffy, C.: Understanding watershed hydrogeochemistry: 2. Synchronized hydrological and geochemical processes drive stream chemostatic behavior, Water Resour. Res., 53, 2346–2367, 2017. a

Liu, F., Williams, M. W., and Caine, N.: Source waters and flow paths in an alpine catchment, Colorado Front Range, United States, Water Resour. Res., 40, 1–16,, 2004. a, b

Loheide, S. P., Deitchman, R. S., Cooper, D. J., Wolf, E. C., Hammersmark, C. T., and Lundquist, J. D.: A framework for understanding the hydroecology of impacted wet meadows in the Sierra Nevada and Cascade Ranges, California, USA, Hydrogeol. J., 17, 229–246,, 2009. a

Lowry, C. S., Deems, J. S., Loheide II, S. P., and Lundquist, J. D.: Linking snowmelt-derived fluxes and groundwater flow in a high elevation meadow system, Sierra Nevada Mountains, California, Hydrol. Proc., 24, 2821–2833,, 2010. a, b, c

Lowry, C. S., Loheide, S. P., Moore, C. E., and Lundquist, J. D.: Groundwater controls on vegetation composition and patterning in mountain meadows, Water Resour. Res., 47, 1–16,, 2011. a

Luce, C. H.: Effects of Climate Change on Snowpack, Glaciers, and Water Resources in the Northern Rockies, in: Climate Change and Rocky Mountain Ecosystems, edited by: Halofsky, J. E. and Peterson, D. L., Springer International Publishing, Cham., 25–36,, 2018. a

Mackay, A.: Climate Change 2007: Impacts, Adaptation and Vulnerability. Contribution of Working Group II to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change, J. Environ. Qual., 37, 2407,, 2008. a

Manciati, C., Villacís, M., Taupin, J.-D., Cadier, E., Galárraga-Sánchez, R., and Cáceres, B.: Empirical mass balance modelling of South American tropical glaciers: case study of Antisana volcano, Ecuador, Hydrol. Sci. J., 59, 1519–1535,, 2014. a

Mark, B. G.: Tracing tropical Andean glaciers over space and time: Some lessons and transdisciplinary implications, Glob. Planet. Change, 60, 101–114,, 2008. a

Mark, B. G. and Mckenzie, J. M.: Tracing Increasing Tropical Andean Glacier Melt with Stable Isotopes in Water, Environ. Sci. Tech., 41, 6955–6960,, 2007. a, b, c, d

Mark, B. G. and Seltzer, G. O.: Tropical glacier meltwater contribution to stream discharge: a case study in the Cordillera Blanca, Peru, J. Glaciol., 49, 271–281,, 2003. a, b, c, d, e

Mark, B. G., McKenzie, J. M., and Gómez, J.: Hydrochemical evaluation of changing glacier meltwater contribution to stream discharge: Callejon de Huaylas, Peru / Evaluation hydrochimique de la contribution évolutive de la fonte glaciaire à l'écoulement fluvial: Callejon de Huaylas, Pérou, Hydrol. Sci. J., 50, 975–988,, 2005. a

Markovich, K. H., Maxwell, R. M., and Fogg, G. E.: Hydrogeological response to climate change in alpine hillslopes, Hydrol. Proc., 30, 3126–3138,, 2016. a, b, c

Maussion, F., Gurgiser, W., Großhauser, M., Kaser, G., and Marzeion, B.: ENSO influence on surface energy and mass balance at Shallap Glacier, Cordillera Blanca, Peru, The Cryosphere, 9, 1663–1683,, 2015. a

McLaughlin, R.: Hydrochemical Signatures of Glacial Meltwater on Volcán Chimborazo, Ecuador, M. Sc. dissertation, University of Minnesota, 2017. a, b, c

Messerli, B., Viviroli, D., and Weingartner, R.: Mountains of the World: Vulnerable Water Towers for the 21st Century, AMBIO Special Report, 13, 29–34,, 2004. a

Minaya, Maldonado, V. G.: Ecohydrology of the Andes Páramo Region, PhD diss, IHE Delft Institute for Water Education, 2016. a, b, c, d, e, f, g, h, i

Ng, G.-H. C., Wickert, A. D., Somers, L. D., Saberi, L., Cronkite-Ratcliff, C., Niswonger, R. G., and McKenzie, J. M.: GSFLOW-GRASS v1.0.0: GIS-enabled hydrologic modeling of coupled groundwater-surface-water systems, Geosci. Model Dev., 11, 4755–4777,, 2018. a

NRCS, U.: Part 630 Hydrology National Engineering Handbook, USDA-NRCS (United States Department of Agriculture- Natural Resources Conservation Service), 2009. a

Omani, N., Srinivasan, R., Smith, P. K., and Karthikeyan, R.: Glacier mass balance simulation using SWAT distributed snow algorithm, Hydrol. Sci. J., 62, 546–560,, 2017. a, b

Ostheimer, G. J., Hadjivasiliou, H., Kloer, D. P., Barkan, A., and Matthews, B. W.: Structural Analysis of the Group II Intron Splicing Factor CRS2 Yields Insights into its Protein and RNA Interaction Surfaces, J. Molec. Biol., 345, 51–68,, 2005. a

Pellicciotti, F., Brock, B., Strasser, U., Burlando, P., Funk, M., and Corripio, J.: An enhanced temperature-index glacier melt model including the shortwave radiation balance: development and testing for Haut Glacier d'Arolla, Switzerland, J. Glaciol., 51, 573–587, 2005. a

Pepin, N., Bradley, R. S., Diaz, H. F., Baraer, M., Caceres, E. B., Forsythe, N., Fowler, H., Greenwood, G., Hashmi, M. Z., Liu, X. D., Miller, J. R., Ning, L., Ohmura, A., Palazzi, E., Rangwala, I., Schöner, W., Severskiy, I., Shahgedanova, M., Wang, M. B., Williamson, S. N., and Yang, D. Q.: Elevation-dependent warming in mountain regions of the world, Nat. Clim. Change, 5, 424–430,, 2015. a

Podwojewski, P., Poulenard, J., Zambrana, T., and Hofstede, R.: Overgrazing effects on vegetation cover and properties of volcanic ash soil in the pàramo of Llangahua and al Esperanza (Tungurahua, Ecuador), Management, 18, 45–55,, 2002. a, b, c, d

Pohl, E., Knoche, M., Gloaguen, R., Andermann, C., and Krause, P.: Sensitivity analysis and implications for surface processes from a hydrological modelling approach in the Gunt catchment, high Pamir Mountains, Earth Surf. Dynam., 3, 333–362,, 2015. a, b

Polk, M. H., Young, K. R., Baraer, M., Mark, B. G., McKenzie, J. M., Bury, J., and Carey, M.: Exploring hydrologic connections between tropical mountain wetlands and glacier recession in Peru's Cordillera Blanca, Appl. Geogr., 78, 94–103,, 2017. a

Pribulick, C. E., Foster, L. M., Bearup, L. A., Navarre-Sitchler, A. K., Williams, K. H., Carroll, R. W., and Maxwell, R. M.: Contrasting the hydrologic response due to land cover and climate change in a mountain headwaters system, Ecohydrology, 9, 1431–1438,, 2016. a

Qu, Y. and Duffy, C. J.: A semidiscrete finite volume formulation for multiprocess watershed simulation, Water Resour. Res., 43, 1–18,, 2007. a

Rabatel, A., Francou, B., Soruco, A., Gomez, J., Cáceres, B., Ceballos, J. L., Basantes, R., Vuille, M., Sicart, J.-E., Huggel, C., Scheel, M., Lejeune, Y., Arnaud, Y., Collet, M., Condom, T., Consoli, G., Favier, V., Jomelli, V., Galarraga, R., Ginot, P., Maisincho, L., Mendoza, J., Ménégoz, M., Ramirez, E., Ribstein, P., Suarez, W., Villacis, M., and Wagnon, P.: Current state of glaciers in the tropical Andes: a multi-century perspective on glacier evolution and climate change, The Cryosphere, 7, 81–102,, 2013. a

Reveillet, M., Vincent, C., Six, D., and Rabatel, A.: Which empirical model is best suited to simulate glacier mass balances?, J. Glaciol., 63, 39–54,, 2017. a

Rodell, M., Houser, P. R., Jambor, U. E. A., Gottschalck, J., Mitchell, K., Meng, C. J., Arsenault, K., Cosgrove, B., Radakovich, J., and Bosilovich, M.: The global land data assimilation system, B. Am. Meteorol. Soc., 85, 381–394, 2004. a

Ryu, J.-S., Lee, K.-S., and Chang, H.-W.: Hydrochemistry and isotope geochemistry of Song Stream, a headwater tributary of the South Han River, South Korea, Geosci. J., 11, 157–164,, 2007. a

Samaniego, P., Barba, D., Robin, C., Fornari, M., and Bernard, B.: Eruptive history of Chimborazo volcano (Ecuador): A large, ice-capped and hazardous compound volcano in the Northern Andes, J. Volcanol. Geoth. Res., 221–222, 33–51,, 2012. a

Sarmiento, L.: Water Balance and Soil Loss under Long Fallow Agriculture in the Venezuelan Andes, Mt. Res. Dev., 20, 246–253, 2000. a

Schmieder, J., Garvelmann, J., Marke, T., and Strasser, U.: Spatio-temporal tracer variability in the glacier melt end-member – How does it affect hydrograph separation results?, Hydrol. Proc., 32, 1828–1843,, 2018. a

Shi, Y., Davis, K. J., Duffy, C. J., and Yu, X.: Development of a Coupled Land Surface Hydrologic Model and Evaluation at a Critical Zone Observatory, J. Hydrometeorol., 14, 1401–1420,, 2013. a, b, c

Sicart, J. E., Hock, R., and Six, D.: Glacier melt, air temperature, and energy balance in different climates: The Bolivian Tropics, the French Alps, and northern Sweden, J. Geophys. Res., 113, D24113,, 2008. a

Smith, J. A., Mark, B. G., and Rodbell, D. T.: The timing and magnitude of mountain glaciation in the tropical Andes, J. Quat. Sci., 23, 609–634,, 2008. a, b

Somers, L. D., Gordon, R. P., McKenzie, J. M., Lautz, L. K., Wigmore, O., Glose, A. M., Glas, R., Aubry-Wake, C., Mark, B., Baraer, M., and Condom, T.: Quantifying groundwater-surface water interactions in a proglacial valley, Cordillera Blanca, Peru, Hydrol. Proc., 30, 2915–2929,, 2016. a, b

Soruco, A., Vincent, C., Rabatel, A., Francou, B., Thibert, E., Sicart, J. E., and Condom, T.: Contribution of glacier runoff to water resources of La Paz city, Bolivia (16 S), Ann. Glaciol., 56, 147–154,, 2015. a, b

Soulsby, C., Petry, J., Brewer, M., Dunn, S., Ott, B., and Malcolm, I.: Identifying and assessing uncertainty in hydrological pathways: a novel approach to end member mixing in a Scottish agricultural catchment, J. Hydrol., 274, 109–128,, 2003. a

Suecker, J. K., Ryan, J. N., Kendall, C., and Jarrett, R. D.: Determination of hydrologic pathways during snowmelt for alpine/subalpine basins, Rocky Mountain National Park, Colorado, Water Resour. Res., 36, 63–75,, 2000. a

Tague, C. and Grant, G. E.: Groundwater dynamics mediate low-flow response to global warming in snow-dominated alpine regions, Water Resour. Res., 45, 1–12,, 2009. a, b, c, d

Tague, C., Grant, G., Farrell, M., Choate, J., and Jefferson, A.: Deep groundwater mediates streamflow response to climate warming in the Oregon Cascades, Clim. Change, 86, 189–210,, 2008. a, b, c

Veettil, B. K., Leandro Bayer Maier, E., Bremer, U. F., and de Souza, S. F.: Combined influence of PDO and ENSO on northern Andean glaciers: a case study on the Cotopaxi ice-covered volcano, Ecuador, Clim. Dynam., 43, 3439–3448,, 2014a. a

Veettil, B. K., Maier, E. L. B., Bremer, U. F., and de Souza, S. F.: Combined influence of PDO and ENSO on northern Andean glaciers: a case study on the Cotopaxi ice-covered volcano, Ecuador, Clim. Dynam., 43, 3439–3448, 2014b.  a

Veettil, B. K., Bremer, U. F., de Souza, S. F., Maier, E. L. B., and Simões, J. C.: Influence of ENSO and PDO on mountain glaciers in the outer tropics: case studies in Bolivia, Theor. Appl. Climatol., 125, 757–768,, 2016. a

Vermote, E.: MOD09A1 MODIS/Terra Surface Reflectance 8-Day L3 Global 500 m SIN Grid V006, NASA EOSDIS Land Processes DAAC, 2015. a

Vuille, M. and Bradley, R. S.: The tropical Andes, Geophys. Res. Lett., 27, 3885–3888, 2000. a

Vuille, M. and Keimig, F.: Interannual Variability of Summertime Convective Cloudiness and Precipitation in the Central Andes Derived from ISCCP-B3 Data, J. Climate, 17, 3334–3348,<3334:IVOSCC>2.0.CO;2, 2004. a, b

Vuille, M., Bradley, R. S., Werner, M., and Keimig, F.: 20th Century Climate Change in the Tropical Andes: Observations and Model Results, in: Climate variability and change in high elevation regions: Past, Present & Future, Springer, 59, 75–99, 2003. a

Vuille, M., Francou, B., Wagnon, P., Juen, I., Kaser, G., Mark, B. G., and Bradley, R. S.: Climate change and tropical Andean glaciers: Past, present and future, Earth-Sci. Rev., 89, 79–96,, 2008. a

Vuille, M., Carey, M., Huggel, C., Buytaert, W., Rabatel, A., Jacobsen, D., Soruco, A., Villacis, M., Yarleque, C., Elison Timm, O., Condom, T., Salzmann, N., and Sicart, J.-E.: Rapid decline of snow and ice in the tropical Andes – Impacts, uncertainties and challenges ahead, Earth-Sci. Rev., 176, 195–213,, 2018. a

Wagnon, P., Ribstein, P., Francou, B., and Sicart, J. E.: Anomalous heat and mass budget of Glaciar Zongo, Bolivia, during the 1997/98 El Niño year, J. Glaciol., 47, 21–28, 2001. a, b

Wang, X., Sun, L., Zhang, Y., and Luo, Y.: Rationalization of Altitudinal Precipitation Profiles in a Data-Scarce Glacierized Watershed Simulation in the Karakoram, Water, 8, 186,, 2016. a

Wickert, A. D.: The ALog: Inexpensive, Open-Source, Automated Data Collection in the Field, The Bulletin of the Ecological Society of America, 95, 166–176, 2014. a

Wickert, A. D., Sandell, C. T., Schulz, B., and Ng, G.-H. C.: Open-source Arduino-derived data loggers designed for field research, Hydrol. Earth Syst. Sci. Discuss.,, in review, 2018. a

Wilson, A. M., Williams, M. W., Kayastha, R. B., and Racoviteanu, A.: Use of a hydrologic mixing model to examine the roles of meltwater, precipitation and groundwater in the Langtang River basin, Nepal, Ann. Glaciol., 57, 155–168,, 2016. a

Wosten, J. H. M., Lilly, A., Nemes, A., and Le Bas, C.: Development and use of a database of hydraulic properties of European soils, Geoderma, 90, 169–185, 1999. a

Zhang, Y., Slingerland, R., and Duffy, C.: Fully-coupled hydrologic processes for modeling landscape evolution, Environ. Modell. Softw., 82, 89–107, 2016. a


The requested paper has a corresponding corrigendum published. Please read the corrigendum first before downloading the article.

Short summary
The relationship among glacier melt, groundwater, and streamflow remains highly uncertain, especially in tropical glacierized watersheds in response to climate. We implemented a multi-method approach and found that melt contribution varies considerably and may drive streamflow variability at hourly to multi-year timescales, rather than buffer it, as commonly thought. Some of the melt contribution occurs through groundwater pathways, resulting in longer timescale interactions with streamflow.