SEHR-ECHO v 1 . 0 : a Spatially-Explicit Hydrologic Response model for ecohydrologic applications

Introduction Conclusions References


Introduction
Hydrological processes result from natural processes that vary strongly in space, such as precipitation, evaporation or infiltration into the subsoil (McDonnell et al., 2007;Beven, 2012).Accordingly, most state-of-the art hydrologic response models have two fundamental components to describe the arrival of water and transported substances at a control section: a component to simulate the temporal evolution of water storage (and possibly of energy or of solutes) and released fluxes in some hydrologically meaningful sub-units and a component to describe the transport of the fluxes between the sub-units along the river network, (e.g.Hingray et al., 2010;Clark et al., 2008;Kunstmann and Stadler, 2005).There are, however, many models without an explicit description of the water flow within the landscape and the river network.This transport component is either completely omitted as in lumped models (Perrin et al., 2003;Merz and Blöschl, 2004), assumed to be negligible at the spatio-temporal scale of interest (Viviroli et al., 2009;Schaefli et al., 2005), or assumed to fall out of the sum of transport processes simulated between small spatial units, without further parameterizing flow in channels (e.g.Tague and Band, 2001;Wigmosta et al., 1994;Liu and Todini, 2002).Some models use an arbitrary (calibrated) discharge routing function to smooth the response computed at the scale of a (sub-)catchment, e.g. with a triangular function as in the HBV model (Bergström, 1995) and derivations thereof (Wrede et al., 2013;Das et al., 2008).Finally, there exist models that can be thought of as having an implicit routing component (Tague and Band, 2001) even if they are applied in a completely lumped manner, i.e. a component that accounts statistically for spatial differences of runoff generation, such as Hymod (Boyle, 2000;Moradkhani et al., 2005) or the well-known Topmodel with its topographical wetness index (Beven and Kirkby, 1979) (which is, however, generally not applied in a lumped manner).
Most existing catchment-scale model applications show an explicit parameterization of channelled flow only for the largest rivers in the analysed system, for which typically a kinematic wave-based routing is used, assuming that at the subcatchment level, the effect of travel times in channels are negligible at the considered spatio-temporal resolution.This might typically hold for hourly discharge simulation with subcatchments of a few 10-100 km 2 in reasonably steep environments (e.g.Hingray et al., 2010).For an example including subcatchment routing parameterization see, e.g., the work of Clark et al. (2008).
Traditionally, the range of aforementioned hydrological models is classified into (semi-)lumped and (semi-)distributed (Reed et al., 2004;Beven, 2012), a classification which refers essentially to the parameterization of the temporal evolution of the water storage within the catchment.The terms distributed or semi-distributed (e.g.Das et al., 2008) generally refer to grid-based models or subcatchment set-ups with different parameter sets for each spatial unit, whereas semi-lumped implies some degree of spatial discretization but with a single parameter set and generally without flow routing through the landscape (e.g.Schaefli et al., 2005).Although not directly related, it is often implied that distributed models are more physics-based.
We propose the term of a spatially explicit hydrologic response (SEHR) model for any model that explicitly parameterizes both, spatial patterns of water storage evolution as well as the effect of geomorphology on the travel time of water having different spatial origins.We believe that the term spatially-explicit is more generic than the often used terms semi-lumped, semi-distributed or distributed model, which refer to specific set-ups in terms of spatial variability of state variables and of parameters.
Hereafter, we describe a simple catchment-scale hydrologic model developed at the Laboratory of Ecohydrology (ECHO) of the Ecole Polytechnique Fédérale de Lausanne, SEHR-ECHO, that explicitly accounts for the spatial variabilities in the runoff generation process and the heterogeneity of the flow-paths within the catchment.The model builds on the geomorphic theory of the hydrologic response (Rodriguez-Iturbe and Valdés, 1979;Rodriguez-Iturbe and Rinaldo, 1997;Rinaldo et al., 2006) pursuing an accurate de-scription of riverine hydraulic processes through the use of the geomorphologic dispersion (Rinaldo et al., 1991), providing a general framework to formulate spatially explicit models (e.g.Nicótina et al., 2008;Tobin et al., 2013).In such an approach, nonlinearities and nonstationarities of the hydrologic response (e.g.McDonnell et al., 2010;Botter et al., 2011;Sivapalan et al., 2002;Hrachowitz et al., 2013) are embedded in the parameterization of the soil moisture dynamics and the related dominant runoff generation processes at the source area scale.
The general model concept is introduced in Sect. 2 and implementation details discussed in Sect.3.For illustration purposes, the model is applied to an example case study from Switzerland (Sect.4), for which we discuss the discharge simulation performance in Sect.5, before summarizing the main conclusions (Sect.6).

