Insights into the water mean transit time in a high-elevation tropical ecosystem

Abstract

This study focuses on the investigation of the mean transit time (MTT) of water and its spatial variability in a tropical high-elevation ecosystem (wet Andean páramo). The study site is the Zhurucay River Ecohydrological Observatory (7.53 km${}^{\mathrm{2}})$ located in southern Ecuador. A lumped parameter model considering five transit time distribution (TTD) functions was used to estimate MTTs under steady-state conditions (i.e., baseflow MTT). We used a unique data set of the ${\mathit{\delta }}^{18}$O isotopic composition of rainfall and streamflow water samples collected for 3 years (May 2011 to May 2014) in a nested monitoring system of streams. Linear regression between MTT and landscape (soil and vegetation cover, geology, and topography) and hydrometric (runoff coefficient and specific discharge rates) variables was used to explore controls on MTT variability, as well as mean electrical conductivity (MEC) as a possible proxy for MTT. Results revealed that the exponential TTD function best describes the hydrology of the site, indicating a relatively simple transition from rainfall water to the streams through the organic horizon of the wet páramo soils. MTT of the streams is relatively short (0.15–0.73 years, 53–264 days). Regression analysis revealed a negative correlation between the catchment's average slope and MTT ($R^{\mathrm{2}}=$ 0.78, $p$ < 0.05). MTT showed no significant correlation with hydrometric variables, whereas MEC increases with MTT ($R^{\mathrm{2}}=$ 0.89,$p$ < 0.001). Overall, we conclude that (1) baseflow MTT confirms that the hydrology of the ecosystem is dominated by shallow subsurface flow; (2) the interplay between the high storage capacity of the wet páramo soils and the slope of the catchments provides the ecosystem with high regulation capacity; and (3) MEC is an efficient predictor of MTT variability in this system of catchments with relatively homogeneous geology.

Introduction

Investigating ecohydrological processes through the identification of fundamental catchment descriptors, such as the mean transit time (MTT), specific discharge, evapotranspiration to precipitation ratios, and others, is fundamental in order to (1) advance global hydrological, ecological, and geochemical process understanding and (2) improve the management of water resources. This is particularly critical in high-elevation tropical environments, such as the wet Andean páramo (further referred to as “páramo”), in which hydrological knowledge remains limited, despite its importance as the major water provider for millions of people in the region (De Bièvre and Calle, 2011; IUCN, 2002). Water originating from the páramo sustains the socio-economic development in this region by fulfilling urban, agricultural, industrial, and hydropower generation water needs (Célleri and Feyen, 2009). Here we focus on the MTT of water, which we define as the average time elapsed since a water molecule enters a catchment as recharge to when it exits it at some discharge point (Bethke and Johnson, 2002; Etcheverry and Perrochet, 2000; Rodhe et al., 1996).

Despite the importance of tropical biomes as natural sources and regulators of streamflow, there are very few studies of MTT in tropical environments (e.g., Farrick and Branfireun, 2015; Muñoz-Villers et al., 2016; Roa-García and Weiler, 2010; Timbe et al., 2014). The majority of MTT studies have been conducted in catchments with strong climate seasonality, i.e., located in the Northern and Southern hemispheres (e.g., Lyon et al., 2010; McGlynn and McDonnell, 2003; McGuire et al., 2005), and considerably less attention has been devoted to tropical environments. Most tracer-based studies conducted in tropical latitudes focused on isotope hydrograph separation at storm event scale (e.g., Goller et al., 2005; Muñoz-Villers and McDonnell, 2012), the isotopic characterization of precipitation patterns (e.g., Vimeux et al., 2011; Windhorst et al., 2013), and the identification of ecohydrological processes (e.g., Crespo et al., 2012; Goldsmith et al., 2012; Mosquera et al., 2016). However, studies focusing on MTTs in order to improve the understanding of rainfall–runoff processes and their dependence on landscape biophysical features in tropical regions are still lacking and urgently needed in order to improve water resource management.

The páramo is widely recognized by its high runoff regulation capacity (i.e., páramo water yield is sustained year-round regardless of precipitation inputs) (Buytaert et al., 2006; Célleri and Feyen, 2009). However, efforts to investigate the processes that control such hydrological behavior are rare. Recent investigations in our study site suggest that runoff originates from the shallow organic horizon of the páramo soils located near the streams (Histosol soils or Andean wetlands), thus favoring shallow subsurface flow. By contrast, deep groundwater contributions to discharge are minimal and saturation excess overland flow (even in the nearly saturated Histosol soils) rarely occurs (Buytaert and Beven, 2011; Crespo et al., 2011). The hydrological importance of shallow subsurface flow to runoff generation has also been demonstrated in a variety of ecosystems around the globe (e.g., Freeze, 1972; Hewlett, 1961; Penna et al., 2011), but yet, MTTs have not been explored in these systems. Our study site provides a unique opportunity to gain understanding of the MTT of a shallow subsurface flow dominated system in a tropical setting. In addition, the study of the MTT in natural wetland systems has been limited to sites located in northern boreal catchments in Sweden (Lyon et al., 2010) and peatlands in Scottish mountainous regions (e.g., Hrachowitz et al., 2009a; Tetzlaff et al., 2014). While streams draining these catchments have significant contributions from spring snowmelt and groundwater, respectively, neither of these allows for the isolation of the effect of wetlands in the shallow subsurface and/or near-surface (i.e., overland flow) transport of the water within the catchments.

Another critical issue is the identification of controls on MTT variability among catchments. As detailed observations of combined hydrometric and isotopic information are not feasible in many regions due to limited funding and site accessibility, identifying controls on MTT variability in nested and paired monitoring systems of streams is fundamental towards regionalization of ecohydrological processes (Hrachowitz et al., 2009a) and prediction in ungauged basins (Tetzlaff et al., 2010). However, investigation of controls on MTT variability is still fairly scarce (Tetzlaff et al., 2013). Most studies have found that MTT scales with topographic and/or hydropedological controls. For instance, topographical controls on MTT variability were found in New Mexico, USA, catchments (slope direction and exposure) (Broxton et al., 2009) and a system of streams in Oregon, USA (ratio between flow path length and flow path gradient) (McGuire et al., 2005), whereas the proportions of wetlands and responsive soils were reported as major MTT controls in Swedish catchments (Lyon et al., 2010) and Scottish streams (Soulsby et al., 2006), respectively.

In the last few decades, MTT modeling has been conducted by applying an approach that assumes steady-state conditions in the hydrologic systems, i.e., the lumped convolution approach (LCA) (Amin and Campana, 1996; Małoszewski and Zuber, 1996). However, this assumption is often violated as a result of precipitation variability, strong climate seasonality, and heterogeneities in the landscape configuration and hydropedological distribution of catchments (Kirchner, 2016a, b). This has led to a growing recognition of the time-variant nature of transit times (e.g., Birkel et al., 2015; Harman, 2015; Hrachowitz et al., 2013), which a number of studies have begun to explore. Most of this initial work has yielded results with a high degree of uncertainty (Harman, 2015; Klaus et al., 2015; McMillan et al., 2012), or uncertainty has not been estimated (Davies et al., 2013; Heidbüchel et al., 2012; van der Velde et al., 2015) as a result of expensive computation costs or high uncertainties related to the spatial variability of the input hydrometric and tracer field measurements. It is clear, however, that given the mathematical limitations in representing systems' non-stationarities (Duvert et al., 2016; Seeger and Weiler, 2014), the high temporal resolution of tracer data (Harman, 2015; Heidbüchel et al., 2012), and the general unavailability of long-term tracer records (Hrachowitz et al., 2010; Klaus et al., 2015) required for hydrological modeling under non-stationary conditions, the LCA remains a useful methodology for MTT estimation. This holds true not only in understudied regions such as the tropics (e.g., Farrick and Branfireun, 2015; Muñoz-Villers et al., 2016; Timbe et al., 2014), but also elsewhere (e.g., Duvert et al., 2016; Hale and McDonnell, 2016; Hale et al., 2016; Hu et al., 2015; Seeger and Weiler, 2014).

Given the relatively homogenous landscape features and very low seasonality in the hydrometeorological conditions in our páramo site, applying the LCA to our unique water stable isotopic data set represents a robust first step to build up improved catchment functioning understanding using hydrometric-tracer-based hydrologic modeling. To our knowledge, this is the first contribution regarding the modeling of MTT in páramo ecosystems, and more generally, in regions with low climate seasonality and catchments with a low degree of heterogeneity. Future efforts building on the monitoring infrastructure and continuously collected data sets will allow for continual improvement in hydrologic interpretation, eventually incorporating alternative modeling techniques that explicitly recognize non-stationarity and storage dynamics in the hydrological behavior of this tropical ecosystem (e.g., Birkel et al., 2015; Harman, 2015; Hrachowitz et al., 2013).

