glacial influence in tropical mountain hydrosystems evidenced by the diurnal cycle in water levels

.


Introduction
In view of the accelerated glacier shrinking worldwide (Lemke et al., 2007;Rabatel et al., 2013;Sakakibara et al., 2013), coupling glacier and glacier-fed hydrosystem dynamics is a timely research thematic (Bradley et al., 2006; Published by Copernicus Publications on behalf of the European Geosciences Union.S. Cauvy-Fraunié et al.: Technical Note: Glacial influence evidenced by the diurnal cycle in water levels Jacobsen et al., 2012).While at the early stages of glacier retreat the reduction in ice volume would yield a significant increase in annual runoff (see the conceptual model presented by Baraer et al., 2012), after a critical threshold (depending on the glacier size) the annual discharge would decrease up to the end of the glacial influence on outflow (Huss et al., 2008).Worldwide, glacial river discharges have shown both increasing and decreasing trends, depending on ice cover in the catchment, the study region, and where the glacier stands along the deglaciation trajectory (Fleming and Clarke, 2003;Stahl and Moore, 2006;Casassa et al., 2009;Moore et al., 2009;Dahlke et al., 2012;Fleming and Weber, 2012).In this context, a growing number of studies have quantitatively explored the potential future impacts of various climate change and glacier recession scenarios upon water resources, using modern glaciological and hydrological modeling techniques (e.g., Schaefli et al., 2007;Huss et al., 2008;Stahl et al., 2008;Villacis, 2008;Jost et al., 2012;Yao et al., 2013).Such studies have demonstrated that glacier change effects are likely to be hydrologically substantial, even in relatively lightly glacierized basins.Modifications in water regimes may have significant consequences on water quality, aquatic biota and water security for human populations (Barnett et al., 2005;Brown et al., 2010;Kaser et al., 2010).
In this context, detecting the influence of glacial runoff to stream discharge has become a key challenge for a broad community of researchers, including glaciologists, hydrologists, water managers and ecologists (Brown et al., 2010;Baraer et al., 2012;Cauvy-Fraunié et al., 2013).To measure glacial meltwater influence on mountain streams, most approaches have focused on determining the water source contribution in glacierized river basins (e.g., glacial melt, snowmelt, rain, and groundwater) using methods ranging from thermal and discharge balances to stable isotope analyses (Huss et al., 2008;Kaser et al., 2010;Dahlke et al., 2012) or hydro-glaciological model (Condom et al., 2012)."Glacier indices" have also been developed such as (1) the glacial index (Jacobsen and Dangles, 2012) calculated from glacier size and distance from the glacier terminus, (2) the percentage of glacier cover in the catchment (Rott et al., 2006;Füreder, 2007;Milner et al., 2009), (3) the Alpine River and Stream Ecosystem classification (ARISE, Brown et al., 2009) based on hydrochemical analyses of water samples and statistical mixing models, and (4) the "glaciality index" (Ilg and Castella, 2006) based on four physico-chemical habitat variables (water temperature, channel stability, conductivity, and suspended sediment concentration).
However, a major challenge for these methodologies is the need to incorporate the high spatiotemporal variability of the different water source contributions in glacierized catchments (Brown et al., 2009).In this respect, existing glacier indices suffer from several limitations.First, although commonly used, the estimation of glacier cover in the catchment may be neither an easy nor a reliable approach.For example, in the upper reaches of mountain catchments where accumu-lation zone of different glacier tongues can be connected, the accurate delimitations of each individual glacier can be difficult, mainly due to the lack of information on the bedrock topography under the glacier and on the ice-flux directions.Likewise, catchment delimitation can be hazardous in places with complex topographies dominated by flats (as in South American plateaux) and short-scale steep altitudinal gradients (Verbunt et al., 2003).Second, it may be complicated to determine glacial influence on stream locations as the apparent absence of glacier cover may not be a reliable indicator of an absence of glacial influence on streamflow (Favier et al., 2008).Indeed, in glaciers located on terrains with complex geology and groundwater reservoirs (e.g., volcanoes, karstic areas), meltwater infiltrations are more often the rule than the exception (Bazhev, 1973;Bengtsson, 1982;Favier et al., 2008;Finger et al., 2013).Third, there is growing evidence that water chemical signatures may not be so reliable in detecting ice melt influence on streamflow as they can be modified by many factors such as climate, bedrock substrates and altitude (Nelson et al., 2011;Zhu et al., 2012).In particular, when glacial meltwater infiltration occurs, water chemistry is likely to be considerably modified during the underground flow routing, depending on the residence time underground, the distance of the underground flow routing and the bedrock substrates (Hindshaw et al., 2011;Nelson et al., 2011).Finally, incorporating the high spatiotemporal variability of the different water sources contributions in glacierized catchments requires extensive measurement campaigns (e.g., glacier area measurement, water sampling, and stream habitat measurements), the building of water monitoring structures (e.g., hydrological and climatological stations) or costly analyses (e.g., water chemistry over long time periods).While these factors may not appear as major constraints in temperate regions where many monitoring field stations have been established over the last 50 yr, most glacierized catchments in the world (e.g., subtropical and tropical mountains) remain poorly studied due to the difficulties in accessing and monitoring costs over long time periods (Baraer et al., 2012).In this context, we were interested in developing a new cost-efficient method for detecting the glacial influence in mountain catchments.
During the ablation period, glacier-fed streams are characterized by diurnal flood events (Milner et al., 2009) with discharge depending on the portion of glacier exposed to melting conditions (Favier et al., 2004).In this article, we proposed using the diurnal cycle amplitude as a quantitative measure of glacial influence in hydrosystems located in glacierized catchments.To identify the diurnal flow variation caused by the diurnal glacier melting, we propose using wavelet analyses on water-level time series.Indeed, contrary to Fourier transform and autocorrelogram approaches commonly used for time-series analyses (Chatfield, 1989;Bloomfield, 2004;Andreo et al., 2006), wavelet analysis is a time-dependent spectral analysis that decomposes a data series in time-frequency space and enables to identify  10.0), and catchment basins were defined using SAGA (2.0.8).The location of the Antisana volcano is indicated on the map of Ecuador by a red triangle.The catchments "Los Crespos" (including the catchments of stream sites 7 to 10, 14, and 15) and "Antisana 14 Glacier" (including the catchments of stream sites 1 to 6, and 11 to 13) are delimited by an orange and brown contour line, respectively.repeated events at different temporal scales (Lafreneire and Sharp, 2003).This method has a long tradition in climatology and hydrology (e.g., Smith et al., 1998;Mathevet et al., 2004;Labat, 2005;Jiang et al., 2007), but, surprisingly, only a handful of studies used wavelet analyses on glacier-fed stream discharge time series, with, to our knowledge, only one study (Lafreneire and Sharp, 2003) using wavelet transforms on discharge time series to identify the seasonal and inter-annual variability in the relative contributions of different water sources (e.g., glacial ice, snow, rain and groundwater).
In our study, we used wavelet analyses on water-level time series from 15 stream sites in two glacierized catchments in the tropical Andes of Ecuador to determine the glacial influence in alpine hydrosystems and to identify the fluctuation of this glacial influence throughout the year.We further propose three metrics to quantify the power, the frequency, and the temporal clustering of the diurnal flow variation based on the scale-averaged wavelet power spectrum at 24 h scale.Our goals were (1) to test whether our method was reliable using water-level time series instead of discharge time series; (2) to verify that the diurnal flow variations were only caused by the glacial meltwater; (3) to describe the glacial influence in two well-studied catchments in the tropical Andes using wavelet analyses on water levels; and (4) to test the information provided by the three wavelet metrics (i.e., power, frequency, and temporal clustering of the diurnal flow variation) with regards to one of the most widely used index, the percentage of glacier cover in the catchment.