Model description
The SEHR-ECHO model is composed of two main components (Fig. 1): (i) a precipitation-runoff transformation module that computes surface and subsurface water fluxes from the source areas (the basic sub-units that describe the spatial structure of the model domain), and (ii) a routing module that computes fluxes in the river network through to the control section (i.e. the outlet).In other terms, the model is composed of a module for unchannelled state processes at the source area scale and one for channelled state transport (Rinaldo et al., 2006).
The source areas are extracted from a digital elevation model (DEM) with the well-known Taudem algorithm (Tarboton, 1997) for subcatchment and river network delineation (see Sect. 4 for further details).The scale of these source areas are selected such as to allow for sufficiently homogeneous hydro-meteorological conditions without losing too much geomorphologic complexity.Relevant geomorphologic issues are discussed in Rodriguez-Iturbe and Rinaldo (1997).

Precipitation-runoff module at source area scale
The precipitation-runoff module solves the mass balance equations at the source area scale.This component is driven by precipitation, temperature and potential evaporation input time series, which need to be properly provided at the source area scale.The choice of methods to interpolate the observed input time series to this scale depends on the variable (precipitation, temperature), on the application and on the simulation time step (e.g.Tobin et al., 2011), and the choice of the general method is largely independent of the exact set-up of the hydrological model.In the remainder of this section we describe the precipitation-runoff transformation for one source area assuming the input variables are provided at the proper spatial and temporal scales.The precipitation-runoff transformation module has the following key elements: interception and re-evaporation of intercepted water, rainfall/snowfall separation, evolution of water stored in the snowpack in solid form, evolution of the liquid water content of the snowpack and equivalent precipitation (rain and meltwater)-runoff transformation.If a source area has partial glacier cover, the runoff resulting from the glacier is computed separately (Fig. 1).
Interception I c (t) is simulated using a constant interception capacity ρ: where P (t) [L T −1 ] is the precipitation and ρ [L T −1 ] is the maximum interception during a time step t (e.g.Feni-cia et al., 2006).No separation between the aggregation state of precipitation (snow, rain) is made for ρ, a simplification which is not advisable for applications to catchments with considerable forest cover (e.g.Gelfan et al., 2004).Part of the intercepted water is assumed to re-evaporate during the same time step, limited by potential evaporation E pot (t) [L T −1 ]: where E i (t) [L T −1 ] is the evaporation flux from intercepted water.This evaporated water is assumed to not be available for the precipitation-runoff generation process, i.e. total incoming precipitation is reduced to net precipitation P n [L T −1 ]: E pot is reduced to potential transpiration E t,pot [L T −1 ] by the amount of E i : (4) It is noteworthy that the above formulation assumes that interception is an instantaneous process, which takes place at timescales smaller than the simulation time step (i.e.subhourly).Only the evaporated water is subtracted from the incoming precipitation, which corresponds to a return of nonevaporated water as throughfall.
The estimation of the aggregation state of precipitation is based on a simple temperature threshold T r [ • C] (Schaefli et al., 2005) that splits net precipitation P n into rainfall P r (t) [L T −1 ] and snowfall P s (t) [L T −1 ] depending on the mean temperature T (t) (for a smooth threshold approach see Schaefli and Huss, 2011).
The evolution of the water equivalent of the snowpack height h s [L] is computed as where M s (t) [L T −1 ] is the snowmelt due to energy input from the atmosphere, F s (t) [L T −1 ] is the flux of refreezing water during periods of negative heat input and G s (t) [L T −1 ] is the snowmelt due to ground-heat flux (all in water equivalent).G s is assumed to be constant in time, G s = G max as long as there is a snowpack: where a s is the degree-day factor for snowmelt [L T −1 • C −1 ] and T m [ • C] is the threshold temperature for melting that is set to 0 • C. Refreezing F s is assumed to occur when T (t) < T m and is linearly related to negative air temperature with a freezing degree-day factor that is proportional to, but smaller than, the melt degree-day factor a s (Kokkonen et al., 2006;Formetta et al., 2014): where a f [0, 1] is the degree-day reduction factor for refreezing.
The snowpack is assumed to have a certain water retention capacity θ •h s [mm] and water is released from the snowpack only if this retention capacity is exceeded.The balance equation for the liquid water h w [L] content of the snowpack is written as where the snowpack outflow M w only occurs if the air temperature is above melting conditions and if the water retention capacity θ is reached: It is noteworthy that rainfall P r only occurs if T > T r , snowmelt M s only if T > T m and refreezing F s only if T < T m .It generally holds that T r > T m = 0 • C, translating the well-known fact that snowfall can occur at temperatures above the melt temperature (for an analysis of observed snowfall at the Davos station, see Rohrer et al., 1994).The general case of Eq. ( 10) holds for all values of the threshold parameters or for fuzzy transitions between rain-and snowfall.
The water fluxes P r , G s and M w are summed up to produce the equivalent precipitation P eq that enters the equivalentprecipitation-runoff transformation: The partitioning of equivalent precipitation into surface runoff, fast and slow subsurface runoff and transpiration resulting from water infiltration and percolation in the subsoil is performed via a minimalist description of the soil moisture dynamics at the source area scale (Laio et al., 2001;Rodriguez-Iturbe and Porporato, 2004;Rodriguez-Iturbe et al., 1999): where η [−] is the soil porosity, Z r [L] is the depth of the soil layer that is active during water redistribution processes, s [−] is the relative soil moisture in the active layer, is the rate of transpiration of water from the root zone and L [L T −1 ] is the water flux (called leakage here) mobilized from the root zone as subsurface flow.
It is noteworthy that this soil moisture dynamics equation, if forced with Poisson infiltration, can be solved exactly for a number of cases and forms the basis of substantial analytic work on the probabilistic properties of stream flow (Botter et al., 2007a, b;Botter, 2010).
The leakage L is parameterized as a nonlinear function of the soil moisture as where K sat [L T −1 ] is the saturated hydraulic conductivity and c the Clapp-Hornberger exponent (Clapp and Hornberger, 1978).
The transpiration rate is a linear function of relative soil moisture between the wilting point, s w [−] (i.e. the moisture content below which the plants cannot further extract water from the soil) and the upper limit of water stress, s m [−], at which it is assumed to reach the limit imposed by the potential transpiration rate (e.g Porporato et al., 2004): The infiltrated water corresponds to the equivalent precipitation from which direct surface runoff is subtracted: where R hort [L T −1 ] is surface runoff occurring if the infiltration capacity is exceeded and R dun [L T −1 ] is surface runoff occurring if the source area is saturated.These two mechanisms of surface runoff, inspired from Hortonian and Dunne overland flow (e.g.Dingman, 2002), enable the model to simulate different timescales of reaction to a precipitation or melt water input.R hort is parameterized with a constant maximum infiltration capacity φ [L T −1 ]: where φ is supposed to be constant in time.If the soil is saturated, s(t) = 1, R dun occurs: The water mobilized from the active layer L is transformed to subsurface runoff at the source area outlet through two linear reservoirs that simulate a fast and a slow subsurface flux, R fast [L T −1 ] and R slow [L T −1 ] (similar to e.g. the formulation in the HBV model Bergström, 1995).The part of L feeding the slow subsurface flux, L slow , is assumed to be a constant flux L max limited by L: The linear reservoir equations for the fast and slow subsurface runoff thus read as Note that s is a relative soil moisture, whereas S slow and S fast have length units.