In this study, we seek to add to the current geographical scope of MTT studies by addressing two questions that remain open in hydrological science and have received little attention in high-elevation tropical ecosystems: “How old is stream water?” (McDonnell et al., 2010) and “How does landscape structure influence catchment transit time across different geomorphic provinces?” (Tetzlaff et al., 2009). Detailed hydrometric observations highlighting subsurface dominated rainfall–runoff response (Crespo et al., 2011; Mosquera et al., 2016) together with information on the landscape biophysical characteristics in our páramo study site will allow for process-based understanding regarding (i) the spatial variability of baseflow MTTs and (ii) the factors controlling such variability. Based on our current knowledge of the hydrology of the ecosystem, in particular, the apparent dominance of shallow subsurface flow to runoff generation, we hypothesize relatively short baseflow MTTs compared to systems dominated by groundwater contributions to discharge. Also, based on the hydropedological and climatic similarities between our páramo site and the peatland-podzol dominated ecosystems in the Scottish highlands (e.g., Soulsby et al., 2006; Tetzlaff et al., 2014), we hypothesize the proportion of wetlands to be a dominant control on the variability of the MTT in this high-elevation tropical ecosystem.

Materials and methods

Study site

The Zhurucay River Ecohydrological Observatory is a basin located within a tropical alpine biome, locally known as wet Andean páramo. It is situated in southern Ecuador (3${}^{\circ }$04${}^{\prime }$ S, 79${}^{\circ }$14${}^{\prime }$ W) on the western slope of the Atlantic–Pacific continental divide and discharges into the Jubones River (Pacific Ocean tributary). The basin has a drainage area of 7.53 km${}^{\mathrm{2}}$ and extends within an elevation range of 3400 to 3900 m a.s.l. Climate is mainly influenced by the Pacific Ocean regime and to a lesser degree by the continental air masses from the Amazon basin. Mean annual precipitation at the observatory is 1345 mm at 3780 m a.s.l. Precipitation shows low seasonality with two relatively drier months (August and September) and primarily falls as drizzle (Padrón et al., 2015). Mean annual temperature is 6.0 ${}^{\circ }$C at 3780 m a.s.l. and 9.2 ${}^{\circ }$C at 3320 m a.s.l. (Córdova et al., 2015).

The geology of the region is characterized by volcanic rock deposits compacted by glacial activity during the last ice age (Coltorti and Ollier, 2000). The Quimsacocha formation, composed of basaltic flows with plagioclases, feldspars, and andesitic pyroclastics, covers the northern part of the basin. The Turi formation covers the southern part of the catchment, and its lithology mainly corresponds to tuffaceous andesitic breccias, conglomerates, and horizontally stratified sands. Both formations date from the late Miocene period (Pratt et al., 1997). The geomorphology of the landscape bears the imprint of glaciated U-shaped valleys. The average slope of the basin is 17 %. The majority of the basin (72 %) has mean gradients between 0 and 20 %, although slopes of up to 40 % are also found (24 %). There is an interesting geomorphic feature on the northeastern side of the basin corresponding to a ponded wetland at a flat hilltop. As indicated by geologists from INV Metals mining company, this structure most likely resulted from the eutrophication of a lagoon due to high accumulation of volcanic material. This area is locally known as “Laguna Ciega” (“Blind Lagoon” in Spanish) and drains towards the outlet of catchment M7 (see Fig. 1). The analysis of the water stable isotopic composition of soil water and streamflow in this area indicated that the hydrologic processes of this site occur in the shallow ponded water that is directly connected to the drainage network, while deeper water stored in the soil profile has little influence for discharge generation, most likely as a result of the eutrophic condition of the wetland (Mosquera et al., 2016).

Location of the study area, and the isotopic monitoring stations in the Zhurucay observatory for streamflow ($M$) and precipitation ($P$). SW is a spring water source upstream of the outlet of catchments M3 and M4.

[Figure omitted. See PDF]

Andosols are the dominant soil type in the study site. They cover approximately 80 % of the total basin area and are mainly found in the hillslopes. Histosols (Andean wetlands) cover the remaining portion of the basin and are mainly found in flat areas where rock geomorphology allows water accumulation (Mosquera et al., 2015). These soils, formed from the accumulation of volcanic ash in flat valley bottoms and low gradient slopes, are black, humic, and acid soils rich in organic matter with low bulk density and high water storage capacity (Quichimbo et al., 2012). The organic fraction of the Histosol soils corresponds to an H horizon (median depth 76.5 cm), while in the Andosol soils it corresponds to an Ah horizon (median depth 40 cm). The mineral fraction of both soils corresponds to a C horizon (median depth of 31 cm in the Histosols and 40 cm in the Andosols). A complete description of soil properties can be found in Mosquera et al. (2015) and Quichimbo et al. (2012). Vegetation coverage is highly correlated with the soil type. Cushion plants (such as Plantago rigida, Xenophyllum humile, and Azorella spp.) grow primarily in Histosols, while tussock grass (mainly Calamagrostis sp.) (Ramsay and Oxley, 1997; Sklenar and Jorgensen, 1999) grows in Andosols. The main landscape characteristics are summarized in Table 1.

Main landscape characteristics of the monitored catchments.

	Area	Altitude	Distribution of soil			Vegetation				Topography${}^{\text{c}}$				Geology${}^{\text{e}}$			EC${}^{\text{f}}$
Catchment	(km${}^{\mathrm{2}})$	(m a.s.l.)	types${}^{\text{a}}$ (%)			cover${}^{\text{b}}$ (%)				(%)				(%)			($\mathrm{\mu }$S cm${}^{-\mathrm{1}}$)
			AN	HS	LP	TG	CP	QF	PF	AS	$L^{\text{d}}$	G	TWI	Qm	Tu	Qd
M1	0.20	3777–3900	85	13	2	85	15	0	0	14	0.9	1.9	9.6	100	0	0	35.7
M2	0.38	3770–3900	83	15	2	87	13	0	0	24	0.8	1.9	12.8	66	1	33	32.0
M3	0.38	3723–3850	80	16	3	78	18	4	0	19	1.0	2.0	9.4	59	41	0	62.4
M4	0.65	3715–3850	76	20	4	79	18	3	0	18	1.3	2.0	12.5	50	48	1	47.9
M5	1.40	3680–3900	78	20	2	78	17	0	4	20	2.5	1.9	11.8	70	1	30	37.0
M6	3.28	3676–3900	74	22	4	73	24	1	2	18	3.3	1.8	8.2	50	30	20	35.6
M7	1.22	3771–3830	37	59	4	35	65	0	0	12	0.4	1.7	10	87	0	13	15.3
M8	7.53	3505–3900	72	24	5	71	24	2	2	17	4.6	1.9	16.7	56	31	13	33.5

${}^{\text{a}}$ AN: Andosol; HS: Histosol; LP: Leptosol. ${}^{\text{b}}$ TG: tussock grasses; CP: cushion plants; QF: Polylepys forest; PF: pine forest. ${}^{\text{c}}$ AS: average slope; L: flow path length; G: flow path gradient; TWI: topographic wetness index (Beven and Kirkby, 1979a). ${}^{\text{d}}$ $L$ units in kilometers. ${}^{\text{e}}$ Qm: Quimsacocha formation; Tu: Turi formation; Qd: Quaternary deposits. ${}^{\text{f}}$ EC: mean electrical conductivity. Data collected weekly for a 3-year period (June 2012–June 2015).

Hydrometric information

Discharge and precipitation were continuously monitored from October 2010. A nested monitoring network was used to measure discharge. The network consisted of seven tributary catchments (M1 to M7) draining to the outlet of the basin (M8). Catchments M1 to M6 comprise the main stream network draining towards the outlet of the Zhurucay basin (M8), whereas catchment M7 is a small catchment originating in a ponded wetland at a flat hilltop (Fig. 1). V-notch weirs were constructed to measure discharge at the outlet of the tributaries M1–M7 and a rectangular weir at the outlet of basin M8. Each catchment was instrumented with pressure transducers (Schlumberger DI500) with a precision of $\pm$5 mm. Water levels were recorded at a 5 min resolution and transformed into discharge using the Kindsvater–Shen relationship (US Bureau of Reclamation, 2001). The discharge equations were calibrated by applying the constant rate salt dissolution technique (Moore, 2004). Precipitation was recorded using tipping buckets with a resolution of 0.2 mm at two stations located at 3780 and 3700 m a.s.l. (Fig. 1).

Collection and analysis of water stable isotopic and electrical conductivity data

We used a 3-year record (May 2011 to May 2014) of ${}^{18}$O and ${}^{\mathrm{2}}$H isotopic compositions of water samples collected in precipitation and streamflow. Data were collected at different resolutions, from event-based to biweekly, given logistic constraints and opportunities. Higher-resolution data were aggregated to biweekly using precipitation weighted means for record consistency. The same nested monitoring network used for measuring discharge was implemented for stable isotopes in streamflow at the seven tributary catchments M1 to M7, and including M8 at the outlet of the basin. Water samples in precipitation were collected using two rain collectors located at 3780 and 3700 m a.s.l. Each collector consisted of a circular funnel and a polypropylene bottle covered with aluminum foil. Evaporation was prevented by placing a plastic sphere (4 cm diameter) in the funnel and a layer of 0.5 cm mineral oil within the polypropylene bottle. Due to the sampling procedure and the local climate, kinetic fractionation by evaporation can be neglected and hence both stable isotopes yield the same results (Mosquera et al., 2016). Therefore only the results using the isotopic composition of ${}^{18}$O are reported. Rainwater samples are cumulative representations of the isotopic signature between sampling dates, while stream grab water samples represent discrete points in time. The collected water samples were stored in 2 mL amber glass bottles, covered with parafilm, and kept away from the sunlight to prevent fractionation by evaporation, as recommended by the International Atomic Energy Agency (Mook, 2000). The isotopic composition of the water samples was measured using a cavity ring-down spectrometer – L1102-i (Picarro, USA) – with a 0.5 ‰ precision for deuterium (${}^{\mathrm{2}}$H) and 0.1 ‰ precision for oxygen-18 (${}^{18}$O). Isotopic concentrations are presented in the $\mathit{\delta }$ notation and expressed in per mill ( ‰) according to the Vienna Standard Mean Ocean Water (V-SMOW) (Craig, 1961).