Climatic, glaciological and hydrological settings
The study was conducted in the Ecological Reserve of Antisana, Ecuador (0 • 29 06 S, 78 • 08 31 W; Fig. 1).From a climatological viewpoint, the Antisana volcano belongs to the inner tropics (sensu Troll, 1941) with more or less continuous precipitation and homogeneous temperature conditions throughout the year (Fig. 2).The Antisana's precipitation regime is complex.Although substantial precipitation is observed year-round, there is always a period with heavy precipitation between February and June.The beginning of this wet season is however quite variable, and another period between September and November generally shows high a amount of precipitation.These features reflect the different origins of precipitation at the Antisana.First, Antisana receives precipitation from the Amazon Basin.The eastern slopes of the Andes are the first obstacles encountered by air masses coming from the east, pushed by the trade winds from the Atlantic (Vuille et al., 2000), creating an ascent of the air and an adiabatic cooling leading to heavy precipitation.Second, the site is located in a border zone with the inter-Andean plateau; thus on Antisana, the precipitation regime of the Amazon regions (a single maximum between June and July and a minimum in February) is mixed with the inter-Andean valley regime (with two wet seasons in February-May and October-November; Vuille et al., 2000).At interannual timescales, there is a general agreement that a significant fraction of the variability of precipitation is related to the El Niño-Southern Oscillation (ENSO) phenomenon, with El Niño years (warm phase of ENSO) tending to be warmer and drier than the average, while La Niña years (ENSO cold phase) are associated with colder and wetter conditions (e.g., Vuille et al., 2000).
From a glaciological point of view, both ablation and accumulation occur all year round on Ecuadorian glaciers (Favier et al., 2004;Francou et al., 2004;Rabatel et al., 2013).Moreover, in such tropical regions, there is neither permanent nor seasonal snow cover outside the glaciers due to a combination of two characteristics.First, the 0 • C isotherm is located around 4950-5000 m a.s.l.throughout the year as there is no seasonality in temperature in Ecuador (see Fig. 2; Vuille et al., 2008).Second, glacier snouts are located at about 4850-4900 m a.s.l.(Rabatel et al., 2013).Thus, precipitation outside the glaciers is almost exclusively liquid, except during exceptionally cold conditions during strong La Niña events.
From a hydrological viewpoint, the three main components of streamflow are direct superficial runoff, snow and ice melting, and groundwater.The mean monthly discharge ranges from 0.04 to 0.1 m 3 s −1 at Los Crespos station (1.6 km from the glacier snout of the "Crespo" glacier) and from 0.25 to 0.3 m 3 s −1 at Humboldt station (8 km from the glacier snout; see Fig. 3).The differences in absolute values of outflows are due to the different drainage areas with 2.4 and 14.2 km 2 , respectively.Two different patterns in the monthly outflows variations can be observed.The mean monthly discharge for the Los Crespos station shows a perennial flow with the lower values observed between June to August and higher values from October to May.The high discharge values are a consequence of important glacier melting resulting from strong shortwave radiation absorption (i.e., low surface albedo, resulting from liquid precipitation over the lower part of the glacier).Low discharge values are a consequence of higher wind velocity that enhances mass losses through sublimation instead of melting.The correlation between the precipitation and the outflows is weak, and the regime is mainly controlled by glacier melting (Favier et al., 2004).The flow at the Humboldt station shows low seasonal variations in accordance to the pluviometric regime with glacier contribution during the months of lower precipitation.