Discharge simulation from glacierized subcatchments
If a source area has a partial glacier coverage, ice is assumed to start melting if h s = 0 (Schaefli et al., 2005): where a i [L T −1 • C −1 ] is the degree-day factor for ice melt.This melt is routed to the subcatchment outlet through a linear reservoir with coefficient k ice .The routing to the catchment outlet follows the general procedure (see hereafter) but the flux is weighted according to the fraction of source area that is glacier covered.No glacier surface dynamics are modelled (e.g Huss et al., 2010) and the glacier cover is assumed to be constant (but it can of course be updated for different simulation periods).

Discharge simulation at catchment outlet
The transport of the runoff components through the river network uses a linear approach, assuming that most relevant nonlinear processes are captured through the source areascale precipitation-runoff transformation.This assumption only holds for systems where flow velocity can be assumed to be relatively constant in time (independent of discharge) and space.The total discharge at the catchment outlet is obtained by convolution of each of the fluxes R (surface runoffs R hort , R dun , subsurface runoffs R fast , R slow and ice melt runoff R ice ) from all source areas with a travel time distribution f γ (t) along its flow path γ to the catchment outlet (Fig. 2).
Pathways from source areas A j to outlet: . Sketch of the flow paths from a catchment with five source areas A j and three channels C j .The notation C * j means that the injection into this channel is not concentrated at the upstream end but, in theory, randomized and integrated over the channel length (Bras and Rodriguez-Iturbe, 1985).In practice, we take half of the channel length.
The probability density functions of travel times, f γ (t) (assumed statistically independent) are obtained from the travel time distributions in all channels C j → C k → C composing them (Gupta et al., 1980;Rinaldo et al., 1991): where * is the convolution operator and is the outlet.For the example illustrated in Fig. 2, the travel path from the source area A 1 to the outlet is made up of two collecting channels C 1 and C 3 .
The travel time distributions within channelled states C j are obtained assuming longitudinal 1-D dispersion, which is a reasonable assumption for open channel flow in low-order rivers (Rinaldo et al., 1991): where D [L 2 T −1 ] is the hydrodynamic dispersion coefficient, j [L] is the channel length and ν [L T −1 ] is the average velocity.
The simulated discharge at the catchment outlet becomes www.where each of the discharge components Q xyz equals 3 Model implementation The model requires temperature, precipitation and potential evaporation time series for each subcatchment and, for model calibration, at least one concomitant discharge time series observed at the catchment outlet.The numerical implementation uses a fixed time step and a fourth-order Runge-Kutta scheme to compute the soil moisture evolution.The other stores (fast and slow subsurface flux stores, solid and liquid snow stores) are solved with explicit time stepping, which is justified given that these stores have only one outflux, linearly dependent on the storage.The provided code (see Sect. "Code availability") is developed in Matlab R2010b.The parameterization of each of the presented hydrological processes can easily be modified.The basic model structure (passing of variables and parameters among functions) has been designed for an easy combination with the now widely used optimization algorithms developed by Vrugt et al. (2003Vrugt et al. ( , 2009)).For the example presented in this paper, the model is, however, calibrated with simple Monte Carlo generation within a priori parameter ranges (details in Sect.4).