Electrical conductivity (EC) was measured directly in-stream simultaneously with the water isotopic data starting in 2012, the second year of the monitoring period. EC was measured using digital conductivity sensor Tetracon 925 (WTW, Germany) with a precision of $\pm$0.5 %.

Mean transit time modeling and transit time distributions

Given the homogeneous landscape and hydrometeorological conditions in the Zhurucay basin, we estimated mean transit times (MTTs) using an inverse solution of the LCA (Amin and Campana, 1996; Małoszewski and Zuber, 1982), which assumes steady-state conditions (i.e., baseflow MTTs). The LCA seeks the parameter set of the model that best describes the hydrologic system represented by a predefined transit time distribution (TTD) function (Małoszewski and Zuber, 1996). The TTD describes the transition of an input signal (e.g., precipitation, snow) of a tracer (e.g., ${\mathit{\delta }}^{18}$O, ${\mathit{\delta }}^{\mathrm{2}}$H) to the signal at an outlet point (e.g., groundwater, streamflow) resulting from the subsurface transport of water molecules within a catchment. Mathematically the TTD is described by a convolution integral that transforms the input signal (${\mathit{\delta }}_{\mathrm{in}})$ into an output signal (${\mathit{\delta }}_{\mathrm{out}})$, considering a time lag between them $\left(t,-,\mathit{\tau }\right)$ through a transfer function $\left(\text{TTDs}\text{or}g\left(\mathit{\tau }\right)\right)$ describing the subsurface transport of tracer as follows: $${\mathit{\delta }}_{\text{out}}\left(t\right)={\int }_{\mathrm{0}}^{\mathrm{\infty }}g\left(\mathit{\tau }\right){\mathit{\delta }}_{\text{in}}\left(t,-,\mathit{\tau }\right)\text{d}\mathit{\tau },$$ where $\mathit{\tau }$ is the integration variable representing the MTT of the tracer. A more robust approximation weights the isotopic concentration of the input by considering recharge mass variation ($w\left(\mathit{\tau }\right))$ so that the outflow composition reflects the mass flux leaving the catchment: $${\mathit{\delta }}_{\text{out}}\left(t\right)={\displaystyle \frac{{\int }_{\mathrm{0}}^{\mathrm{\infty }}g\left(\mathit{\tau }\right)w\left(t,-,\mathit{\tau }\right){\mathit{\delta }}_{\text{in}}\left(t,-,\mathit{\tau }\right)\text{d}\mathit{\tau }}{{\int }_{\mathrm{0}}^{\mathrm{\infty }}g\left(\mathit{\tau }\right)w\left(t,-,\mathit{\tau }\right)\text{d}\mathit{\tau }}},$$ where $\mathrm{w}\left(t,-,\mathit{\tau }\right)$ can be described in terms of rainfall magnitude, intensity, or effective precipitation (McGuire and McDonnell, 2006). Precipitation intensity was used to volume weight the isotopic composition of precipitation in our study. Recharge was represented by the rainfall isotopic composition weighted by precipitation rate and accounted for relatively small recharge (i.e., lower precipitation inputs) during the less wet months (August and September).

MTT was estimated by adjusting the response function or TTD to fit the measured and simulated stream-water isotopic composition. Five TTDs were considered to describe the subsurface transport of water molecules in the Zhurucay basin. We used the exponential model (EM), the exponential-piston flow model (EPM), the dispersion model (DM) (Małoszewski and Zuber, 1982), the gamma model (GM) (Kirchner et al., 2000), and the two parallel linear reservoir model (TPLR) (Weiler et al., 2003). Each model is briefly described below and Table 2 summarizes their equations, fitting parameters, and the range of initial parameters used in this study.

Models considered to describe water mean transit time (MTT) in the study area and their transit time distribution (TTD) functions, parameters, and range of initial parameters.