Stream sites
The study was conducted in 15 stream sites belonging to two glacier-fed catchments in the Ecological Reserve of Antisana.These sites were selected because detailed long-term climatic and hydrological data were available in this area.For example, work on the geology and water chemistry of these catchments (see Favier et al., 2008;Villacis, 2008) has detected glacial meltwater infiltrations and resurgences, which were of great interest to test the relevance of our method (see Introduction).
To assess the wavelet signals of a broad range of glacial contributions, 10 of the 15 stream sites (no.1-10; Fig. 1) were localized along two glacier-fed streams and five on their respective tributaries (no.11-15; Fig. 1).Among the studied tributaries, four stream sites (no.11-14) were considered non-glacial as they had no glacier cover in their catchment (see Fig. 1) and did not present physico-chemical features generally observed in glacier-fed streams (Brown et al., 2003), e.g., high turbidity (> 30 NTU), low conductivity (< 10 µS cm −1 ) (see Table 1), and one site (no.15) was partially fed by glacial meltwater (glacier cover in the catchment = 1.0 %).The two glacier-fed streams originated at 4730 m a.s.l.: one from the snout of the "Crespos" Glacier, which covered an area of about 1.82 km 2 at the time of the study in 2010, and the second from the snout of the Glacier "14", which covered an area of about 1.24 km 2 (Rabatel et al., 2013).Stream sites were all located between 4040 and 4200 m a.s.l., between 5.9 and 9.6 km away from the glacier snouts.At the glacier snout, both streams had high turbidity (> 285 NTU) and low conductivity (< 9 µS cm −1 ), which decreased and increased downstream, respectively, in particular after confluences with tributaries (see Jacobsen et al., 2010;Kuhn et al., 2011, for details).Contrastingly, non-glacial tributaries presented high conductivity (> 60 µS cm −1 ) and low turbidity (< 10 NTU, see Table 1 for details on the physicochemical characteristics of each stream site).
Table 1.Physico-chemical attributes of the study stream sites (see location of the sites on Fig. 1).Conductivity, turbidity, and temperature are means with min-max values given in brackets (n = 2 to 10 for conductivity and temperatures; n = 1 to 3 for turbidity).Stream sites no.11 to 15 have no visible connection to the glacier.GCC = glacier cover in the catchment.UTM coordinates refer to stream site coordinates in UTM-WGS84 zone 17S expressed in meters.

Field measurements
The location of each stream site was measured using the UTM-WGS84 coordinate system with a GPS (Garmin Oregon 550, Garmin International Inc., Olathe, USA).In January 2010, 15 water pressure loggers were installed (Hobo water pressure loggers, Onset Computer Corp., USA) in the water at each stream site and recorded water pressure every 30 min over 10 months (i.e., from January to October 2010).Water pressure loggers were previously protected in plastic tubes placed vertically on the stream side where the sections were deep enough to avoid overflowing during the glacial flood and with homogeneous shapes among stream sites.
Water level and height between the stream bottom and the Hobo sensor were measured twice, when the loggers were installed and removed.One more logger was fixed on a rock at 4100 m a.s.l. to measure the atmospheric pressure and the air temperature every 30 min over the same 10-month period.
In addition, at the Los Crespos hydrological runoff gauging station (close to site 7), discharge was recorded every 30 min during the entire year of 2010 (Fig. 4).Precipitation was also recorded every 30 min at the weather station, located on the glacier foreland of Antisana 15 Glacier, within our catchment area.
3 Materials and methods