Identification of model parameter patterns
The physical parameters of SEHR-ECHO, which describe the physiographic catchment characteristics and can be extracted from topographic data, are listed in Table 1.The model parameters that require calibration or a relevant method of a priori estimation are summarized in Table 2 (along with a range of a priori values).Depending on the application, all these parameters might be made variable in space, especially the ones for which there are known spatial patterns.
It is, in particular, recommended to relate the mean residence time k −1 fast of subcatchment-scale subsurface fluxes to the source area A γ , as where the scaling coefficient ξ can be set to values of around 1/3 (Alexander, 1972;Pilgrim et al., 1982).In practical terms, such a scaling for k fast is obtained by calibrating a generic reservoir coefficient k cal such that where A γ is the mean subcatchment area.The coefficient k slow is then related to k fast through the calibration of a multiplication parameter m k such that It might be tempting to derive the spatial variability of the active soil depth from soil production theory (Heimsath et al., 1997).Such an approach might e.g.assume that the soil depth of a source area is proportional to the mean topographic curvature in topographically convex areas.Another idea could be to identify topographically concave areas that can be assumed to be saturated at all times.For applications similar to the one presented here, different model tests showed, however, that the effect of spatially variable soil depth on simulated discharge can be compensated by the other model parameters (Nicótina et al., 2011); this approach is therefore not further pursued here.
The saturated hydraulic conductivity can easily be distributed in space according to observed land use and soil types; an example is discussed in Sect. 4 for the present case study.Imposing additional spatial parameter patterns related to directly observable physiographic characteristics is readily possible, but beyond the scope of the present work.