Model	Transit time distribution (g($\mathit{\tau }\left)\right)$	Parameter(s) range
Exponential model (EM)	$\frac{\mathrm{1}}{\mathit{\tau }}\exp \left(\frac{-t}{\mathit{\tau }}\right)$	$\mathit{\tau }[\mathrm{0}-65]$
Exponential-piston model (EPM)	$\frac{\mathit{\eta }}{\mathit{\tau }}\exp \left(-,\frac{t\cdot \mathit{\eta }}{\mathit{\tau }},+,\mathit{\eta },-,\mathrm{1}\right)\mathrm{for}t\ge \mathit{\tau }(\mathrm{1}-{\mathit{\eta }}^{-\mathrm{1}}$	$\mathit{\tau }[\mathrm{0}-65]$
		$\mathit{\eta }[0.5-\mathrm{4}]$
Dispersion model (DM)	${\left(\frac{\mathrm{4}\mathit{\pi }{D}_{\text{p}}t}{\mathit{\tau }}\right)}^{-\mathrm{1}/\mathrm{2}}{t}^{-\mathrm{1}}\exp \left[-,{\left(\mathrm{1},-,\frac{t}{\mathit{\tau }}\right)}^{\mathrm{2}},\left(\frac{\mathit{\tau }}{\mathrm{4}{D}_{\text{p}}t}\right)\right]$	$\mathit{\tau }[\mathrm{0}-65]$
		$D_{\text{p}}[0.5-\mathrm{4}]$
Gamma model (GM)	$\frac{{\mathit{\tau }}^{\mathit{\alpha }-\mathrm{1}}}{{\mathit{\beta }}^{\mathit{\alpha }}\mathrm{\Gamma }\left(\mathit{\alpha }\right)}{\exp }^{-\mathit{\tau }/\mathit{\beta }}$	$\mathit{\alpha }\left[0.01,-,\mathrm{4}\right]$
		$\mathit{\beta }[\mathrm{0}-500]$
		$\mathit{\tau }=\mathit{\alpha }\cdot \mathit{\beta }$
Two parallel linear reservoir (TPLR)	$\frac{\mathit{\phi }}{{\mathit{\tau }}_{\text{f}}}\exp \left(\frac{-t}{{\mathit{\tau }}_{\text{f}}}\right)+\frac{\mathrm{1}-\mathit{\phi }}{{\mathit{\tau }}_{\text{f}}}\exp \left(\frac{-t}{{\mathit{\tau }}_{\text{s}}}\right)$	${\mathit{\tau }}_{\text{s}}[\mathrm{0}-65]$
		${\mathit{\tau }}_{\text{f}}[\mathrm{0}-\mathrm{8}]$
		$\mathit{\phi }[\mathrm{0}-\mathrm{1}]$

$\mathit{\tau }$: tracer's mean transit time (MTT) (months); $\mathit{\eta }$: parameter that indicates the percentage of contribution of each flow type ($-$); $D_{\text{p}}$: dispersion parameter ($-$); $\mathit{\alpha }$: shape parameter ($-$); $\mathit{\beta }$: scale parameter ($-$); ${\mathit{\tau }}_{\text{f}}$ and ${\mathit{\tau }}_{\text{s}}$: transit time of fast and slow flows in months; $\mathit{\phi }$: flow partition parameter between fast and slow flow reservoirs (%). ($-$) indicates that parameters are unitless.

The EM represents a well-mixed system and assumes contributions from all flow paths. It assumes a relatively simple transition of the tracer towards the stream network. The EPM is an extension of the EM in which a delay in the shortest flow paths is assumed by the piston flow portion of the system. In addition to the MTT, it has an additional fitting parameter ($\mathit{\eta }$), which represents the ratio of the total volume to the volume represented by the exponential distribution. The DM arises from the solution of the one-dimensional advection–dispersion equation (Kreft and Zuber, 1978) and assumes that there is influence of hydrodynamic dispersion in the system's flow paths. It also has two fitting parameters, the MTT and the dispersion parameter ($D_{\text{p}}$), which relates to the tracer transport process. The GM is a more flexible and general version of the exponential model in which the product of two parameters provides an estimation of the MTT of the system. These parameters are the shape parameter ($\mathit{\alpha }$) and the scale parameter ($\mathit{\beta }$) (Kirchner et al., 2000). The TPLR represents two parallel reservoirs, each one represented by a single exponential distribution. It has three fitting parameters, the MTT of the slow (MTT${}_{\text{s}})$ and fast (MTT${}_{\text{f}})$ reservoirs and a parameter representing the fraction of each of them with respect to total flow ($\mathit{\phi })$ (Weiler et al., 2003).

The MTT approach is based on the following assumptions: (1) the tracer concentration is conservative (i.e., the tracer does not react with other elements present in the system); (2) the input tracer concentration is input in flux mode (i.e., volumetrically weighted); (3) the tracer enters the system only once and uniformly; (4) a representative tracer input can be identified; (5) transport of solute is one-dimensional and represented by a single TTD; and (6) there is a uniform storage of water within the catchment (i.e., steady-state flow in the system) (Małoszewski and Zuber, 1982). The steady-state assumption is assumed to be accomplished in humid environments during baseflow conditions (McGuire et al., 2002). In order to comply with the latter assumption, streamflow water samples collected during extreme rainfall events were excluded for the MTT simulations (McGuire and McDonnell, 2006; Muñoz-Villers et al., 2016). As a result, estimates correspond to baseflow MTTs, hereafter simply referred to as MTT. To obtain more stable results, we looped the available 3 years of isotopic data ten times during calibration in order to extend the data series for 30 years as a warm-up period following Hrachowitz et al. (2011) and Timbe et al. (2014).

Model performance and uncertainty analysis

The model performance was evaluated using the Kling–Gupta efficiency coefficient (KGE) (Gupta et al., 2009). KGE ranges from $-\mathrm{\infty }$ to 1, where unity indicates an ideal optimization. KGE can be viewed from a multi-objective perspective because it accounts for correlation (i.e., balancing dynamics, $r$), variability error ($\mathit{\gamma }$), and bias error ($\mathit{\beta }$) within a single objective function. The efficiency is mathematically represented by the Euclidean distance (ED) in each of the three dimensions ($r$, $\mathit{\gamma }$, and $\mathit{\beta })$ to an ideal point where all of them are maximized (i.e., where ideally the three factors are set to one). Efficiencies lower than 0.45 were considered poor predictions (Timbe et al., 2014).

Depending on the TTD function used, one to three parameters were fitted during the simulations. Models were built using parameter sets generated through a uniform Monte Carlo sampling procedure (Beven and Freer, 2001). Each model was first run 10 000 times within a wide range of parameter values (Table 2). Once a parameter value that yielded the best KGE was clearly identified, the model was run again within a narrowed range of parameters until at least 1000 behavioral solutions were obtained (i.e., solutions corresponding to at least 95 % of the highest KGE) (Timbe et al., 2014), and their 5 and 95 % limit bounds (i.e., 90 % confidence interval) were estimated using the generalized likelihood uncertainty estimation (GLUE) methodology (Beven and Binley, 1992). The Akaike information criterion (AIC, Akaike, 1974) was used as a parsimoniousness metric for model selection that penalizes model performance based on the number of fitted parameters used to calibrate each model. The model with the lowest AIC is the most efficient at fitting the observed values. A visual inspection of the parameter space (Fig. 4) was also conducted to qualitatively analyze the models' parameter identification. The best model describing the hydrologic conditions of the system was selected using the following criteria: (1) best goodness of fit using the KGE criterion, (2) lowest AIC, (3) results that yielded the lower uncertainty estimations, and (4) strong parameter identification.

Correlation analysis of MTT and catchment characteristics

We used linear regression to investigate relations between landscape characteristics and hydrological behavior with the MTT of the catchments. For this analysis, we included the catchments that comprise the main drainage network (i.e., catchments M1 to M6) and the catchment outlet (M8) given that they possess comparable hydropedological and geomorphological characteristics. That is, catchments situated at the valley bottom have well-defined interconnections between wetlands in the riparian areas and the surrounding Andosol soils at the slopes. Catchment M7, on the other hand, is located at a flat hilltop at the outlet of a wetland area that remains ponded throughout the year. The geomorphology of this concave (lagoon-shaped) structure and its ponded eutrophic condition has allowed for the hydrologic processes to mainly occur in the shallowest ponded portion of the water directly connected to the stream network (with little influence of the most likely immobile water, which remains stored in the deeper soil fraction) (Mosquera et al., 2016). Therefore, its hydrological response is not comparable to the other catchments, where hydrologic processes mainly occur in the soils, and consequently it was excluded from the regression analysis. The statistical significance of the correlations was tested using the $F$ test at a 95 % confidence level (i.e., $p$ < 0.05).

The landscape and hydrometric variables tested for correlation were obtained from previous studies at the site (Mosquera et al., 2015) and from detailed soil, vegetation, and topographic information provided by INV Metals. The landscape features considered were soil type, vegetation, geology, catchment size, slope, flow path length and gradient, and topographic wetness index (TWI) (Beven and Kirkby, 1979b) (Table 1). The hydrometric variables considered were annual runoff, annual precipitation, runoff coefficient, and streamflow rates (Table 3). The weekly collected EC for 3 years (June 2012 to June 2015) was averaged and also tested for correlation with MTT.

Main hydrometric variables of the catchments.

	Precipitation	Total runoff	Runoff	Average specific	Flow rates as frequency of non-exceedance (L s${}^{-\mathrm{1}}$ km${}^{-\mathrm{2}})$
Catchment	(mm yr${}^{-\mathrm{1}})$	(mm yr${}^{-\mathrm{1}})$	coefficient${}^{\text{*}}$	discharge (L s${}^{-\mathrm{1}}$ km${}^{-\mathrm{2}})$
					$Q_{\text{min}}$	$Q_{10}$	$Q_{30}$	$Q_{50}$	$Q_{70}$	$Q_{90}$	$Q_{\text{max}}$
M1	1300	729	0.56	23.1	0.7	2.7	6.6	14.3	26.4	50.1	1039.0
M2	1300	720	0.55	22.8	1.2	4.8	7.9	14.9	26.7	49.0	762.9
M3	1293	841	0.65	26.7	2.3	7.3	10.8	17.7	28.1	52.4	894.2
M4	1294	809	0.62	25.6	4.2	6.2	9.8	16.6	27.3	52.1	741.2
M5	1267	766	0.6	24.3	1.5	4.1	8.3	15.3	26.9	50.8	905.7
M6	1254	786	0.63	24.9	1.2	3.7	8.2	15.9	27.5	53.2	930.4
M7	1231	684	0.56	21.7	0.3	1.8	5.2	11.0	23.3	53.9	732.0
M8	1277	864	0.68	27.4	1.9	4.0	8.7	15.2	29.2	60.8	777.9

${}^{\text{*}}$ Total runoff as a proportion of precipitation.

Results

Hydrologic and isotopic characterization in rainfall and streamflow

Precipitation in the Zhurucay basin is evenly distributed throughout the year (Fig. 2a), except for 2 months with relatively lower precipitation inputs (i.e., August and September), both accounting for less than 8 % of total annual precipitation. Spatially, annual precipitation ($P$) is evenly distributed across the basin, with an average of 1275 $\pm$ 9 mm. Total annual runoff ($Q$) and runoff coefficients ($Q$ $/$ $P$) are spatially more heterogeneous, varying between 684 and 864 mm per year and 0.55 and 0.68, respectively (Table 3). The hydrograph at the outlet of the basin (M8) also depicts a flashy response to precipitation inputs, even during these less humid months (see zoom in Fig. 2a). Similar behavior is observed at all catchments.

(a) Hourly precipitation and unit area streamflow; (b) ${\mathit{\delta }}^{18}$O isotopic composition in precipitation and streamflow for 3 years (May 2011 to May 2014); and (c) electrical conductivity for 2 years (May 2012 to May 2014) at the catchment outlet (M8; see location in Fig. 2). The size of the bubbles in plot (b) indicates the relative cumulative rainfall in millimeters for each collected sample.

[Figure omitted. See PDF]

The ${\mathit{\delta }}^{18}$O isotopic composition in rainfall is highly variable throughout the year (e.g., average $-$10.2 $\pm$ 0.32‰ at the upper station) (Fig. 2b) and follows a seasonal pattern with isotopically enriched values during highest precipitation rates (April–May) and isotopically depleted values in the less humid period (August–September). The ${\mathit{\delta }}^{18}$O isotopic composition in streamflow collected during low flows on the other hand, is much more damped (average $-$10.0 $\pm$ 0.06 ‰, at M8) than the isotopic composition in precipitation (Table 4).

Statistics of the ${\mathit{\delta }}^{18}$O isotopic composition in precipitation and streamflow used as input data for the MTT modeling.

Sampling	Altitude	${\mathit{\delta }}^{18}$O (‰)
station	(m a.s.l.)	$n^{\text{a}}$	Average	SE${}^{\text{b}}$	Max	Min
M1	3840	123	$-$10.6	0.06	$-$9.0	$-$12.6
M2	3840	124	$-$10.4	0.07	$-$8.8	$-$12.6
M3	3800	121	$-$10.7	0.05	$-$8.8	$-$12.1
M4	3800	122	$-$10.6	0.05	$-$8.7	$-$11.9
M5	3800	118	$-$10.5	0.06	$-$9.1	$-$12.8
M6	3780	121	$-$10.3	0.06	$-$8.9	$-$12.2
M7	3820	121	$-$8.9	0.15	$-$6.2	$-$13.9
M8	3700	118	$-$10.0	0.06	$-$8.3	$-$11.6
Upper precip.	3779	137	$-$10.2	0.32	$-$1.2	$-$25.0
Middle precip.	3700	134	$-$10.1	0.32	$-$2.7	$-$20.0

${}^{\text{a}}$ $n$: number of samples collected. ${}^{\text{b}}$ SE: standard error.

TTD evaluation and selection

In order to identify the TTD that best describes the hydrologic system in the Zhurucay basin, we tested and evaluated the performance of all TTDs at all catchments. Considering that similar results were obtained for all catchments and for brevity, only results for M8, the basin outlet, are shown.

All TTDs reproduce the ${\mathit{\delta }}^{18}$O isotopic composition at the outlet of the basin (M8), with efficiencies varying between 0.50 and 0.76, i.e., above the threshold of model acceptance (KGE > 0.45) (Table 5). The more flexible models, GM and TPLR, yield the highest performances, with KGEs of 0.75 and 0.76, respectively. The EM and the EPM yield similar efficiencies (KGE $=$ 0.63), while the DM yields the lowest efficiency among all (KGE $=$ 0.50). The models associated with the highest KGEs yield the highest uncertainty bounds according to their threshold of behavioral solutions, whereas the EM shows the lowest uncertainty (Fig. 3). The lowest AIC value was determined for the EM (AIC $=$ 2.92), indicating that this model is the most parsimonious for the system under investigation. Higher AIC values were determined for the other models that fit more than one parameter (GM, 4.58; EPM, 4.92; DM, 5.39; TPLR, 6.55). The models' visual parameter identification inspection indicates that even though the TPLR model yields the highest KGE, the identification of its parameters is poor, particularly the one representing the MTT of the slow reservoir (MTT${}_{\text{s}}$, Fig. 4). For the DM and GM models one parameter is well identified (MTT and $\mathit{\alpha }$, respectively), while the others are not well identified by the models. For the EPM, both parameters are not well identified. The single parameter that defines the EM is very well identified. This coupled analysis of model efficiency and parsimoniousness in combination with parameter identifiability indicates that although models with a higher number of fitting parameters provide higher efficiencies, their parameters are more uncertain.

Statistical parameters of observed and modeled ${\mathit{\delta }}^{18}$O for the stream at the outlet of the basin (M8).

Model	Observed ${\mathit{\delta }}^{18}$O		Simulated ${\mathit{\delta }}^{18}$O						Model parameters${}^{\text{b}}$
	Mean	${\mathit{\sigma }}^{\text{a}}$	Mean	${\mathit{\sigma }}^{\text{a}}$	KGE${}^{\text{a}}$	RMSE${}^{\text{a}}$	Bias	AIC${}^{\text{a}}$
	(‰)	(‰)	(‰)	(‰)	($-$)	(‰)	(‰)	($-$)
EM	$-$10.05	0.45	$-$10.02	0.52	0.63	0.59	0.03	2.92	$\mathit{\tau }$ (days)	191 (166–224)
EPM	$-$10.05	0.45	$-$10.02	0.52	0.63	0.59	0.03	4.92	$\mathit{\tau }$ (days)	120 (90–205)
									$\mathit{\eta }$ ($-$)	0.63 (0.53–0.92)
DM	$-$10.05	0.45	$-$9.95	0.64	0.50	0.72	0.10	5.39	$\mathit{\tau }$ (days)	666 (480–755)
									$D_{\text{p}}$ ($-$)	3.99 (2.74–3.93)
GM	$-$10.05	0.45	$-$10.02	0.45	0.75	0.47	0.03	4.58	$\mathit{\tau }$ (days)	392 (296–1478)
									$\mathit{\alpha }$ ($-$)	0.70 (0.51–0.76)
									$\mathit{\beta }$ ($-$)	40 (26–108)
TPLR	$-$10.05	0.45	$-$10.02	0.45	0.76	0.47	0.03	6.55	${\mathit{\tau }}_{\text{f}}$ (days)	1664 (356–1738)
									${\mathit{\tau }}_{\text{s}}$ (days)	45 (12–90)
									$\mathit{\phi }$ ($-$)	0.31 (0.10–0.48)

${}^{\text{a}}\mathit{\sigma }$: standard deviation; KGE: Kling–Gupta efficiency (Gupta et al., 2009); AIC: Akaike information criterion (Akaike, 1974). Statistical parameters of the simulated results correspond to the best-matching value of the objective function KGE. ${}^{\text{b}}\mathit{\tau }$: tracer's mean transit time; $\mathit{\eta }$: parameter that indicates the ratio between the contribution of piston and exponential flow; $D_{\text{p}}$: dispersion parameter; $\mathit{\alpha }$: shape parameter; $\mathit{\beta }$: scale parameter; ${\mathit{\tau }}_{\text{f}}$ and ${\mathit{\tau }}_{\text{s}}$: transit time of fast and slow flows; $\mathit{\phi }$: flow partition parameter between fast and slow flow reservoirs. ($-$) : dimensionless parameter. Uncertainty bounds (5–95th percentiles) of simulated parameters shown in parentheses were estimated using the generalized likelihood uncertainty estimation (GLUE, Beven and Binley, 1992).

Fitted results of the five lumped parameter models used to simulate the temporal variability in the ${\mathit{\delta }}^{18}$O streamflow composition at the outlet of the basin (M8). (a) Exponential model (EM); (b) exponential-piston model (EPM); (c) dispersion model (DM); (d) gamma model (GM); and (e) two parallel linear reservoir model (TPLR). The open circles represent the observed isotopic composition in streamflow; the black line represents the best simulated isotopic composition in streamflow according to the KGE (Gupta et al., 2009) objective function; and the blue shaded area corresponds to the 5–95 % confidence limits of the possible solutions from the parameter sets within the range of behavioral solutions, i.e., solutions that yield at least 95 % KGE.

[Figure omitted. See PDF]

Taking into account the goodness of fit, parsimoniousness, uncertainty bounds, and the identifiability of the models' parameters, the EM is the model that best describes the temporal variability of the baseflow ${\mathit{\delta }}^{18}$O isotopic composition across the Zhurucay basin. These results are supported by previous findings of Mosquera et al. (2015, 2016) at the study site, indicating a fast response system with minimal impact of old water sources, and the further discussion in the final part of this paper. Although the contribution of old water component/s cannot be definitely neglected, it seems unfeasible that given the biophysical features of the landscape and its pedological and geomorphological structure, baseflow MTTs in the range of those yielded by the GM and TPLR TTDs are reliable. The EM represents a well-mixed reservoir with relatively simple transition of the water (i.e., tracer) in the subsurface towards the stream network. In the Zhurucay basin, the organic and porous páramo soils allow for the efficient mixing of water within the whole profile of these poorly developed soils. This effect particularly occurs in the Histosols (Andean wetlands) that are directly connected, hydrologically, to the drainage network (Mosquera et al., 2016). These factors result in a relatively simple transition of infiltrated precipitation from a well-mixed reservoir towards the catchment outlet. This process-based analysis of physical characteristics of the system further supports the EM as the model that best describes the transport of water across the Zhurucay basin. The EM was also found to describe the subsurface transport of water in another system of catchments in eastern Mexico, where soils have predominantly formed by volcanic ash accumulation (Muñoz-Villers et al., 2016).

Baseflow MTT

Results of the EM for selected catchments with the longest (M3), intermediate (M6), and lowest (M7) MTTs are shown in Fig. 5, and statistics of the EM simulations at all catchments are summarized in Table 6. The EM overcomes the modeling acceptance criterion of KGE > 0.45 at all catchments, with KGE values ranging between 0.48 and 0.84. The longest MTT is found in catchment M3 (0.73 years, 264 days), and the shortest at M7 (0.15 years, 53 days). The MTT for the other catchments varies within this range. On average, within the 90 % confidence level for the catchments forming the main drainage network (M1–M6 and M8), MTT estimations show small variations (25 days at the lower confidence bound and 35 days at the upper confidence bound) with small standard deviation (4 days for the upper bound and 6 days for the lower bound). For catchment M7, variations are even smaller (9 days at the lower confidence bound and 11 days at the upper confidence bound). In addition, the model performs best for catchments with high variability in their isotopic composition during the monitoring period. For instance, catchment M3 (Fig. 5a) shows the smallest amplitude in isotopic variation for the observed and simulated data (Table 6), coupled with the lowest KGE (0.48) and the highest MTT. On the other hand, catchment M7 (Fig. 5c) shows the highest amplitude in isotopic variation for the observed and simulated data, coupled with the highest KGE (0.84) and the shortest MTT. Similarly, catchment M6 (Fig. 5b), which has a MTT shorter than the one in M3 and longer than the one in M7, has an amplitude and KGE varying between the ones in M3 and M7. The Monte Carlo simulations for the fitted parameter MTT (Fig. 5) clearly depict how the MTTs that yield the highest KGE in each catchment decrease as the variation in their isotopic composition increases as described above. Results from all the catchments are also described by this trend.

Monte Carlo simulations of the fitted parameters of the five lumped parameter models used to simulate the ${\mathit{\delta }}^{18}$O streamflow composition at the outlet of the basin (M8). (a) EM; (b) EPM; (c) DM; (d) GM; and (e) TPLR. The ($-$) symbol on the $x$ axes denotes that the fitting parameter is dimensionless. Horizontal red lines indicate the threshold of behavioral solutions (at least 0.95 of maximum KGE).

[Figure omitted. See PDF]

Fitted results and Monte Carlo simulations of the fitted parameters of the exponential model (EM) used to simulate the ${\mathit{\delta }}^{18}$O streamflow composition in the catchments: (a) M3; (b) M6; and (c) M7. The open circles represent the observed isotopic composition in streamflow; the red crosses represent the isotopic composition in precipitation; the black line represents the best simulated isotopic composition in streamflow according to the KGE objective function; and the blue shaded area corresponds to the 5–95 % confidence limits of the possible solutions from the MTT fitting parameters within the range of behavioral solutions, i.e., solutions that yield at least 95 % KGE. Panels on the right represent the explored parameter range for the MTT parameter and the KGEs associated with each of them.

[Figure omitted. See PDF]

The MTT probability density functions (PDFs), which indicate the distribution of MTTs in the hydrologic system, and cumulative density functions (CDFs), which express the tracer “mass recovery from an instantaneous, uniform tracer mass addition” (McGuire et al., 2005) based on the fitted parameter distributions, not shown for brevity, exhibit a dominance of relatively short MTTs in the hydrology of the Zhurucay basin. The CDFs also indicate that the tracer is completely recovered ($\sim$ 100 %) in all catchments at around 160 weeks, except for M7, where the tracers are even more rapidly recovered ($\sim$ 38 weeks). As we used a stable isotopic record of 156 weeks (3 years), these results indicate that a 3-year record of tracer data is sufficient to estimate the MTT of waters using the LCA in this páramo basin.

Correlations of baseflow MTT with landscape and hydrometric variables

Correlation analysis showed no statistically significant correlations ($p$ values > 0.05) between MTT and landscape features and hydrometric variables of the nested monitoring system when all catchments were included. This lack of correlation is likely related to the previously reported distinct responsiveness of catchment M7 to precipitation inputs due to its different geomorphologic configuration (i.e., ponded eutrophic wetland disconnected from the slopes) and catchments M3 and M4 (Fig. 1) driven by a spring water contribution during low flow generation. In general, deep subsurface flow and groundwater contributions to discharge seem to be minimal, and geology has not been found to directly control the hydrology in this páramo ecosystem (Mosquera et al., 2015). However, the existence of this shallow spring sourced at the interface between the soil mineral horizon and the shallow bedrock upstream of the outlet of M3 and M4 favors the generation of higher low flows (Mosquera et al., 2016) and increases MTTs in these catchments. The latter indicates that geology (fractures in the shallow bedrock) influences the hydrology of these small headwater catchments, thus masking relationships between landscape features and the MTT of the whole system. Therefore, we tested the MTT correlations without including these small catchments (M3 and M4) and M7.

Statistical parameters of observed and simulated ${\mathit{\delta }}^{18}$O for all catchments using the exponential model (EM).

	Observed ${\mathit{\delta }}^{18}$O		Simulated ${\mathit{\delta }}^{18}$O			Simulated MTT${}^{\text{b}}$
Catchment	Mean	${\mathit{\sigma }}^{\text{a}}$	Mean	${\mathit{\sigma }}^{\text{a}}$	KGE${}^{\text{a}}$
	(‰)	(‰)	(‰)	(‰)	($-$)	(Years)	(Days)
M1	$-$10.63	0.37	$-$10.51	0.44	0.48	0.54 (0.48–0.63)	194 (171–227)
M2	$-$10.46	0.53	$-$10.51	0.63	0.61	0.43 (0.38–0.51)	156 (137–183)
M3	$-$10.64	0.23	$-$10.51	0.26	0.48	0.73 (0.64–0.86)	264 (232–310)
M4	$-$10.63	0.27	$-$10.51	0.3	0.48	0.67 (0.59–0.78)	240 (212–280)
M5	$-$10.51	0.39	$-$10.51	0.46	0.53	0.52 (0.46–0.61)	188 (165–219)
M6	$-$10.37	0.42	$-$10.43	0.5	0.59	0.52 (0.46–0.61)	188 (164–220)
M7	$-$8.93	2.92	$-$10.02	2.93	0.84	0.15 (0.12–0.18)	53 (45–64)
M8	$-$10.05	0.45	$-$10.02	0.52	0.63	0.53 (0.46–0.62)	191 (167–224)

${}^{\text{a}}\mathit{\sigma }$: standard deviation; KGE: Kling–Gupta efficiency. Statistical parameters of the simulated results correspond to the best-matching value of the objective function KGE. ${}^{\text{b}}$ Uncertainty bounds (5–95th percentiles) of the simulated mean transit time (MTT) shown in parentheses were estimated using the generalized likelihood uncertainty estimation (GLUE).

The reanalysis with the modified data set revealed significant relations of MTT to topographical indexes (Fig. 6). The relations between MTT and average slope (Fig. 6a, $R^{\mathrm{2}}=$ 0.78, $p=$ 0.047) and percent area having slopes in the range 20–40 % are negative (Fig. 6b, $R^{\mathrm{2}}=$ 0.90, $p=$ 0.015). Conversely, the relation between the percent area having slopes 0–20 % and MTT is positive ($R^{\mathrm{2}}=$ 0.85, $p=$ 0.026). No significant correlations ($p$ values > 0.05) between MTT and vegetation, soil types, geology, flow path length, topographic wetness index, and hydrometric variables were found (Table 7). However, a relatively strong relation between MTTs and low flows with smaller significance was also found. That is, catchments with higher MTTs yielded lower low flows ($R^{\mathrm{2}}=$ 0.62, $p=$ 0.11 for $Q_{\mathrm{5}}$ and $R^{\mathrm{2}}=$ 0.61, $p=$ 0.12 for $Q_{10})$. The regression analysis including all catchments also showed that mean electrical conductivity (MEC) of the waters explains 89 % ($p=$ < 0.001) of the catchments' MTT variability (Fig. 7). Streams with higher MEC have longer MTTs.

Correlations between mean transit time (MTT) and topographic indexes of the catchments: (a) average catchment slope; (b) catchment area with slopes between 0 and 20 %; and (c) catchment area with slopes between 20 and 40 %. Catchments M3 and M4 (additional spring water source; see Fig. 1) and M7 (at a flat hilltop disconnected from the hillslopes) are not included in the regressions. Solid lines are linear regressions and dashed lines are the 90 % confidence intervals of the regressions. ${}^{*}$ indicates variables are normalized by their means.

[Figure omitted. See PDF]

Discussion

General hydrometric and isotopic characterization

The rainfall–runoff process evidences a rapid response of discharge to precipitation inputs in the Zhurucay basin. This rapid response occurs even during the less humid periods (August–September) in which relatively small rainfall events result in peak flow generation (Fig. 2a). This high responsiveness results from the combined effect of the relatively uniform distribution of precipitation year-round – common in tropical regions – and the unique properties of the Histosol soils or Andean wetlands located near the streams. The high storage capacity of wetlands was also highlighted by Roa-García and Weiler (2010) after the comparison of three paired catchments in the growing coffee region of Colombia at lower elevations (2000–2200 m a.s.l.). Similarly, Histosol soils in our study site are rich in organic matter content (mean 86 % by volume), allowing for high water storage capacity. In addition, due to their relatively low saturated hydraulic conductivity (0.72–1.55 cm h${}^{-\mathrm{1}})$, these soils remain near saturation throughout the year. These factors, in combination with the local climate, allow páramo soils to regulate and maintain a sustained discharge throughout the year. Moreover, as these processes occur in the shallow organic horizon of the soils, the hydrology of the Zhurucay basin páramo ecosystem is dominated by shallow subsurface water flow. This is supported by the similar isotopic composition between streams and soil waters in the organic layer of the Histosols in the Zhurucay basin (Mosquera et al., 2016).

Although the ${\mathit{\delta }}^{18}$O isotopic composition of stream waters is damped and lagged with respect to that of precipitation, streamflow samples in the Zhurucay basin still reflect the variability of the ${\mathit{\delta }}^{18}$O composition of rainfall (Fig. 2b), as expected in a system dominated by shallow subsurface flow. However, catchments are differently influenced by precipitation. M7, located at the outlet of a wetland that remains constantly ponded, shows a faster response to rainfall (Fig. 5c), most likely as a result of the rapid mixing of rainfall water with the shallow water moving in the shallow organic horizon of the soils and the ponded water above it. All of the other catchments show considerably less influence of rainfall (Fig. 5a and 5b) due to the mixing of rainfall water with the water stored within the whole profile of the Histosol soils.

Coefficient of determination ($R^{\mathrm{2}})$ between the mean transit time (MTT) and (i) landscape features and (ii) hydrological variables for each of the catchments. Catchments M3 and M4 (additional spring water source; see Fig. 1) and M7 (at a flat hilltop disconnected from the hillslopes) are not included in the regressions, except for electrical conductivity; i.e., all catchments are considered (Fig. 7).

	Vegetation				General features
Landscape features	Cushion plant	0.29	(0.35)	Hydrologic variables	Runoff coefficient	0.25	(0.39)
	Tussock grass	$-$0.31	(0.33)		Total runoff	0.21	(0.44)
					Precipitation	$-$0.17	(0.48)
	Soil type				Average specific discharge	0.21	(0.44)
	Histosol	0.13	(0.56)
	Andosol	– 0.13	(0.55)		Streamflow rates
					$Q_{99}$	0.42	(0.24)
	Geologic formation				$Q_{90}$	0.18	(0.48)
	Quimsacocha	0.04	(0.75)		$Q_{80}$	0.06	(0.68)
	Turi	0.12	(0.57)		$Q_{70}$	0.09	(0.63)
	Quaternary deposits	$-$0.51	(0.18)		$Q_{60}$	0.06	(0.68)
					$Q_{50}$	0.01	(0.86)
	Topographic features				$Q_{40}$	0.10	(0.61)
	Average slope	$-$0.78	(0.05)		$Q_{30}$	$-$0.02	(0.82)
	Slope 0–20 %	0.85	(0.03)		$Q_{20}$	$-$0.14	(0.53)
	Slope 20–40 %	$-$0.90	(0.01)		$Q_{10}$	$-$0.61	(0.12)
	Area	0.13	(0.56)		$Q_{\mathrm{5}}$	$-$0.62	(0.11)
	TWI	$-$0.03	(0.79)
	Flow path length ($L$)	0.23	(0.42)		Water intrinsic properties
	Flow path gradient ($G$)	$-$0.20	(0.45)		Electrical conductivity	0.89	(< 0.001)
	$L/G$	0.23	(0.45)

Signs indicate positive (no sign) or negative ($-$) correlation between parameters. Values in bold are statistically significant to a 95 % level of confidence ($p$ < 0.05). Values in parentheses are the $p$ values of the correlations. ${}^{*}$ TWI: topographic wetness index (Beven and Kirkby, 1979).

What is the baseflow MTT?

The high performance (KGE > 0.48, Table 6), low uncertainty (Fig. 3), high AIC (Table 5), and strong parameter identification (Fig. 4) of the exponential model (EM) indicate that this model best mimics the subsurface transport of water in all catchments within the Zhurucay basin (Fig. 5). In addition, the model captures some particularities in the functioning of each catchment. For instance, results indicate relatively long MTTs in two of the headwater catchments, M3 and M4 (0.73 and 0.67 years, respectively). This likely results from a shallow spring water contribution to these catchments during low flow generation (Mosquera et al., 2015). The model seems to capture the effect of the shallow spring contribution by yielding the longest MTT estimations in these catchments, and an intrinsic influence of geology on MTT variability. In addition, the performance of the model in these two catchments is the lowest within the basin, the latter most likely because of less efficient mixing of water due to the influence of the spring water source, suggesting that this effect is also captured by the model that assumes a well-mixed reservoir. This observation led us to consider that another model representing an additional slow flow reservoir (e.g., TPLR or GM) might better represent the subsurface movement of water in these catchments. Results, however (Figs. 3 and 4, Table 5), indicate that the contribution from this additional water source is small and an additional reservoir is not well distinguished by these TTDs as their parameters are not well identified. Recently, Muñoz-Villers et al. (2016) also identified the EM as the model that best mimics subsurface flow in 7 of 12 nested catchments underlain by volcanic soils (Andosols) in a tropical montane cloud forest (TMCF) located in Veracruz, Mexico. These authors estimated even longer MTTs (1.2–2.2 years) due to deeper groundwater contributions to discharge.

On the other end, M7, dominated by the contribution from the shallowest part of the organic horizon of the soils and the ponded fraction of water accumulated in a ponded wetland – which is directly connected to the stream channel – presents the shortest MTT of all catchments (0.15 years, 53 days), linked to the highest model performance. Our results support the hypothesis that this catchment presents a shorter MTT, indicating that the ponded condition of the wetland allows for a rapid and efficient mixture of precipitated water with ponded water and water stored in and released from the shallow organic horizon of the soil. The latter results in a rapid delivery of event (new) water to the stream, whereas water stored deeper in the soil seems to remain mostly immobile, with minimal influence in the hydrology of the catchment. The efficient mixing and simpler delivery of water to the stream is also captured by the high model performance. The MTTs estimated for the rest of the catchments lie in between these two extremes and their values and efficiencies vary depending on the amplitude of the isotopic tracer variation, with longer MTTs in catchments where the amplitude of the signal is more damped – evidencing lower influence of precipitation and less efficient mixing with the soil storage – and vice versa. The MTTs in these catchments vary between 0.43 and 0.53 years (156 to 191 days). Overall, the MTTs are relatively short, further supporting previous evidence that shallow subsurface flow dominates the hydrology of the ecosystem.

In other tropical latitudes, MTTs higher than 300 days were found in three paired Colombian catchments applying the TPLR model (Roa-García and Weiler, 2010). These basins show higher MTTs than the catchments in the Zhurucay basin, most likely as a result of the higher development of the volcanic ash soils (> 10 m), which allow the water to be stored for longer periods in the subsurface. MTTs longer than 2 years were also found in a TMCF in southern Ecuador (Timbe et al., 2014), evidencing that, differently from our findings, this lower elevation ecosystem is dominated by deep groundwater contributions. Preliminary MTT estimations in another TMCF biome located in central Mexico (Muñoz-Villers and McDonnell, 2012) yielded a MTT of 3 years. Although the ecosystem is dominated by soils formed by volcanic ash accumulation, as the páramo soils are, a combination of deeper hillslope soils (1.5–3 m depth) with highly fractured and permeable geology allows for the formation of longer flow paths of water and longer MTTs. Therefore, the relatively young and little weathered geology in the Zhurucay basin allows for a dominance of shallow subsurface flows. The results of these studies suggest that the particular shallow development of the rich organic soils with low saturated hydraulic conductivities, in combination with a homogeneous and low permeable geology, provides the páramo basin of the Zhurucay River with a high water retention capacity and relatively long transit times and flow paths considering the low development of the organic horizon of the soils. Hrachowitz et al. (2009b) reported MTTs (135–202 days) around the ones found in the Zhurucay basin catchments in a montane catchment in Scotland dominated by peatland soils and relatively little weathering geology. Nevertheless, the models that provided the best fit were the GM and the TPLR, as opposed to the EM in our study site. As in the Zhurucay basin, these authors attributed this short transit time to the dominance of ecohydrological processes occurring in the upper horizon of the peat soils. Therefore, we can conclude that in these two ecosystems, located at different latitudes but with similar hydropedological conditions, the hydrology is dominated by shallow subsurface flows. Nevertheless, the soil development of the shallow peaty soils in Scotland is lower (40 cm) in comparison to the soil development of the Histosols (80 cm) in the Zhurucay basin. These factors, in combination with differences in underlying geologies, suggest that their overall hydrologic functioning might differ, as evidenced by different TTDs describing the subsurface transport of solutes.

Although we used a methodology that assumes stationary conditions in the hydrologic system (LCA), it is relevant to note here the results of a recent application of conceptual modeling for the investigation of non-stationary conditions in a hydrologically similar region (Birkel et al., 2015). These authors detected non-stationary effects in water age distributions only during extreme weather conditions (extensive dry or wet periods) and attributed this behavior to the large mixing capacity of the Histosol soils. Although future investigations of the non-stationary nature of MTTs are needed at the páramo, based on the dominance of flow generation in the Histosols at the Zhurucay basin, in combination with low annual changes in the environmental conditions, we consider that our results from the LCA provide robust MTT estimates in our study site.

Controls on baseflow MTT variability

We found significant correlations ($R^{\mathrm{2}}\ge$ 0.78, $p$ < 0.05) between catchment slope dependent indexes and MTT using a subset of the main stream catchments (Table 7 and Fig. 6). Results of the correlation analysis indicate that (1) the higher the average slope of the catchments, the shorter their MTTs; (2) the higher the percent of area corresponding to slopes between 0 and 20 %, the longer the MTT; and (3) the higher the percent of area corresponding to slopes between 20 and 40 %, the shorter the MTT. These results indicate a clear control of the catchments' slopes in the Zhurucay basin's MTTs. Locally, the same topographical features were found to control low flow generation. Mosquera et al. (2015) attributed the latter to expected contributions from the water originating in the slopes (Andosol soils) during low flow generation as a result of the gravitational potential of the water that drains downslope from these soils. These authors also found that wetlands (Histosols soils located near the streams) control the generation of moderate and high flows. Although we did not find significant correlations with other landscape features, vegetation showed expected trends in relation to MTT. That is, catchments with a higher proportion of cushion plants (wetlands) ($R^{\mathrm{2}}=$ 0.29, $p=$ 0.35) have longer MTTs and an inverse relation to tussock grass vegetation ($R^{\mathrm{2}}=$ 0.31, $p=$ 0.33). In another tropical system of catchments in Colombia, a catchment with a higher areal proportion of wetlands was found to prolong the water MTTs, but appeared to reduce water yield (Roa-García et al., 2011). Although these authors did not report the slope of the catchments, we can infer that the catchment with the highest proportion of wetlands – as they form in flat areas – is also the catchment with the lowest gradients. Therefore, their observations might result from the combination of the deeper soil development (> 10 m) with high water retention capacity and low saturated hydraulic conductivity, perhaps in combination with low slope gradients. This would support the result of our study, where the catchments with the lower slopes and higher proportion of wetlands present the longer MTTs.

In other latitudes, in 20 Scottish catchments with different geomorphologies and climates, MTT variability was controlled by the areal proportion of peat soils, and no influence of catchment slopes was found (Hrachowitz et al., 2009a). As such, and given the similarities between these soils and our Histosol soils (Andean wetlands), we hypothesized MTT variability to be controlled by the areal extent of wetlands. Even though we found that MTT variability is rather mainly controlled by topography in our tropical alpine site, a small trend of wetland cover to increase MTT was also identified. Although the later relation is not statistically significant, the lattermost likely results from the influence of topography on Histosol soil (wetland) formation, where the formation of this soil mainly occurs in catchments with lower slopes where water accumulation is favored. This finding indirectly suggests that wetlands influence MTT spatial variability to a lesser extent. Therefore, it appears that although relatively similar processes control the ecohydrology of both ecosystems, controls on MTT variability cannot be extended from one ecosystem to the other. MTT variability was also found to be controlled by the proportion of wetlands in cold snow dominated boreal catchments in Sweden for the MTT of spring snowmelt water (Lyon et al., 2010). These authors attributed this effect to the formation of shallow ice acting as impermeable barriers above the wetlands, thus changing the flow paths of water. Nevertheless, because of the different climate and geological features between their catchments and ours, we did not find wetlands to be major controls on MTT variability.

Other slope topographic indexes – e.g., flow path length ($L$), flow path gradient ($G$), and the ratio between both ($L$ $/$ $G$) (e.g., McGuire et al., 2005; Tetzlaff et al., 2009) – have been identified as controls of MTT variability in catchments in other latitudes. Although these landscape features did not significantly explain MTT variability in the Zhurucay basin, the $L$ $/$ $G$ ratio was reported as the major control of MTT variability ($R^{\mathrm{2}}=$ 0.91) in steep temperate catchments in the central western Cascades of Oregon (McGuire et al., 2005), suggesting that this relation “reflects the hydraulic driving force of catchment-scale transport (i.e., Darcy's law)”. Similarly to our study site, they also found the average slope of these catchments to be one of the most important individual controls on MTT, explaining 78 % of the MTT variability. Recently, topography was also identified as a major control on the MTT of 12 TMCF catchments in eastern Mexico (Muñoz-Villers et al., 2016). Results from these studies reflect the fact that the integrated effect of catchment slope on MTT variation can be identified in distinct geological and hydropedological provinces. The latter also suggests that rather than using a predictor that indicates more local effects of hydraulic force driving in the stream channel (e.g., $L$ $/$ $G$), catchment slope might be a better measure to compare catchment functioning, as it integrates the hydrologic connectivity of hillslope, riparian, and stream areas.

The topographic controls on MTT in the Zhurucay basin indicate that water resides for a longer time in catchments having lower slope gradients. These results also indicate that in catchments having higher areal proportions of low gradients and lower areal proportions of steeper gradients coupled with higher wetland coverage, water resides longer in the shallow reservoir of the soils. The control of the proportion of steeper gradients in MTT variability suggests that the gravitational potential of water draining downslope in the Andosol soils also indirectly influences MTTs. Therefore, it is our interpretation that the hydrology of this ecosystem is mainly dominated by the interplay of two factors: (1) the high storage capacity in the shallow organic horizon of the porous páramo soils and (2) the catchment slope. Factor 1 drives the high water retention capacity and factor 2 controls the high regulation capacity of the ecosystem, and thus help maintain a sustained delivery of water to the streams along the year. Without the storage–slope interplay, water would remain stored in the soils, and perhaps the delivery of water towards the streams would be dominated by saturated overland flow (SOF), affecting the regulation capacity of the ecosystem. Nevertheless, SOF rarely occurs in the Zhurucay basin (Mosquera et al., 2015).

Mean electrical conductivity (MEC) was also found to be significantly correlated with MTT using all catchments of the nested system in the basin (Fig. 7). The regression analysis showed strong correlation, with MEC increasing as MTT increases. As EC is an intrinsic property of water, due to the time it spends in contact with the surrounding pore space, rather than a control on MTT variability, this result indicates that this property might be used as a proxy to estimate MTT spatial variability. The well-defined connection between MTT and MEC most likely results from the relatively homogenous geology of the Zhurucay basin. To our knowledge, there are no studies that have identified similar (or different) relations between MEC and MTT in other biomes.

Correlation between mean transit time (MTT) and mean electrical conductivity (MEC) for weekly measurements of stream-water samples collected during 3 years (June 2012 to June 2015). The solid line is the linear regression and the dashed lines are the 90 % confidence intervals of the regression. ${}^{*}$ indicates variables are normalized by their means.

[Figure omitted. See PDF]

Given that estimating MTT using isotope tracer methodologies is financially expensive due to the logistical setup of a monitoring network and the processes of data collection and analysis, finding proxies (i.e., predictors) that allow one to infer MTTs at lower operational costs is critical to improve water resource management. In this sense, the strong relation between MEC and MTT indicates that MEC could be used as a relatively inexpensive and directly measurable proxy for MTT in this wet Andean páramo catchment. Therefore, although this result cannot be expanded beyond páramo areas, perhaps not even beyond the study site, it seems that it is worth evaluating whether or not MTTs can be inferred from MEC in other hydrologic systems. Nevertheless, one should be careful that EC measurements can be relatively variable over time. As a result, a single measurement of EC is most likely not enough to provide robust MTT estimates. Therefore, longer EC measurement records are recommended to reduce MEC variability and provide more robust MTT estimates.

Conclusions

The baseflow MTT evaluation using a LCA indicated that the EM best describes the subsurface transport of water in the basin. This result indicates efficient mixing in the high organic and porous wet Andean páramo soils and a simple subsurface transition of rainfall water towards the streams. MTT estimations showed relatively short MTTs linked to relatively short subsurface flow paths. Therefore, we confirm that the hydrologic system of the tropical alpine biome of the Zhurucay basin is dominated by shallow subsurface flow. MTT estimations showed that catchment M7, located at a flat hilltop at the outlet of a wetland that remains ponded year-round and disconnected from the slopes – most likely as a result of the eutrophication of a lagoon – showed a particularly low MTT (0.15 years, 53 days) in relation to the MTT in all of the other catchments (0.40–73 years, 156–250 days) in which the morphology corresponds to U-shaped valleys, with the wetlands located at the valley bottoms near the streams and connected to the slopes. Two headwater catchments, M3 and M4, showed the longest MTTs, related to a small contribution from a spring shallowly sourced. These results indicate that in this páramo ecosystem, the geomorphology of the wetlands and geology, to a lesser extent, influence the responsiveness of the streams to precipitation inputs. Correlation analysis between landscape variables and MTT indicates that MTT variability is mainly explained by the slope of the catchments, and a related influence of vegetation to a lesser extent. Catchments with the steepest average slopes and lower proportion of wetlands have the shortest MTTs. The lack of significant correlations between the MTTs and hydrological response variables (runoff coefficient and specific discharge rates) indicates that neither water yield nor streamflow rates control the time water resides in the subsurface of the páramo soils. These results indicate that the interplay between the high storage capacity of the páramo soils and the slope of the catchments defines the ecosystem's high regulation capacity. Mean electrical conductivity (MEC) of stream waters – with the oldest waters presenting the highest MECs – seems to be a promising proxy of MTT in systems of catchments under homogeneous geological conditions. Finally, we want to highlight the usefulness of a nested monitoring system for acquiring better process-based hydrologic functioning understanding. For instance, if M3, M4, and/or M7 catchments would not have been monitored, the influence of geology and/or geomorphology on catchment hydrological response could not have been identified, and important information about the whole ecosystem functioning would remain unknown.

Data availability

The isotope data are property of the Central Research Office at the University of Cuenca (DIUC) and the Ecuadorian National Secretariat of Higher Education, Science, Technology and Innovation (SENESCYT). These data will become publicly available in accordance with the rules and embargo regulations of the funding agencies, but it is not yet known where the data will be hosted (please contact the corresponding author for updates). The codes and model outputs used in this study can be requested by contacting the corresponding author ([email protected], [email protected]).

Acknowledgements

The research was funded by the Central Research Office at the University of Cuenca (DIUC) via the project “Desarrollo de indicadores eco-hidrológicos funcionales para evaluar la influencia de las laderas y humedales en una cuenca de páramo húmedo”, and the Ecuadorian National Secretariat of Higher Education, Science, Technology and Innovation (SENESCYT). The authors express gratitude to INV Metals S.A. (Loma Larga Project) for the logistic support. Special thanks to Alicia Correa and Patricio Lazo for their fieldwork assistance, Irene Cardenas for the laboratory analysis, and Steven Perakis, Christian Birkel, and an anonymous reviewer for their valuable comments on earlier versions of the manuscript. Edited by: C. Stumpp Reviewed by: C. Birkel and one anonymous referee