Wavelet transform analyses on water-level time series
Our method proposes using wavelet transform analyses on water-level time series to detect the hydrological signal originating from glacier melting.As glacial runoff exhibits repeated cyclic fluctuations at the daily timescale during the ablation period (Hannah et al., 1999(Hannah et al., , 2000)), we aim to detect corresponding variations in water level at 24 h scale.Previous work has reviewed in detail the concepts of wavelet analysis for different applications (Daubechies, 1990;Torrence and Compo, 1998;Cazelles et al., 2008).Here we list some important concepts with special attention to properties used in this study.The wavelet transform analysis is a time-dependent spectral analysis that decomposes a data series in time-frequency space.The wavelet transforms therefore express a time series in a three-dimensional space: time (x), scale/frequency (y), and power (z).The power matches the magnitude of the variance in the series at a given wavelet scale and time.Various types of wavelet functions (e.g., Morlet, Mexican hat, Paul) can be used for the signal transform, depending on the nature of the time series and the objectives of the study.Here, we chose the Morlet wavelet, a nonorthogonal, continuous, and complex wavelet function (with real and imaginary parts), because it is particularly well adapted for hydrological time-series analyses (Torrence and Compo, 1998;Labat et al., 2000;Lafreneire and Sharp, 2003).Nonorthogonal continuous wavelet transforms are indeed more robust to noise than other decomposition schemes and are robust to variations in data length (Cazelles et al., 2008).Moreover, complex wavelet functions are well suited for capturing oscillatory behavior, whereas real wavelet functions do better at isolating peaks or discontinuities (Torrence and Compo, 1998).Finally, the Morlet wavelet function has a high resolution in frequency compared to other continuous wavelets (Cazelles et al., 2008), which was fundamental in our method as we intended to detect the repeated water-level variations at 24 h scale.
The continuous wavelet transform W n (s) of a discrete time series x n (n being the time position) at scale s is defined as the convolution of x n with a scaled and translated version of the wavelet function ψ(t): where N is the number of points in the time series, ψ * (t) the complex conjugate of wavelet function (the Morlet wavelet in our case) at scale s and translated in time by n, and δt the time step for the analysis (Torrence and Compo, 1998).The Morlet wavelet is defined as where t is a nondimensional "time" parameter, and w 0 the nondimensional frequency.w 0 must be equal to or greater than 5 to satisfy the wavelet admissibility condition (Farge, 1992;Cromwell, 2001): the function must have zero mean and be localized in both time and frequency space to be "admissible" as a wavelet (Santos et al., 2001).In the present application, we use w 0 = 6, a value commonly used for geoscience applications (Torrence andCompo, 1998, Labat, 2005;Si and Zeleke, 2005;Schaefli and Zehe, 2009).For w 0 = 6, the Morlet wavelet scale was almost identical to the corresponding Fourier period of the complex exponential, and the terms "scale" and "period" may conveniently be used synonymously (Torrence and Compo, 1998;Torrence and Webster, 1999;Maraun and Kurths, 2004).Thus the left axis in Figs. 5 and 6 is the equivalent Fourier period corresponding to the wavelet scale.
To visualize the magnitude of the variance in the series at a given wavelet scale and location in time, we determined the local wavelet power spectrum (Torrence and Compo, 1998), defined as the squared absolute value of the wavelet transform (|W n (s)| 2 ) and calculated as follows: where W n * (s) is the complex conjugate of W n (s).
To examine fluctuations in power over a range of scales (a band), we calculated the scale-averaged wavelet power spectra W 2 n , defined as the weighted sum of the wavelet power spectrum (over two scales s 1 to s 2 ), over the whole measurement period.
where δj is the spacing between discrete scales, δt the time step of the time series, and C δ the reconstruction factor (Anctil and Coulibaly, 2004;Coulibaly and Burn, 2004;Markovic and Koch, 2005;White et al., 2005).

Glacial signal determination
As glacier-fed streams are characterized by diurnal flood events, with discharge depending on the portion of glacier exposed to melting conditions (Favier et al., 2004;Villacis, 2008), we calculated the scale-averaged wavelet power spectrum over scales around 24 h over the whole time series, which permitted visualizing the fluctuation of the diurnal flow variation power throughout the year.The significance of the wavelet power spectrum was tested against a background (or noise) spectrum, which is either white noise (constant variance across all scales or frequencies) or red noise (increasing variance with increasing scale or decreasing frequency; Schiff, 1992).When the wavelet power of the time series exceeds the power of the background (at the chosen confidence level), the time-series variance can be deemed significant (see Torrence and Compo, 1998;Lafreneire and Sharp, 2003, for background spectrum calculation details).Assuming a random process, we chose here the red-noise spectrum (at 95 % confidence level; see Lafreneire andSharp 2003, Schaefli et al., 2007, for a justification on red-noise choice).
To express the main features of the scale-averaged wavelet power spectrum quantitatively over time, allowing a comparison of glacial influence among sites, we defined three metrics: (1) the diurnal variation power, (2) the diurnal variation frequency, and (3) the diurnal variation temporal clustering.The scale-averaged wavelet power at 24 h scale is considered significant when above the 95 % confidence level curve (see Torrence and Compo, 1998;Markovic and Koch, 2005).We calculated the diurnal variation power as the integration of the corrected scale-averaged wavelet power curve (i.e., divided by its corresponding 95 % confidence level value).Note that we calculated the area under the curve above the "y = 1" line (see Torrence and Compo, 1998).Higher values of diurnal variation power correspond to higher power of the variance in water levels at 24 h scale.The diurnal variation frequency was calculated as the frequency of days with significant diurnal flow variations in the time series.A diurnal variation frequency equal to 0 means no significant diurnal flow variation over the study period while when equal to 1 it corresponds to a significant diurnal flow variation every day.Concerning the diurnal variation temporal clustering (sensu De Vos et al., 2010;Hsu and Li, 2010), we first defined two "hydrological states" corresponding to days with and without significant diurnal flow variation.We then calculated the number of hydrological state changes and divided it by the total number of days in the time series minus one (the maximum number of possible state modifications).If the diurnal variation temporal clustering is equal to 1, there is no hydrological state shift (i.e., there is never/or every day significant diurnal flow variations).On the contrary, if it is equal to 0, the hydrological state of the stream changes every day.