Case study
The Dischmabach catchment, located in the southeast of Switzerland near Davos (Fig. 3), has a size of 43.3 km 2 at the Kriegsmatten gauging station, for which long discharge time series are available from the Swiss Federal Office for the Environment.Its elevation ranges from 1668 up to 3146 m a.s.l.(mean altitude 2372 m a.s.l.) with around 2.1 % of the catchment area covered by glaciers.The annual mean temperature at mean elevation is around −0.5 • C. The discharge regime shows a strong seasonal pattern due to accumulation and melting of snow.The relatively steep hillslopes are covered with pasture (38 %) and forest (10 %); around 16 % of the catchment are bare soil, rock outcrops cover 24 % (Verbunt et al., 2003).The geology is crystalline composed of gneiss and amphibolites, overlain by shallow soils (Verbunt et al., 2003).The nearby meteorological station of Weissfluhjoch (2690 m a.s.l.) records around 1450 mm year −1 of precipitation (period 1981-1999), which is relatively low compared to other Alpine locations at the same altitude.The discharge over the same period was around 1350 mm year −1 .The mean  evaporation in this catchment is around 300 mm year −1 for the period 1973-1992 (Menzel et al., 1999).Part of the discharge is due to net glacier ice melt.The subcatchments as well as the river network characteristics required to run the model (network topology, river reach lengths) are identified with TauDEM Version 5 (Tarboton, 1997), a hydrologic terrain analysis tool which is freely available under the terms of the GNU General Public License version 2 at http://hydrology.usu.edu/taudem/taudem5/.The 23 subcatchments as well as the river network are shown (Fig. 3).
The temperature time series for each subcatchment is obtained through a linear interpolation of the temperature ob-served at the Davos weather station (1594 m a.s.l.) to the mean subcatchment altitude, using the average temperature lapse rate between this and the Weissfluhjoch station (which equals −0.50 • C/100 m).The precipitation for each subcatchment is the one recorded at Weissfluhjoch.The potential evaporation is evaluated with the Priestley-Taylor method (Maidment, 1993;Priestley and Taylor, 1972) S1).Based on the distribution of r j within each subcatchment γ , we compute the following scaling parameter: where γ identifies a given subcatchment, f j is the relative frequency of occurrence of the land use class j within the subcatchment and r D is the surface runoff coefficient of the dominant land use class.
For each subcatchment, the saturated hydraulic conductivity is then obtained as where K sat is obtained through calibration.Impervious subcatchment areas are accounted for by setting the soil depth of the corresponding portions to 0 (in total 1.2 km 2 ).

Model calibration
For the purpose of this paper, the model is calibrated on daily and hourly discharge with simple Monte Carlo simulation: we draw a high number of random parameter sets in the a priori parameter ranges and retain the best simulations with respect to the well-known Nash-Sutcliffe performance criterion (Nash and Sutcliffe, 1970), which evaluates how much better the simulated discharge Q s fits the observed discharge Q o than the simplest possible model, the mean of the observed discharge over the entire period Q o : where N is the Nash-Sutcliffe efficiency, NSE, and t i the ith time step, i = 1, . .., n t .In addition, we analyse the NSE-log value computed on log-transformed discharges (N L ) and the relative bias between the simulated and the observed mean discharge.
For hydrological regimes with a strong annual discharge cycle, the above NSE value is not very meaningful since any model that reproduces the annual cycle more or less will have a high NSE value (Schaefli and Gupta, 2007).We therefore compute the benchmark NSE value, N (Q b ), for a benchmark model which corresponds to the average of all observed discharges on a given time step k of the year y (either a Julian day or an hour of a Julian day) (Schaefli and Gupta, 2007): where k y is the kth time step of year y and Y the total number of years.For the observed discharge of the Dischmabach, N(Q b ) equals 0.74 at the daily time step and 0.73 at the hourly time step.For N L (Q b ) the values of the benchmark is 0.87 for the daily and the hourly time step.Given the insignificant role of Horton direct runoff in this environment, this runoff mechanism is deactivated here (assuming infinite infiltration capacity).For all other processes, the a priori parameter ranges are obtained based on existing literature (see Table 2).The upper limit of the percolation flux feeding the slow subsurface flux, L max , is chosen equal to the mean daily precipitation.The upper limit of the average root zone depth is fixed to 1.5 m, which is a conservative estimate given the type of vegetation present.The generic fast subsurface residence time, k cal , Eq. ( 27), is assumed to be of the order of magnitude of days and the slow subsurface residence to be substantially longer (maximum scaling factor m k = 50).
A similar scaling approach is also used to ensure that the degree-day factor for ice a i is higher than the one for snow a s and that the retention capacity s m is higher than the wilting point s w .

Results
For the Dischma catchment, a total of 12 model parameters have to be calibrated, seven for the water input-runoff transform and five for the glacier and snowmelt simulation.Here, these calibration parameters have been estimated through simple Monte Carlo simulation to illustrate the main features of the SEHR-ECHO model.Figure 4 shows the discharge simulation along with the simulated series of evapotranspiration and observed meteorological time series at the catchment outlet for the best NSE parameter set obtained with 35 000 parameter sets sampled uniformly within the prior distributions of Table 2.With a NSE value of 0.82 during calibration, this parameter set performs better than the benchmark model (Sect.4).It furthermore gives reasonable estimates for evapotranspiration, which indicates that the observed precipitation time series is an acceptable proxy for catchment-scale area-average precipitation input.
The splitting between the three hillslope scale runoff generation processes corresponds to the expected pattern: Fig. 5 illustrates that the slow subsurface component contributes essentially to base flow and that the direct surface runoff is activated only occasionally.It is noteworthy, however, that this pattern results partially from the imposed subsurface residence time scaling.(The corresponding subcatchment scale state variables are given in Fig. S1 of the Supplement along with an example of the simulated snowpack evolution, Fig. 2

.)
This parameter set comes along with a number of sets that lead to equally good discharge simulations for the reference performance criteria for the calibration period.The plots of NSE vs. NSE-log and of NSE vs. the relative bias (Fig. 6) illustrate that these performance criteria can be well optimized simultaneously, which is not always the case for hydrologi-  3. cal models.Such models typically show a strong tradeoff between NSE and NSE-log because the same set of processes cannot reproduce high flows and recessions.This problem is reduced in the analysed hydrologic regime where low flows are dominated by the long winter recession, relatively simple to simulate (see also the log discharge plot in Fig. 4).
A notable aspect of this SEHR-Echo application is that a large number of the best-performing parameter sets at the daily time step perform equally well during the calibration period (1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992) and the validation period (1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000) and at the hourly time step (Fig. 6).This is also illustrated in Fig. 7 which shows the ensemble of discharge simulations obtained for the 100 best daily parameter sets (see Fig. 8) applied at the daily and the hourly simulation time step.The corresponding prediction limits span the observed discharge equally well at both time steps.The simulation bias increases however at the hourly time step (Table 3): this is most likely due to the assumption of a constant within-day relation between air temperature and snowmelt and refreezing and ice melt; the related parameters might require a specific calibra- tion to the hourly time step or a more appropriate formulation for the hourly time step (Tobin et al., 2013).
A comparable assessment of model performance at different timescales without re-calibration is rarely reported in the literature.For an example, see the work of Schaake et al. (1996).
The above evidence, timescale independence and splitting between the runoff generation processes, are important hints that the model works well for the right reasons: the parameters play the role they have been designed for, rather than trying to mimic omitted runoff generation processes or to compensate for the lack of spatial differentiation of travel times, which might typically occur for lumped models.
This conclusion is also supported by the good identifiability of some of the model key parameters, illustrated by the relatively peaky distributions in Fig. 8, which shows histograms of the best 100 parameter sets (in terms of NSE) at the daily time step.Albeit not providing a formal assessment of model and parameter uncertainty, this simple analysis illustrates that the model has a relatively well defined range of optimal parameters, which might be further refined for realworld applications and specific questions (e.g.extreme event analysis).It is noteworthy, however, that a flat distribution of the best-performing parameter sets for a given specific parameter (e.g.here the wilting point) does not point towards model insensitivity with respect to this parameter, since its relation to other parameters might simply not be visible in the marginal distribution (Bardossy, 2007).The question arises as to how far the geomorphologybased set-up and the routing scheme influence the results.The model simulations are insensitive to hydrodynamic dispersion since its effect is overruled by advection.Given the relatively short distances in this catchment (the longest travel paths from a subcatchment to the outlet is 11 km), the velocity of in-stream discharge routing has only a minor effect at the daily timescale and a notable influence of the velocity on the simulated discharge would be obtained only with unrealistically low velocities (e.g.Yochum et al., 2012) (Fig. 9).For the hourly time step, any variation of flow velocity affects the discharge simulation and including this parameter in the calibration process might be a possible option, keeping in mind, however, that it might interact with the recession coefficients.3 and velocity of 0.5 m s −1 ; the sensitivity is expressed as the relative difference to the reference simulation for each time step.