Local and scale-averaged wavelet power spectrum
To determine the glacial influence based on water-level time series, we first transformed water pressure values, obtained from the 15 stream sites, into water-level values by subtracting the atmospheric variations from the water pressure data.Time series were centered on their means and normalized by their standard deviations prior to wavelet transform calculation to allow across-site comparison of our results.We then developed a code inspired by C. Torrence and G. Compo (available at http://paos.colorado.edu/research/wavelets)that we ran in MATLAB, version R2009a (The MathWorks Inc., Natick, MA, USA).This code allowed producing three types of figures: (1) the averaged normalized water-level time series, which presents the water-level variations throughout the year; (2) the local wavelet power spectrum (normalized by their standard deviations), which gives the magnitude of the variance in the series at a given wavelet scale and location time; and (3) the scale-averaged wavelet power spectrum, which presents the fluctuation in power over 24 h scale over a whole year.For illustrative purposes, Fig. 6 presents the outputs of our wavelet analyses for three stream sites with contrasting time-series patterns resulting from different glacial influences: a glacier-fed stream site without tributaries (no.7, Fig. 6a, d and g) and two groundwater-fed stream sites (no.13 and 14, Fig. 6b, c, e, f, h, and i).For all stream sites we then calculated from the scale-averaged wavelet power spectrum the three metrics defined above (diurnal variation power, diurnal variation frequency, and diurnal variation temporal clustering).
To test whether our method was reliable using water-level time series instead of discharge time series, we compared wavelet outputs of water-level time series at site 7 and those from discharge time series at the Los Crespos hydrological station (Fig. 4).Moreover, to exclude the rainfall as contributor of diurnal flow variation, we plotted the precipitation and air temperature time series (Fig. 5a and b), their corresponding local wavelet power spectrum (Fig. 5c and d), and scale-averaged wavelet power spectrum at 24 h scale (Fig. 5e  and f).

Glacier metrics vs. glacier cover in the catchment
As the percentage of glacier cover in the catchment is commonly used to estimate the potential influence of a glacier on a stream (e.g., Hari et al., 2005;Thayyen and Gergan, 2010;Jacobsen et al., 2012), we were interested in assessing how our three wavelet metrics at each stream site would behave as a function of glacier cover in the catchment.
To measure the percentage of glacier cover in the catchment, we first created the channel network, and the catch-ment area was calculated with the SAGA GIS software (System for Automated Geoscientific Analyses, version 2.0.8).Briefly, SAGA derives a channel network based on gridded digital elevation model (DEM) with the specification of the target cells (gauge station), for which the upslope contributing area is identified.The catchment delimitation is based on the multiple flow direction model (Tarboton, 1997), and the extraction of the drainage network uses the algorithm described in O' Callaghan and Mark (1984).We created the hydraulic channel network of our two study catchments using a 40 m resolution DEM in SAGA GIS.The DEM was created using 40 m resolution contour line from the Ecuadorian Military Geographical Institute (available at http://www.igm.gob.ec/site/index.php) in ArcGIS (10.0).We verified each created channel with our GPS point measurements and field observations, and determined for each stream site the corresponding catchment using SAGA GIS (Olaya and Conrad, 2009).Hereafter we named "Los Crespos" and "Antisana 14 Glacier" catchments as the area of all small catchments they contain: Los Crespos catchment includes catchments of stream sites 7 to 10, 14, and 15; and Antisana 14 Glacier catchment includes catchments of stream sites 1 to 6, and 11 to 13 (see Fig. 1).
Glacier outlines were first automatically extracted from LANDSAT satellite images (30 m pixel size) using the Hydrol.Earth Syst.Sci., 17, 4803-4816, 2013 www.hydrol-earth-syst-sci.net/17/4803/2013/ common Normalized Difference Snow Index (NDSI).The glacier outlines were then manually checked and adjusted by overlaying the glacier outline shapefile on the satellite images for which a spectral band combination associating the shortwave infrared, the near infrared, and the green bands had been applied (such a combination facilitates the distinction between snow, ice and rock; see Fig. 4 in Rabatel et al., 2012).We finally merged the glacier outlines and catchment contours shapefiles using ArcGIS 10.0.This enabled computing the percentage of the glacier cover in the catchment basin of each stream site by dividing the glacier area by the total catchment basin area.The three metrics (diurnal variation power, diurnal variation frequency, and diurnal variation temporal clustering) were plotted against the percentage of glacier cover in the catchment.The strength of the glacier cover in the catchment vs. the three metrics correlation was measured using Spearman correlation coefficients and associated p values.

Water-level vs. discharge time series
To verify whether our method was reliable using water-level time series, we compared our wavelet analysis outputs on discharge time series (in "Los Crespos" hydrological station) vs. those obtained with water-level time series at site 7.We found a very good match between both curves (Fig. 4), and all significant diurnal flow variations on discharge data were also detected on water-level data.This indicates that our method for detecting and characterizing glacial influence can be applied on both stream discharge and water-level time series.This is an interesting result as using water-level data instead of stream discharge presents at least two main advantages.First, water-level data are easily obtained using staff gauges or pressure transducers, whereas discharge data additionally require detailed velocity measurements under a range of conditions to generate a rating curve.Second, discharge is usually inferred from stage measurements using a rating curve, which is subject to error.However, waterlevel data can be used only if the stage-discharge relationship is not too strongly non-linear, and if the shape of the cross sections remains unchanged throughout the measurement period, which was the case in our study.

Sorting out the contribution of rainfall and snowmelt to diurnal flow variation
To verify that diurnal flow variations were mainly caused by glacier or snow melting, we first performed wavelet analyses on air temperature time series.The local wavelet power spectra of air temperature time series show that most of the variance in temperature time series indeed occurred at 24 h scale (Fig. 5d and f).We then performed wavelet analyses on precipitation time series to rule out the potential influence of diurnal rainfall events (e.g., convective storm activity) on diurnal flow variation.We only found seven significant peaks over 2010 in the scale-averaged wavelet power spectrum at 24 h scale (Fig. 5e), which corresponded to a total of only 19 rainy days with significant diurnal variation.This supports the fact that, in our case, the diurnal variation of the streamflow could be employed as an indicator of meltwater influence (Favier et al., 2004).Note that if the local wavelet spectrum for precipitation data had significant power at 24 h scale all year round (which was not the case in our study), a wavelet coherence spectrum analysis on the two time series could then be run (Torrence and Compo, 1998;Maraun and Kurths, 2004), and would allow determining whether or not the water-level-precipitation cross spectrum mimics the general pattern observed in the wavelet spectrum of water-level time series.
In the inner tropics, the absence of snow cover outside the glaciers (see study site section) ensures that glacier melting is the main cause of diurnal flow variations.Note, however, that in other regions, our method would not allow identifying the relative contribution of ice melt and snowmelt on streamflow variation, and should therefore be employed as an indicator of seasonal or multi-annual storage of water in the form of snow or ice.

Local and scale-averaged wavelet power spectrum to detect glacial influence
Most of the variance in the local wavelet spectra of stream sites 1 to 10 was concentrated at the 24 h period (Fig. 6d).This diurnal wavelet power represented the diurnal glacial flood (see also Lafreneire and Sharp, 2003), resulting from the diurnal fusion of the ice (see above).At all sites along the two glacier-fed streams, the power of the local wavelet spectra at the 24 h period was statistically significant over the 10 months (Fig. 6d).Indeed, glacial floods occur all year round in equatorial glacier-fed streams due to the lack of thermal seasonality in the inner tropics (Favier et al., 2008;Vuille et al., 2008;Jacobsen et al., 2010).However, we found that, for all stream sites from 1 to 10, the scale-averaged wavelet power at 24 h scale was continuously significant between January and May while it was seldom significant after May (e.g., site 7, Fig. 6g).This phenomenon is related to strongest ablation rates of the glaciers during this period (January to May).Indeed, Rabatel et al. (2013) showed that during the period from January to April 2010 the Antisana glaciers experienced high ablation rates related to El Niño conditions (see Fig. 9 in Rabatel et al., 2013), which favored glacier melting.Contrastingly, ablation is generally reduced under La Niña conditions (Francou et al., 2004).As expected, the three non-glacial tributaries (sites 11, 12, and 13) did not present any significant power at 24 h scale (Fig. 6f).The scale-averaged wavelet power spectrum was below the 95 % confidence level curve during the entire study period (Fig. 6i).Surprisingly, significant glacial signal was identified at one supposedly non-glacial site (site 14, Fig. 6e).While this site presented no glacier cover in its catchment as well as non-glacial characteristics (turbidity < 9 NTU, conductivity > 90 µS cm −1 ; Fig. 1 and Table 1), our wavelet analysis revealed a significant glacial influence (Fig. 6e and h).However, the unexpected spectrum at site 14 confirms the presence of glacial water resurgence, as previously detected by Villacis (2008) and Favier et al. (2008) (see study site section).Site 14 was a spring with a substantial discharge (around 20 Ls −1 , similar to that of site 7) revealing the presence of an aquifer with meltwater alimentation and underground flow transfer.The scale-averaged wavelet power spectrum at 24 h scale allowed identifying such infiltrations using only data of water levels, making the diurnal flow variation analysis a good method for understanding glacial influence.This is an important result as catchments with complex geological structure containing groundwater reservoirs and meltwater infiltrations are common throughout the world, in particular in volcanic regions (Favier et al., 2008).

Relationship between wavelet metrics and glacier cover in the catchment
As expected, we found an overall significant positive relationship between the diurnal variation power and the percentage of glacier cover in the catchment (Spearman rank test, r = 0.71, F = 13.12,p < 0.01).The diurnal variation power generally decreased as moving downstream, and this decrease was more pronounced at sites located after a confluence with a groundwater tributary (except for site 14), a pattern already observed by Jacobsen and Dangles (2012) using their glacier index.Figure 7a also highlights the particularity of site 14, whose diurnal variation power was far above the regression line.Note that removing site 14 from the relationship between the diurnal variation power and percentage of glacier cover in the catchment increased markedly the correlation coefficient (r = 0.93, F = 76.14, p < 0.001).This confirms that analyzing the diurnal flow variation may be a much better alternative to the percentage of glacier cover in the catchment as the latter does not permit the detection of glacial meltwater reemergence.We found a high diurnal variation frequency (> 0.45, i.e., the diurnal flow variation was significant about half of the year) for all sites that had at least part of their catchment covered by glacier.Interestingly, while we found a significant positive relationship between the diurnal variation frequency and percentage of glacier cover in the catchment in the Antisana 14 Glacier catchment (r = 0.98, F = 194.98,p < 0.001, Fig. 7b), this was not the case for the Los Crespos catchment (no significant relationship).This suggests that the two catchments have different hydrological behaviors, which may be related to the different origin of the tributaries of the main glacier-fed stream.Indeed, while the Los Crespos stream had several tributaries with glacial influence (e.g., site 14, 15), the Antisana 14 Glacier catchment had only non-glacial tributaries (see Fig. 1).
Stream sites with no glacier cover in their catchment had diurnal variation temporal clustering values equal to 1 as they were always in the same hydrological state of no significant diurnal flow variation (full clustering; Fig. 7c).All other stream sites (e.g., with glacier cover in the catchment > 0 %) had also high temporal clustering values (> 0.85) meaning that the two possible hydrological states (days presenting either significant or non-significant diurnal flow variation) did not alternate frequently but were rather clustered over the year (i.e., the ablation process concentrates at specific periods of the year).For these stream sites, we found a significant positive relationship between the diurnal variation temporal clustering and the percentage of glacier cover in the catchment (r = 0.76, F = 12.33, p < 0.01, see Fig. 7c).This suggests a decreasing number of switches between the two hydrological states as percentage of glacier cover in the catchment increases.Due to the high diurnal variation frequency measured in glacier-fed sites (see above), the probability of having alternative hydrological states decreases closer to the glacier, where significant diurnal flow variations are more frequent.Note that both diurnal variation frequency and temporal clustering do not have an absolute meaning, but they give substantive additional temporal information to characterize the glacial influence, and are used to compare signals among sites.
The three metrics obtained with our wavelet analyses are complementary and do not convey the same information on the hydrology of the watersheds.Indeed each metric presented a different relationship with the percentage of glacier cover in the catchment (Fig. 7a, b, and c).Even though the diurnal variation power was sufficient to detect, quantify and compare the glacial influence among sites, it did not give any information about the variation over time of this glacial influence.Both diurnal variation frequency and temporal clustering revealed that the temporal dynamics of water regimes differed between the two study catchments, and that the high flow/low flow regime shift could be a relevant index of glacier behavior.As the variance in hydrological regime of glacier-fed streams in the tropics has been proposed to be a good indicator of the state of glacier retreat (Baraer et al., 2012), further applications of the diurnal variation temporal clustering index may be relevant.

Conclusions
Our study is the first to propose the analysis of diurnal flow variation as an indicator of glacial influence in mountain catchments.Opting for wavelet analyses as a methodological framework had several advantages as it allowed a full description of the glacial influence over time.We showed that our method (1) is relatively simple and cost-effective as it can be used with water-level data instead of discharge (3) reliably quantifies the glacial influence and can therefore be used as a substitute of the percentage of glacier cover in the catchment; (4) allows detecting the influence of glacial meltwater in resurgence; and (5) allows providing a detailed description of temporal patterns in flow variations over time thereby making possible quantitative comparisons of the glacial influence among sites.Taking all these aspects in consideration, we think that our method represents a significant improvement when compared to existing methods, as it overcomes most of their limitations (see Introduction).
Several issues however require further investigation.In particular, it would be interesting to study whether our method could be used on continuous stream discharge time series to quantify (in terms of volume) the glacier contribution of different water sources at a given stream site.Another important issue is to what extent this method may be applicable in other regions of the world.On the one hand, as diurnal glacial floods always occur during the ablation season in any glacier-fed streams worldwide -e.g., in the tropical Andes (Rabatel et al., 2013), the Himalayas (Sorg et al., 2012), the European Alps (Schutz et al., 2001), the North American Rockies (Lafreneire and Sharp, 2003), and the Arctic (Dahlke et al., 2012) -the core of our approach should be valid.On the other hand, several refinements would be needed to account for the specificities of temperate and arctic regions such as the presence of snow cover outside the glacier or the potential interference of rainfall with diurnal cycles, and strong annual cycle.Despite these limitations, our hope is that our method may provide a testable and applicable methodological framework to understand better the complex interactions between glacier and glacier-fed hydrosystems in a warming world.

Fig. 1 .
Fig. 1.Study area at the Antisana volcano, Ecuador.Stream sites are represented by orange circles.Catchment basins of all stream sites are represented by blue (glacier-fed stream sites) and green (supposedly non-glacial tributaries) polygons.Map was made using ArcGIS (10.0), and catchment basins were defined using SAGA (2.0.8).The location of the Antisana volcano is indicated on the map of Ecuador by a red triangle.The catchments "Los Crespos" (including the catchments of stream sites 7 to 10, 14, and 15) and "Antisana 14 Glacier" (including the catchments of stream sites 1 to 6, and 11 to 13) are delimited by an orange and brown contour line, respectively.

Fig. 2 .
Fig. 2. Mean monthly precipitation (bars) and temperatures (dots) over 1 yr at the Antisana volcano, Ecuador.The weather station is located on the proglacier margin of Glacier 15 at 4850 m.Mean values and standard deviations were calculated over 6 yr (2005-2010).

Fig. 3 .
Fig. 3. Mean monthly discharge data at two gauging stations with high (A, Los Crespos station, 54.5 % of glacier cover in the catchment) and low glacial influence (B, Humboldt station, 8.6 % of glacier cover in the catchment) over 1 yr in the Antisana volcano catchment, Ecuador.Mean values and standard deviations for (A) and (B) were calculated over 5 (2006-2010) and 11 yr (2000-2010), respectively.

Fig. 4 .
Fig. 4. Comparison of the wavelet analysis output (scale-averaged power spectrum over 24 h scale) between discharge (black line) and water-level time series (blue line) at site 7 (Los Crespos station).Dashed black and blue lines indicate the corresponding 95 % confidence levels for the red-noise spectrum for discharge and waterlevel time series, respectively.

Fig. 5 .
Fig. 5. Wavelet analysis outputs for rainfall (A, C, and E) and air temperature time series (B, D, and F): (A-B) averaged normalized rainfall (A), and air temperature (B) time series.(C-D) local wavelet power spectrum normalized by their standard deviations.The black line delineates the areas where the power is significant (i.e., exceeds the 95 % confidence level of a red-noise process).The dashed black line delineates the cone of influence that delimits the region not influenced by edge effects.(E-F) normalized scale-averaged wavelet power spectra at 24 h.The dashed blue line shows the corresponding 95 % confidence level for the red-noise spectrum.In each panel, day one corresponds to the 1 January 2010.The color bar shows the legend for the different colors: blue and red for low and high wavelet power, respectively.

Fig. 6 .
Fig.6.Wavelet analysis outputs at three stream sites(7, 13, 14)  with contrasting glacial influence (site 7 is glacier-fed while site 13 and 14 have no superficial connection to the glacier -but site 14 has glacial resurgence as indicated by significant diurnal flow variation).(A), (B), (C): averaged normalized water-level time series.(D), (E), (F): local wavelet power spectrum normalized by their standard deviations.The black line delineates the areas where the power is considered significant (i.e., exceeds the 95 % confidence level of a red-noise process); the dashed black line delineates the cone of influence that delimits the region not influenced by edge effects.(G), (H), (I): normalized scaleaveraged wavelet power spectra at 24 h.The dashed blue line shows the corresponding 95 % confidence level for the red-noise spectrum.In each panel, day one corresponds to the 1 January 2010.The color bar shows the legend for the different colors: blue and red for low and high wavelet power, respectively.

Fig. 7 .
Fig. 7. Scatter plot of (A) diurnal variation power, (B) diurnal variation frequency, and (C) diurnal variation temporal clustering vs. the percentage of glacier cover in the catchment.In (B) black dots represent stream sites located in the Los Crespos catchment and open dots stream sites in the Antisana 14 Glacier catchment.The full regression line excludes site 14 in (A), sites localized in the Los Crespos catchment (i.e., excludes black dots) in (B) and sites with no glacier cover in the catchment (i.e., excludes sites 11 to 14) in (C).Sites numbers are indicated in the three panels.