Comparison of the subcatchment set-up to elevation bands
Comparable precipitation-runoff models often either use a grid-based spatial discretization (e.g.Huss et al., 2008) or an elevation band approach (e.g.Schaefli et al., 2005;Stahl et al., 2008) to account for temperature gradients as the strongest spatially variable driver.Such an approach can be assumed to yield a more reliable representation of the snow accumulation and melt processes but it necessarily leads to a description of the equivalent precipitation-runoff transformation that cannot properly account for the spatial origin of flow.We compared the performance of a model set-up where the catchment is subdivided into 10 elevation bands (Fig. 3) to the subcatchment based set-up.A look at the model performance in terms of NSE and NSE-log for both set-ups (Fig. 10, top left) demonstrates that the elevation band approach marginally outperforms the subcatchment approach for the calibration period at the daily time step.If the best 100 of these parameter sets (in terms of NSE efficiency) are applied to the hourly time step (Fig. 10, top row, centre column), the elevation band set-up does a noticeably better job.This higher performance however disappears for the validation period (Fig. 10, top right), which is a strong hint of overfitting during the calibration period, where the calibration parameters might compensate for the lack of a proper accounting for the spatial origin of flows.This hypothesis is supported by the fact that if we use subcatchments divided into elevation bands, we obtain a consistent improvement of the discharge simulation performance at the hourly time step for the calibration and the validation period (Fig. 10, bottom row).

Conclusions
This paper presents a precipitation-runoff model that computes spatially explicit water fluxes at the ecosystem level and that can, thus, be used as a simulation tool for ecohydrologic applications requiring distributed discharge information.The model formulates the hydrologic response of a catchment as a convolution of the subcatchment-scale hydrologic flow processes with the river network, where the kernels account for the spatial arrangements of the subcatchments linked by the river network.The hydrologic response accommodates directly any direct information on observable physiographic catchment characteristics such as in-stream flow path lengths or subcatchment area as a proxy for subsurface residence time scaling.Remaining model parameters are calibrated on observed discharge.This spatially explicit parameterization confers the model transferability across timescales, as has been demonstrated in this paper based on a cryosphere-dominated catchment from the Swiss Alps where, due to the steep topography, travel times in unchannelled areas are dominating in-stream travel times.The main focus of this paper was on discharge simulation.Including appropriate formulations of subcatchment-scale mass transformation processes, the general modelling framework can be extended to transport processes.The Supplement related to this article is available online at doi:10.5194/gmd-7-2733-2014-supplement.

Figure 1 .
Figure 1.Flow diagram of the precipitation-discharge computation.The grey boxes highlight the model output time series.
where S fast[L]  and S slow[L]  are the water storage in the fast and the slow reservoirs.R fast [L T −1 ] and R slow [L T −1 ] are the fast and slow reservoir outflows, which are supposed to linearly depend on the storage, i.e.R fast = k fast S fast and R slow = k slow S slow , where k −1 fast [T] and k −1 slow [T] are the mean residence times.

Figure 3 .
Figure 3. Location of the Dischmabach catchment within Switzerland(source: SwissTopo, 2005(source: SwissTopo,  , 2008) )  and the 23 subcatchments identified with TauDEM Version 5(Tarboton, 1997).The 10 elevation bands in the right plot are used for comparison purposes in Sect.5.1.The latitude and longitude are indicated in the Swiss coordinate system (in km).
using the Weissfluhjoch meteorological data.Given the important heterogeneity of land use in this catchment, we distribute the saturated hydraulic conductivity according to land use types, which are available from the Swiss land use database at a resolution of 100 m (Swiss Federal Office for Statistics, 2001).We assign each relevant land www.geosci-model-dev.net/7etal.: The hydrologic response model SEHR-ECHO v1.0 use class j a surface runoff coefficient r j (see Supplement, Table

Figure 4 .
Figure 4. Observed and simulated hydro-meteorological time series of model state variables for the parameter set of Table3.

Figure 5 .
Figure 5.Time series of the streamflow components (direct surface flow, fast subsurface flow and slow (deep) subsurface flow) corresponding to the parameter set of Table 3 and the discharge plot of Fig. 4.

Table 3 .Figure 6 .
Figure 6.NSE vs. NSE-log for all simulated parameter sets during model calibration with daily time step (period 1981-1994), values for 100 parameter sets with best NSE values, values for the same parameter sets simulated over validation period (1993-2003) and values for the same parameter sets simulated at hourly time step over the validation period.The lines indicate the benchmark values for daily and hourly (broken line) time steps.

Figure 7 .
Figure 7. Predictions resulting from the 100 best parameter sets obtained under NSE for a daily time step (of the total 35 000 Monte Carlo simulations); top: simulation at the daily time step, bottom: simulation at the hourly time step (same parameter sets).The first half of the plot shows that the prediction limits are reasonably narrow, the reserved plotting order in the second half shows that the observations are well spanned -79 % (daily) respectively 75 % (hourly) of the observed time steps fall into the range spanned by the simulations.

Figure 8 .
Figure 8. Parameter distributions obtained for the best 100 parameter sets under NSE at a daily time step.The y-axis shows the relative frequency of parameter values in each bin.s m , k slow respectively a i result from a multiplication of s w , k fast respectively a s with the distribution of the corresponding calibrated scaling parameters.

Figure 9 .
Figure 9. Sensitivity of the discharge simulation with respect to the in-stream flow velocity plotted against the reference simulation with parameters of Table3and velocity of 0.5 m s −1 ; the sensitivity is expressed as the relative difference to the reference simulation for each time step.

Figure 10 .
Figure 10.Comparison of model performance of the subcatchmentbased model set-up (23 subcatchments) to elevation band set-ups.Top row: comparison to 10 elevation bands (see Fig. 3), bottom row: comparison to a combination of elevation bands and subcatchments (5 bands for each of the 23 subcatchments).Left column: distribution of the Nash-Sutcliffe (NSE) efficiency computed at a daily time step for the calibration period (shown are parameter sets with NSE > 0); centre column: NSE distribution of the 100 sets with the highest NSE values at the daily time step run at hourly time step for the calibration period; right column: same 100 sets run at an hourly time step for the validation period.The same 35 000 parameter sets are run for all three model set-ups.

Table 1 .
Model parameters describing the physiographic catchment characteristics.The units are given in generic terms (L: length).

Table 2 .
Model parameters that have to be calibrated or estimated otherwise with indication of a priori values, a reference value for parameters that are not calibrated here and key references for further details.A total of 12 parameters are calibrated here.The maximum value of the reservoir coefficient is the numerical time step.