EcH 2 O-iso 1.0: water isotopes and age tracking in a process-based, distributed ecohydrological model

. We introduce EcH 2 O-iso, a new development of the physically based, fully distributed ecohydrological model EcH 2 O where the tracking of water isotopic tracers ( 2 H and 18 O) and age has been incorporated. EcH 2 O-iso is evaluated at a montane, low-energy experimental catchment in northern Scotland using 16 independent isotope time series from various landscape positions and compartments, encompassing soil water, groundwater, stream water, and plant xylem. The simulation results show consistent isotopic ranges and temporal variability (seasonal and higher frequency) across the soil proﬁle at most sites (especially on hillslopes), broad model–data agreement in heather xylem, and consistent deuterium dynamics in stream water and in groundwater. Since EcH 2 O-iso was calibrated only using hydrometric and energy ﬂux datasets, tracking water composition provides a truly independent validation of the physical basis of the model for successfully capturing catchment hydrological functioning, both in terms of the celerity in energy propagation shaping the hydrological response (e.g. runoff generation under prevailing hydraulic gradients) and ﬂow velocities of water molecules (e.g. in consistent tracer concen-trations at given locations and times). Additionally, we show that the spatially distributed formulation of EcH 2 O-iso has the potential to quantitatively link water stores and ﬂuxes with spatiotemporal patterns of isotope ratios and water ages. However, our case study also highlights model–data discrepancies in some compartments, such as an over-dampened variability in groundwater and stream water lc-excess, and over-fractionated riparian topsoils. The adopted minimalistic framework, without site-speciﬁc parameterisation of isotopes and age tracking, allows us to learn from these mismatches in further model development and benchmarking needs, while taking into account the idiosyncracies of our study catchment. Notably, we suggest that more advanced conceptualisation of soil water mixing and of plant water use would be needed to reproduce some of the observed patterns. Bal-ancing the need for basic hypothesis testing with that of improved simulations of catchment dynamics for a range of applications (e.g. plant water use under changing environmental conditions, water quality issues, and calibration-derived estimates of landscape characteristics), further work could also beneﬁt from including isotope-based calibration.

dynamics -either measured or modelled -characterising isotopic composition provides powerful insights into water pathways at scales ranging from the pedon (Sprenger et al., 2018) to the catchment landscape (McGuire and McDonnell, 2006;. At larger scales, such approaches can yield global estimates of terrestrial water flux partitioning (Good et al., 2015), where recent scrutiny has been brought upon separating plant transpiration from other sources of evaporative losses (e.g. Jasechko et al., 2013;Coenders-Gerrits et al., 2014;Schlesinger and Jasechko, 2014;Wei et al., 2017). At catchment and watershed scales, an understanding of landscape functioning in turn helps designing robust models to predict the impact of climate extremes and environmental changes in society-relevant issues such as water resources management, flood forecasting, and impact assessment of land cover-land use change (e.g. Troy et al., 2015;Zhang et al., 2017).
Tracers have been of particular importance in understanding catchment functioning, as they highlight pore velocities of water molecules (i.e. how fast does a given parcel of water move) in a way that distinguishes this from the celerity (i.e. how fast energy propagates via the hydraulic gradient) of the rainfall-runoff response (McDonnell and Beven, 2014). Historically, isotopic transport models were initially developed at the plot scale (∼ 1-100 m 2 ) to represent 1-D isotope transfers in the soil profile and at the surface-atmosphere interface (Mathieu and Bariac, 1996;Melayah et al., 1996;Braud et al., 2005Braud et al., , 2009Haverd and Cuntz, 2010;Soderberg et al., 2012;Sprenger et al., 2018). Process-based simulation of isotopic trajectories has also been considered in larger-scale studies using land surface models (Haverd et al., 2011;Henderson-Sellers, 2006) where couplings with atmospheric isotopic circulation can be captured (Haese et al., 2013;Risi et al., 2016;Wong et al., 2017). While the simulation of energy budgets and biogeochemical cycles is increasingly detailed in these land surface models -sometimes including vegetation dynamics -the hydrology has, however, remained somewhat simplistic (or even absent) regarding the simulation of lateral transfers as overland flow, shallow and deeper subsurface flows, and channel routing (Fan, 2015). This makes it difficult to take advantage of isotope tracking to characterise the role of cascading downstream water redistribution in the spatial patterns of catchment functioning.
In parallel, isotopes have been used to explore water velocities, travel times, and ages in catchments using analytical and conceptual models (e.g. Neal et al., 1988;Barnes and Bonell, 1996;Weiler et al., 2003;Sayama and McDonnell, 2009;McGuire and McDonnell, 2015). These numerical tools allow testing hypotheses regarding how catchment storage relates to hydrological fluxes via mixing (or the relative absence thereof), and extending insights to spatiotemporal scales and variables inaccessible to current observation methods. An example of the latter is the estimation of water age, for which such models hold great promise McGuire and McDonnell, 2006;Sayama and McDonnell, 2009), with a more recent focus on the statistical properties of water transit time with timevarying and/or spatially distributed conceptualisations (Botter et al., 2010;Birkel et al., 2012;Heidbüchel et al., 2012;Harman, 2015;Rinaldo et al., 2015;Benettin et al., 2017;Heße et al., 2017). Additionally, the distinct information content of tracer observations, compared to more traditional hydrometric data, dictates that the integration of the two offers a strong hypothesis-testing framework for catchment model development (Uhlenbrook and Sieber, 2005;Fenicia et al., 2008;McDonnell and Beven, 2014). This opportunity is reinforced by decreasing costs of stable isotope analysis, now allowing for collection of daily (or more frequent) time series over several years (Kirchner and Neal, 2013) to inform simulations.
Yet, applications of a velocity-celerity framework in model-data fusion for catchment-scale hydrology remain relatively rare . Such studies are urgently needed at this scale where the emphasis is mainly on the characterisation of water pathways from precipitation to streamflow generation and/or evaporative losses. Recent efforts have nonetheless provided insights, either into wholecatchment dynamics with conceptual rainfall-runoff models (Birkel et al., 2011;Stadnyk et al., 2013;Hrachowitz et al., 2013;van Huijgevoort et al., 2016;Smith et al., 2016;Ala-aho et al., 2017;Knighton et al., 2017), or at finer detail  and using process-based 2-D hillslope models (Windhorst et al., 2014). We argue that extending tracer-aided approaches to physically based models could resolve both intra-and whole-catchment dynamics of stable water isotopes and bridge perspectives at multiple and process-specific scales, as recently shown in hydrometric-based studies (e.g. Endrizzi et al., 2014;Pierini et al., 2014;Niu and Phanikumar, 2015;Manoli et al., 2017). This process-oriented characterisation could also include non-conservative isotope behaviour such as evaporative fractionation, whereby lighter isotopes ( 1 H and 16 O) preferentially evaporate (Gat, 1996), and whose impact on downstream water signatures has been highlighted even in energy-limited landscapes (Sprenger et al., 2017a). Birkel et al. (2014) and Knighton et al. (2017) are amongst the rare attempts to include fractionation in catchment-scale studies, albeit with conceptual rainfall-runoff models. Investigation of internal catchment heterogeneity, marked in some geographical settings , is facilitated by spatially distributed resolutions of the catchment domain. In previous tracer-aided catchment modelling, however, this aspect is either indirectly considered -e.g. a semi-distributed separation of non-saturated/saturated domains  -or simply absent. Where spatial distribution has been taken into account in the model structure (van Huijgevoort et al., 2016;Ala-aho et al., 2017), fractionation processes were not included.
The two former pathways can result in evaporative fractionation, and root uptake for transpiration is usually considered non-fractionating (e.g. Wershaw et al., 1966;Dawson and Ehleringer, 1991;Harwood et al., 1999), although whether this is the case has recently been subject of debate (Lin and da SL Sternberg, 1993;Zhao et al., 2016;Vargas et al., 2017). While these different isotopic dynamics are of key importance in disentangling ecohydrological couplings in tracer-aided modelling, previous approaches generally lack a process-based conceptualisation of vegetation. Knighton et al. (2017) separately distinguished E t from other E components in catchment-wide isotopic model-data fusion. However, their spatially lumped approach was parsimonious, using empirical partitioning of potential evapotranspiration which has high uncertainty in natural ecosystems (Kool et al., 2014).
Here, we implement isotope and age tracking in the physically based, fully distributed model EcH 2 O (Maneta and Silverman, 2013). This model was chosen because it provides a physically based, yet computationally efficient representation of energy-water-ecosystem couplings where intracatchment connectivity (both vertical and lateral) is explicitly resolved. In addition, EcH 2 O separately solves the energy balance at the top of the canopy and at the soil surface, allowing a process-based separation of E s , E t , and E c . The novel isotopic and age-tracking module is designed in a manner directly consistent with the original model structure, assuming full mixing in each model compartment, and crucially without catchment-specific parameterisation. The conceptualisation of evaporation fractionation uses the well-known Craig-Gordon approach (Craig and Gordon, 1965).
We ask the following questions: -To what extent can a hydrometrically calibrated, physically based hydrologic model correctly reproduce internal catchment dynamics of isotopes?
-What are the limitations of these isotopic simulations? Do they relate to the underlying model physics and/or to the tracking approach adopted?
-How useful and transferrable is this model framework for simulating spatiotemporal patterns of isotopes and water ages?
These questions are here addressed by testing the new tracer-enhanced model (EcH 2 O-iso, Sect. 2) in a small, lowenergy montane catchment (Sect. 3). This site has previously been modelled applying EcH 2 O for calibration, using multiple datasets of long-term ecohydrological fluxes and storage variables (Kuppel et al., 2018a). We take advantage of this earlier work as a reference ensemble of calibrated model parameterisations, and no additional isotopic calibration is conducted. In addition to using long-term, high-resolution isotopic datasets for rainfall and runoff ( 2 H and 18 O), we assess the spatiotemporal variations of model-data agreement in soil water, groundwater, and plant xylem at different locations (Sect. 4.1). Following this generic evaluation, the model is used to infer seasonally varying patterns of water fluxes and isotopic signatures (Sect. 4.2), and water ages (Sect. 4.3). Model strengths and weaknesses, insights in processes, and potential ways forward are discussed in Sect. 5, before drawing conclusions in Sect. 6.
2 Model description

The EcH 2 O model
The ecohydrological model EcH 2 O combines a land surface module for calculating vertical energy balances (canopy and understory), with a kinematic-wave-based scheme for lateral and vertical water transfers, while vegetation productivity, allocation, and growth are derived from plant transpiration (Maneta and Silverman, 2013). Energy fluxes, water fluxes and storage, and vegetation state are explicitly coupled to capture the feedbacks between ecosystem productivity, hydrology, and local climate, at time steps larger than or equal to that of the meteorological inputs (precipitation P , incoming longwave and solar radiation, air temperature T a (maximum, average, and minimum), relative humidity, and wind speed). In addition, the flexible definition of the spatial domain in EcH 2 O allows for applications at a range of scales: from the plot (Maneta and Silverman, 2013), to small catchments (1-10 km 2 ; Lozano-Parra et al., 2014;Kuppel et al., 2018a), to larger watersheds (10 2 -10 3 km 2 ; Maneta and Silverman, 2013;Simeone, 2018). Despite some potential limitations due to the absence of diffusion-driven water redistribution or an explicit biogeochemical cycle providing ecosystem respiration, the model yielded satisfactory results and insights across the diversity of climatic settings (semiarid to humid) and scientific foci (e.g. water balance, energy balance, or plant hydraulics) covered by the aforementioned studies. A comprehensive description of EcH 2 O can be found in Maneta and Silverman (2013), and subsequent developments in Lozano-Parra et al. (2014) and Kuppel et al. (2018a).
We provide here a brief step-wise overview, focused on the different hydrological compartments and transfers simulated in EcH 2 O at the grid cell level (Fig. 1). For each vegetation cover present in a grid cell, a linear bucket approach is used for canopy interception. The capacity-excess P (i.e. belowcanopy throughfall) and P over bare soil are partitioned between liquid and snow components using a snow-rain temperature threshold (fixed to 2 • C) together with the minimum and maximum air temperatures at each time step. The canopy energy balance then separately yields plant transpiration (E t ) and evaporation of intercepted water (E c  . Water compartments (black rectangles) and fluxes (coloured arrows) as represented in EcH 2 O-iso, with the dashed arrows indicating processes where isotopic fractionation is simulated. The numbers between brackets reflect the sequence of calculation within a time step. Note that water routing (steps 8 to 13) differs between cells where a stream is present ( ) or not ( * ). of E t uses, for each vegetation type, the canopy conductance at each time step based on a Jarvis-type multiplicative model accounting for environmental limitations from incoming solar radiation, T a , vapour pressure deficit at the leaf surface, and soil water potential (see Silverman, 2013 andappendices in Kuppel et al., 2018a for a more detailed description). Infiltration of surface water in the topsoil layer is computed using a Green-Ampt approximation of the Richards equation. Subsequent soil water content above field capacity (gravitational water) percolates to the underlying soil layers, with fixed bedrock seepage -out of the systemas a lower boundary condition (Fig. 1). Soil evaporation (limited to the topsoil layer) and snowmelt (resulting in surface ponding; Fig. 1) under each vegetation type are calculated by solving the energy balance at the surface. Following a local drainage direction derived from the input elevation map, lateral water routing is simulated at three levels: in the deepest soil layer, groundwater seeps in the channel (if present) while the remainder is transferred laterally using a 1-D kinematic wave, and can result in saturation return flow in downstream cells. All remaining ponded water becomes overland flow, reinfiltrating further downstream or running off until it reaches an outlet or a cell within the stream network; stream water routing is also computed within a 1-D kinematic wave approximation.

Implementation of isotopic and age mixing
The conceptualisation of water mixing equally applies for all the tracked quantities implemented in the model (isotopes and age), so that a generic notation C is in this section used to designate both isotopic tracer composition ( 2 H and 18 O) and water age. The only specific conceptualisation of isotope dynamics in EcH 2 O-iso relates to fractionation (see Sect. 2.3), while precipitation inputs have a fixed age of zero and the water age in all compartments is increased at the end of each simulation time step by the length of the latter. The delta notation (δ) for isotopic composition quantifies, for a given water sample, the difference in the mass ratio of heavy to light isotopes (R) as compared to the Vienna Standard Mean Ocean Water (VSMOW): δ = R sample R VSMOW − 1 ×10 3 . First, the instantaneous mass balance for water signature is where V res and C res are, respectively, the volume and signature (δ 2 H, δ 18 O, or age) of the water in the reservoir, t is time, q out is the flux of water exiting the reservoir, and q in,k and C in,k are, respectively, the flux and signature of water entering the reservoir from each the N in adjacent upstream loca- tions. An implicit first-order finite-difference scheme is used to compute mixing during a given time interval t: where V t+ t res and C t+ t res are, respectively, the volume and water signature in the reservoir after mixing, V t res and C t res are the volume and water signature in the reservoir before mixing, and C t+ t in,k is the signature of the kth input source after mixing in the latter. Replacing V t+ t (2) finally yields In practice, Eq.
(3) is applied in EcH 2 O-iso at every sub-time step where water transfers are computed, in the sequence shown in Fig. 1. Note that C t+ t (2). Full mixing was used as a simplifying approximation because this model is to be first evaluated in a wet environment Sprenger et al., 2017a) with relatively long time steps (i.e. daily; see Sect. 3.3). Note that because of its representation of a single, fully mixed pool in each soil layer, EcH 2 O-iso essentially provides bulk water values for isotopic content and water age. This needs to be kept in mind when comparing the simulations with soil isotopic datasets (see Sects. 3.2 and 4) and for the discussion (Sect. 5).
One exception to immediate mixing is the snowpack, where the snowmelt flux (q melt out ) is assumed to tap first into the snow throughfall of the same time step (q snow in ) if present, before mobilising older snow, fully mixed in the snowpack. Consequently, the signatures of the snowpack (C t+ t pack ) and snowmelt water (C t+ t melt ) which goes into the surface reservoir in EcH 2 O-iso at step 7 ( Fig. 1) are calculated as follows: Only the same-time-step precipitation can contribute to throughfall in the EcH 2 O model, whenever the resulting canopy storage exceeds the maximum canopy storage capacity (Maneta and Silverman, 2013), the latter being constant in our simulations. As a result, all intercepted water eventually evaporates from the canopy and does not interact with the surface/subsurface. Therefore, throughfall water (liquid and snow) is assumed to have the isotopic composition of sametime-step precipitation and age zero. This simplification is reasonable for our study site where vegetation interception has only a trivial effect on the isotopic partitioning of rainwater , yet further developments could be implemented for model application in different ecoclimatic settings.
Finally, transpiration is considered as a non-fractionating process. This is based on previous work (Wershaw et al., 1966;Dawson and Ehleringer, 1991;Harwood et al., 1999), and the fact that non-steady-state effects cancel out at the daily time step (Farquhar et al., 2007). However, this simple conceptualisation is increasingly questioned (Lin and da SL Sternberg, 1993;Zhao et al., 2016;Vargas et al., 2017), and the implications for our study will be discussed later.
Here, during the canopy energy balance (step 2 in Fig. 1), the signature of transpired water C E t is taken as the weighted sum of the signature in the three soil layers: where f L1 , f L2 , and f L3 are the respective fractions of roots in each layer, as described in Eq. (A8) in Kuppel et al. (2018a).

Isotopic fractionation from soil evaporation
The change in isotopic composition of the first soil layer during soil evaporation (step 6 in Fig. 1) is simulated using the Craig-Gordon model Craig and Gordon (1965); Gat (1995), without any empirical parameterisation specific to our study. In this section, generically refers to the standardised isotopic ratio in for either 2 H or 18 O. For each time step t, where f is the remaining fraction of water after evaporation (f = V t+ t soil L1 /V t soil L1 ), while δ * is the limiting isotopic composition given the local atmospheric conditions in ‰ (Gat and Levy, 1978) and m is the dimensionless enrichment slope (Welhan and Fritz, 1977;Allison and Leaney, 1982). Their formulation is generalised following Good et al. (2014): The different terms in Eqs. (8) and (9) are sequentially defined as follows: δ a is the stable isotope composition of the ambient air moisture in ‰, derived from that of the precipitation by assuming isotopic equilibrium (Gat, 1995;Gibson and Reid, 2014): ε + is a factor (in ‰) derived from the equilibrium fractionation α + of water between the liquid and vapour phases (Skrzypek et al., 2015): with α + taken as temperature dependent following Horita and Wesolowski (1994), here using the air temperature T a : ε k accounts for the diffusion-controlled fractionation in air (Craig and Gordon, 1965): where D i /D is the diffusivity ratio of the gaseous water molecules bearing an isotope i to that of lighter isotopic water. We use literature values given by Vogt (1976) n translates the dominant mode of transport of water molecule at the surface. We adopted a time-varying formulation taking into account soil water content θ (Braud et al., 2005;Mathieu and Bariac, 1996): where φ and θ r are, respectively, the soil porosity and residual water content. n increases from 0.5 in a saturated soil to 1 for a dry soil where diffusion is the dominant mode of transport.
h a is the relative humidity of the atmosphere (measured at the weather stations; see Sect. 3.2) after being normalised to the saturated vapour pressure e * at the soil surface (Gat, 1995): where the surface temperature T s is given at each time step from solving the surface energy balance equation (Maneta and Silverman, 2013).
-Finally, h s is the relative humidity of the air of the soil pores, following the formulation of soil evaporation flux E s in EcH 2 O (Eqs. 9-10 in Maneta and Silverman, 2013): where β is adjusted as a growing function of the volumetric water content θ , equal to 1 whenever θ is superior or equal to field capacity θ fc (Lee and Pielke, 1992): 3 Data and methods

Study site
Simulations were conducted for the Bruntland Burn (BB) catchment in the Scottish Highlands (57 • 8 N, 3 • 20 W) ( Fig. 2a and b). It is a small (3.2 km 2 ) headwater catchment of the Dee, a major Scottish river providing freshwater resources for ∼ 250 000 people in the Aberdeen urban area, having EU conservation designations, and hosting ecosystem services (e.g. Atlantic salmon fishery). Annual precipitation averages around 1000 mm, with a mild seasonal cycle (Fig. 2c). The water balance is energy-limited, given the northern latitude, with 400 mm of annual evapotranspiration with pronounced seasonality in daily losses: from 0.5 mm in winter to 4 mm in summer (Birkel et al., 2011). Mean annual temperature is 7 • C and no monthly averaged temperatures fall below 0 • C; the climate qualifies as temperate to boreal oceanic; less than 5 % of precipitation usually occurs as snowfall.
The topography of the BB reflects glacier retreat, with elevation ranging from 220 m a.s.l. in the wide valley bottom to 560 m a.s.l. above the steeper slopes ( Fig. 2a). Glacial drift deposits cover 60-70 % of the catchment bedrock (granite, schist, and other metasediments) and form the dominant soil parent material . Mostly saturated, these deposits are important reservoirs of groundwater, sustaining base flows in the stream and maintaining persistent wet conditions across the valley bottom . Thin Regosols (rankers) dominate the pedology of the catchment above 400 m a.s.l., where drift deposits are marginal (Fig. 2a). Freely draining shallow podzols (< 0.7 m deep) dominate steeper hillslopes, overlying moraines and marginal ice deposits. Finally, deep (> 1 m) soils with high organic matter content (Histosols: peat and gley) characterise the riparian area (Fig. 2a). The Histosols are saturated most of the time, so that rainfall events generate runoff mostly via rapid saturation overland flow, with a surface connectivity in the podzols limited to the wettest periods . Spatial patterns of land cover reflect these hydropedological units (Fig. 2b). Heather shrublands (Calluna vulgaris and Erica spp.) are the dominant cover over podzols Table 1. Local datasets used in this study, grouped by location and purpose: model evaluation ( ), model calibration ( ), and model inputs ( ). For soil isotopes, a and b , respectively, indicate suction-lysimetric sampling (2013) and direct equilibration from soil sampling (2015)(2016). Other notations: Srf -surface water, GW -groundwater, P -precipitation, SWC -soil water content, E t -transpiration, NR -net radiation, * -relative air humidity, precipitation, air temperature, and wind speed, • -collection campaign at 92 to 94 locations (see text).

Water isotopes
Water fluxes and stores Meteorology  and rankers. Such a land use results from red deer (Cervus elaphus) and sheep overgrazing, at the expense of naturally occurring Scots pine trees (Pinus sylvestris L.), which are now mostly found in the steep sections of the northern hillslopes and in the plantation areas neighbouring the stream outlet. Finally, grasses (Molinia caerulea) cover the riparian gley soils, while the peat is dominated by bog mosses (Sphagnum spp.).

Datasets
We used the wealth of diverse and often multi-year time series available at different locations in the BB catchment (Fig. 2a). These measurements capture numerous ecohydrological processes and observables, used either for model inputs or calibration/evaluation of simulations (Table 1). A brief description follows.

Isotopic measurements
At the catchment outlet, rainfall and stream water have been sampled daily for isotope analysis from June 2011 to the present, providing an isotope time series of unusual high frequency and longevity. Samples have been collected using an ISCO 3700 automatic sampler (Teledyne ISCO, Lincoln, USA). The auto-sampler bottles were emptied at fortnightly frequency or higher, while paraffin was added to each bottle to prevent evaporation. Stable isotope measurements in the soil fall into two categories, differing in the sampling method and time period. Between 2011 and 2013, soil water was extracted at 0.1, 0.3, and 0.5 m depths at four locations: peat, gley, and podzol sites (weekly), and Forest site A (fortnightly) (Fig. 2a). Since MacroRhizon suction lysimeters were used (Rhizosphere Research Products, Wageningen, the Netherlands) , isotopic characterisation represents the mobile water held under lower tensions (Sprenger et al., 2015). From September 2015 to August 2016, near-monthly soil water sampling was carried out at four locations (Forest sites A and B, and Heather sites A and B) using soil samples collected with a spade from four layers (0-0.05, 0.05-0.1, 0.1-0.15, and 0.15-0.2 m) with five replicates for each. Isotopic analysis followed on water extracted by the direct equilibration method (Wassenaar et al., 2008), thus fully accounting for bulk pore water, as described by Sprenger et al. (2017a). Conceptually, the lysimeters can be viewed as sampling the "fast domain" of soil water held under low tension, whilst direct equilibration characterises the "bulk" soil water which also includes the "slow domain" of water held under higher tensions.
Groundwater samples were collected monthly between August 2015 and September 2016, at four wells (> 1.6 m) covering a representative transect from the hillslope to valley bottom (Scheliga et al., 2017) encompassing the main hydropedological units: peat (two wells), gley, and podzol sites (Fig. 2a). Vegetation xylem water was collected between autumn 2015 and spring-summer 2016, using cryogenic extraction from Scots pine xylem cores at 1.5 m height (Forest sites A and B) and from heather twigs (Heather sites A and B) (Fig. 2a). Sampling was made at near-monthly resolution (n = 7) with five replicas for each extraction (Geris et al., 2017). We also used isotopic measurements from a synoptic sampling campaign conducted in the drainage network of pools and channels across the valley bottom of the Bruntland Burn on 20 February (92 locations) and 24 May (94 locations) of the year 2013, covering contrasting catchment wetness states. On those days, water was also sampled along the perennial stream network at 10 locations .
Air-tight vials were used to store all water samples, which were kept refrigerated until they were analysed. The soil samples were equilibrated and extracted water analysed within a week of collection (Sprenger et al., 2017a).

Hydrometric and meteorological data
Daily soil moisture data were derived from 15 min retrievals at four locations: three along the peat-gley-podzol transect presented in Sect. 3.2.1 , and in a Scots pine stand (Forest site B; Fig. 2a). Time domain reflectometry (TDR) soil moisture probes (Campbell Scientific, Inc. USA) were located at different depths (0.1, 0.3, and 0.5 m -only 0.1 and 0.3 m in the peat), and replicated ∼ 2 m apart. During two growing seasons, Scots pine transpiration was measured at Forest site A (July-September 2015) and at Forest site B (April-September 2016) (Fig. 2a), by installing Granier-type thermal dissipation sap flow sensors (Dynamax Inc., Houston, TX, USA) on 10 and 14 trees, respectively. Depending on its stem diameter (10 to 32 cm), each tree had two to four sensors. At the end of each study period, incremental wood core sampling in surrounding trees provided sapwood-area -tree-diameter relationships, used to derive stand-scale transpiration estimates (Wang et al., 2017a), which were then daily averaged. At the catchment outlet (Fig. 2a), 15 min stage height records (Odyssey capacitance probe, Christchurch, New Zealand) were used to generate daily discharge observations, with a rating curve previously calibrated for a stable stream section.
Finally, meteorological observations used as model inputs (P , T a , relative humidity, and wind speed) and for calibration (net radiation) were primarily daily averaged from 15 min measurements at the three meteorological stations installed at different landscape positions (valley bottom, bog, and hilltop; Fig. 2a) and operated from July 2014. Prior to that period, a square elevation inverse distance-weighted algorithm was applied to interpolate local precipitation values from five Scottish Environment Protection Agency (SEPA) rain gauges located within 10 km of the Bruntland Burn catchment (Birkel et al., 2011). Daily mean T a , relative humidity, and wind speed values were then available from the Centre for Environmental Data Analysis (CEDA) at the Balmoral station (∼ 5 km NW) (Met Office, 2017). The ERA-Interim climate reanalysis dataset (Dee et al., 2011) was used to retrieve daily minimum and maximum T a (prior to July 2014), and incoming solar and longwave radiation (whole study period). Finally, we applied altitudinal effects on P and T a were accounted for; we applied a 5.5 % increase of P every 100 m a.s.l. (Ala-aho et al., 2017), and a 0.6 • C decrease per 100 m a.s.l., from the moist adiabatic temperature lapse rate (Goody and Yung, 1995).

Model set-up and calibration
The methodology closely follows the approach detailed in Kuppel et al. (2018a). Here, we only provide a brief summary and highlight the modifications adopted for this study.
All simulations were performed on daily time steps, at a 100 × 100 m 2 resolution. This choice of coarser grid cellsfrom 30 × 30 m 2 in Kuppel et al. (2018a) -was made to decrease computation time while preserving reasonable spatial variability across the catchment. The simulation period extends from February 2013 to August 2016, for which a continuous record of daily δ 2 H and δ 18 O in precipitation input was available (see Sect. 3.2.1). For all simulations, a 3-year spin-up period was added using the first 3 years of isotopic and climatic model inputs, as preliminary sensitivity tests combined with visual inspection of simulated hydrometric and isotopic time series at the location used in this study (see Sect. 3.2) indicated it was sufficient to remove transient dynamics.
Based on the soil classes defined by the Hydrology of Soil Types (HOST), four hydropedological units were defined ( Fig. 2a)  to map soil hydrological properties in the modelled domain. Physical soil characteristics relating to the energy balance were considered uniform across the catchment, similar to Kuppel et al. (2018a) (see Table S1 therein). Land cover was divided into five classes, four of them vegetated: Scots pine, heather shrubs, peat moss, and grasslands. From extensive land use mapping (Fig. 2b), the cover fraction of each vegetation type was estimated by combining 1×1 m 2 resolution lidar canopy cover measurements , aerial imagery, and typical vegetation patterns in the different soil units Kuppel et al., 2018a). As in Kuppel et al. (2018a), the dynamic vegetation allocation module is switched off, so that leaf area index remains equal to initial values of 2.9, 1.6, 3.5, and 2 m 2 m −2 for Scots pines, heather shrubs, peat moss, and grasslands, respectively (Albrektson, 1984;Calder et al., 1984;Bond-Lamberty and Gower, 2007;Moors et al., 1998). To avoid an overestimation of local soil evaporation and resulting isotopic fractionation in grid cells of exposed rock/scree, for simplicity, we fixed the depth of the first soil layer to 0.001 m wherever the fraction of bare soil was larger than 0.5 -after performing a sensitivity analysis showing little effect on catchment water balance and downstream isotopic budgets.
Finally, the calibrated model parameters (Table S1 in the Supplement), and associated sampling ranges, are those presented in Kuppel et al. (2018a). The parameter space was sampled using a uniform Monte Carlo approach. The corresponding 150 000 simulations were jointly constrained combining measurements of stream discharge, soil moisture (four sites), pine E t (two sites), and net radiation (three sites) (Table 1) whenever the observation periods overlapped with the current simulations (February 2013-August 2016; Fig. S1 in the Supplement). For soil moisture observations, a b-spline curve was fitted to the measured profile (to account for non-monotonic variations) on each sampling date, followed by a vertical integration. It enabled a consistent comparison against simulations in each of the upper two hydrological layers of EcH 2 O (see Fig. 1), while profile-averaged values were used for calibration in Kuppel et al. (2018a). Constraints were combined in a multi-criteria objective function based on the cumulative distribution functions (CDFs) of dataset-specific goodness of fit (GOF) (Alaaho et al., 2017): mean absolute error for stream discharge and root mean square error for all others observations. This method allows retention of model parameter sets that give most behavioural simulations simultaneously across different variables (Kuppel et al., 2018a). We retained the 30 "best" of these parameterisations as a test bed for ensemble simulations of stable water isotopes and water age dynamics presented in this study.

Analysis
Daily, seasonal, and interannual climate variability result in changing isotopic composition of precipitation inputs. Equilibrium isotopic fractionation processes result in a strong correlation between rainfall δ 2 H and δ 18 O across the globe, defining a global meteoric water line (GMWL; Dansgaard, 1964). At the BB catchment, there is a seasonal trend of more enriched values in summer and depleted in winter (e.g. Fig. 2c). A local meteoric water line (LMWL) was defined, using daily values from February 2013 to August 2016 and weighting by precipitation inputs (r 2 = 0.96, p < 0.001): The line-conditioned excess (hereafter, lc-excess) was defined as the residual from the LMWL (Landwehr and Coplen, 2006): with a LMWL = 7.8, b LMWL = 4.9 ‰ (Eq. 20). As oxygen has a higher atomic weight, non-equilibrium fractionation during the liquid-to vapour-phase change will preferentially evaporate (in terms of statistical expectation) 1 H 2 H 16 O molecules rather than the heavy isotopologue 1 H 218 O (Craig et al., 1963). The isotopic signature of a water sample affected by evaporation thus shows negative lc-excess values, as δ 18 O in non-evaporated water enriches faster than δ 2 H, and plots under the LMWL in the dual-isotope space (Landwehr et al., 2014). For these reasons, we preferred combining δ 2 H and lc-excess in our analysis (over separately looking at both δ 2 H and δ 18 O), to simultaneously highlight absolute isotopic dynamics and evaporative fractionation. Note that lc-excess was also preferred over the oft-used deuterium-excess, which translates the deviation of δ 2 H from the GMWL (Dansgaard, 1964). While the two quantities are mathematically similar, lc-excess displays much smaller seasonal dynamics from the near-0 ‰ value of precipitation inputs; thus, it advantageously allows separation of fractionation impacts from overall isotope dynamics (Sprenger et al., 2017a). Similar to soil moisture observations, measured and simulated isotopic values in the soil are conceptually different: datasets are collected at specific depths (see Sect. 3.2.1), whilst model outputs provide average values for the different hydrological layers (Fig. 1). While original quantities were preserved for temporal analysis of the results, we additionally provided a formal quantification of model-data agreement. To do so, we reconstructed layer-integrated isotopic datasets at each soil sampling site, following the same interpolation-integration methodology used for soil moisture for computing model-data goodness of fit during calibration (Sect. 3.3).
As outlined in Sect. 1, our model evaluation is meant to test the ability of EcH 2 O-iso to generically simulate isotope dynamics across compartments. We used mean absolute error (MAE) to quantify model-data fit for all isotopic outputs, some of which present low temporal variability, have skewed distributions, or have a relatively lower sampling record, resulting in typical hydrograph-oriented efficiency metrics (e.g. Nash-Sutcliffe or Kling-Gupta;Nash and Sutcliffe, 1970;Kling et al., 2012) being less applicable. The median values are shown on corresponding time series (Figs. 3-7), and normalised by each dataset range and used in conjunction with Pearson's correlation factor in Fig. 8 as a summary of model performance. The correlation coefficient axis in this dual model performance space represents the quality of the model in representing the variation of the data, while the normalised MAE axis provides information on the accuracy (bias) of the model.

Time series
Seasonal dynamics of soil water isotopes were well captured on the hillslopes, as exemplified at two sites in Fig. 3: one located in the shrub-dominated moorland (podzol site), the other in a Scots pine plantation (Forest site B), noting that the graphs cover different hydrological years dictated by data availability. Model-data agreement was consistent for δ 2 H, keeping in mind that while measurements were depthspecific, simulated values were averaged over the first and second hydrological layers (Fig. 1). As a result of model calibration (see Sect. 3.3), the thickness of the first (topsoil) and second layers spans 0.10-0.19 and 0.02-0.39 m in the simulated podzol soil unit, respectively (not shown). Still, at the podzol site the model captured well the vertical variability δ 2 H across the summer of 2013 but overestimated topsoil enrichment during the following winter (Fig. 3a). Lc-excess was generally underestimated in the topsoil layer there, with negative simulated values indicating evaporative influence generally not found in the data. At Forest site B, both δ 2 H and lc-excess dynamics showed modelled ranges consistent with measurements. Note, however, that EcH 2 O-iso simulated a vertical profile during the winter of 2016 with richer δ 2 H in the deeper layer, a condition that was only occasionally found in δ 2 H measurements (November 2015 and January 2016). At both sites, the temporal dynamics of soil moisture were Geosci. Model Dev., 11, 3045-3069, 2018 www.geosci-model-dev.net/11/3045/2018/ well captured by the model (bottom rows). We note however that the observed decrease of moisture content with depth -especially marked at Forest site B -was generally not reproduced, as the vertically constant parameterisation of soil hydrology in EcH 2 O (Brooks-Corey conceptualisation; Maneta and Silverman, 2013) does not allow sufficient water retention in highly organic upper soil layers. Isotopic consistency was also found further in the valley bottom, as shown at the peat site in the riparian area (Fig. 4). The bimodal summer δ 2 H enrichment measured was well captured in the topsoil layer of the model (thickness: 0.02-0.25 m), as were the mildly negative lc-excess values. In addition, the weak variability and range of measured δ 2 H and lc-excess at greater depths were consistently reproduced. As for the podzol site, we noted in the peat higher peak enrichment values from the model than for the available data. As for other elements of the analysis, we remind here that soil isotopic data were sampled differently at the three sites: lysimeter extraction was used for the podzol and peat sites, therefore characterising mobile water in the fast domain, while the direct equilibration analysis conducted at Forest site B effectively applies to bulk water including water held under higher www.geosci-model-dev.net/11/3045/2018/ Geosci. Model Dev., 11, 3045-3069, 2018 tension. The model essentially gives a bulk isotopic composition of stored water (Sect. 2), which might also explain why results were comparatively better at Forest site B. We also explored the accuracy of simulated spatial patterns of isotopic signatures in the riparian zone, using two synoptic sampling surveys of surface water and stream water (Sect. 3.2.1) on separate days in late winter and late spring of 2013 ( Fig. 4d-g). A good agreement was found for δ 2 H in the main branch of the stream network on both dates, while there was a tendency to overestimate δ 2 H values in the northwest part of the riparian area. Model-data fit of lc-excess mostly oscillated between good and a few per mille underestimation, depending on the location. For both δ 2 H and lcexcess there was fine-scale spatial variability in the modeldata fit, especially marked in late May and in the main channel. This reflects both the spatial variability of measurements Sprenger et al., 2017b) and the different resolution of sampling (∼ 10 m intervals) and the much coarser grid of the simulations (100 × 100 m 2 ). Figure 5 shows EcH 2 O-iso's simulation of the isotopic imprint of plant water uptake in the transpiration flux. The isotopic composition of xylem water samples was directly compared to that of root water uptake simulated from the canopy energy balance (sub-step 2b in Fig. 1). At the heather sites, the simulated ranges were consistent from model to data, with an excellent model-data fit for lc-excess despite the low 90 % spread of simulations outputs ( Fig. 5a and b). The seasonal cycle of simulated δ 2 H conversely seemed opposite to that of xylem samples, which showed gradual enrichment in winter followed by depletion at beginning of the growing season, but the lack of data from January to April limits gen- eral seasonal interpretation. At the forest sites, simulation results were very similar, noting that Forest site A corresponds to the same model grid cell as Heather site A (Fig. 5c and d). However, the measured isotopic composition in xylem was quite different for Scots pine compared to the heather, in two ways. First, the seasonal trends of δ 2 H were reversed, resulting in a good agreement with the modelled seasonality. Second, measured δ 2 H and lc-excess showed consistently lower values as compared to the heather sites, by 5 ‰-24 ‰ for δ 2 H and 4 ‰-13 ‰ for lc-excess (δ 18 O was only slightly positively biased, not shown). As a consequence, simulations showed a permanent positive offset for Scots pine water use despite consistent seasonality.
Isotopic variability was comparatively much lower for the groundwater both in time and across monitored wells, and a general agreement was found in the simulations (Fig. 6). ble. However, the model tended to slightly underestimate lcexcess, with simulated values near zero while measurements were mostly centred on 3 ‰. In addition, the short-term lcexcess variability was somewhat underestimated in the riparian area. Figure 7 shows the model-data comparison at the catchment outlet. The overall signal of stream water δ 2 H (Fig. 7a) and discharge values (Fig. 7c) were well reproduced by the model, with consistent "transition" periods of progressive enrichment when atmospheric demand increased and the catchment got drier. Most behavioural models in the ensemble did not completely capture the full extent of winter δ 2 H depletion, and the seasonal minimum of δ 2 H generally fell below the 90 % spread of the ensemble. However, seasonal variations of modelled lc-excess in the stream were in phase with the datasets throughout the study period: minimum in summer and maximum in winter, although simulated variability was more damped than for δ 2 H, with a slight negative bias (Fig. 7b).
A summary of model performance is shown in Fig. 8 for all sites/compartments, using the dual space of normalised MAE (using each dataset range, x axis) and Pearson's linear correlation factor (y axis). The vast majority of median normalised MAE values were below 1, and more than half of evaluated datasets showed values below 0.5. Values above 0.7 were mostly found for groundwater and xylem compartments, a clustering especially marked for δ 2 H. In addition, most median model-data correlations were significantly positive between 0.4 and 0.85, noting a tighter clustering around high values for δ 2 H than lc-excess. Insignificant or negative correlations were mostly found where only a few data points were available (xylem) or where seasonal variability was low (e.g. groundwater). Interestingly, median model-data agreement in topsoil at Forest site A significantly differed between 2013 (mobile water sampling via lysimeters) and the 2015-2016 period (bulk water sampling via direct equilibration). This was notable in the dramatic increase of model-data correlation (0.17 to 0.8) and decrease of normalised MAE (0.5 to 0.25) for topsoil δ 2 H in the latter case, which is consistent with our interpretation that the simulated soil water composition represents that of bulk water. Figure 9 provides a spatially distributed, seasonal view of the ensemble median of outgoing water fluxes across the catchment over the simulation period. Lateral connectivity was markedly higher during the wetter first half of the hydrological year (October-March; Fig. 9a and b). During this colder, most energy-limited period, surface runoff -cumulative along the flow path, as runoff can cross several grid cells within one time step -was significant in many cells where the slopes transition to the valley bottom (up to 53 mm day −1 , i.e. 0.006 m 3 s −1 ), as well as on some surrounding hillslopes in the southern/southwestern part of the catchment (Fig. 9a). Throughout the spring-summer, very few of these overland flow corridors were usually hydrologically active, only in response to larger storm events. In parallel, lateral subsurface connectivity in autumn-wintertime was quite widespread across the catchment, particularly concerning the two southernmost stream tributaries where subsurface flux largely exceeded surface runoff (up to 90 mm day −1 ; Fig. 9b). Some of these subsurface connections were still active during the growing season, albeit weaker (< 40 mm day −1 ). Given the predominance of subsurface flow near the channel, return flow dominated the vertical water budget (exfiltration minus infiltration > 0) throughout the year at junctions with the main stream and further downstream, especially in the winter (Fig. 9c). The rest of the catchment was dominated by infiltration, with average net rates of a few mm day −1 . Evaporative losses of soil water were much smaller and had a different seasonality than infiltration and throughflow ( Fig. 9d  and e). In autumn-winter, soil evaporation (E s ) was similar in magnitude to ecosystem transpiration E t (integrated over all vegetation cover for each grid cell), although at local scales both fluxes remained below a few tenths of mm day −1 (catchment average: 0.11 mm day −1 for both E s and E t ). Conversely, ecosystem transpiration clearly dominated during the rest of hydrological year, with a catchmentaveraged rate almost 4 times higher than that of soil evap-oration (0.61 mm day −1 vs. 0.16 mm day −1 , respectively). In both cases, the highest values were found in the riparian area, although the spatial contrast was more marked for soil evaporation.

Simulated hydrometric and isotopic spatial patterns
This spatiotemporal variability in water fluxes was somewhat reflected in that of isotopic patterns (δ 2 H in Fig. 10a, and lc-excess in Fig. S2). δ 2 H in the topsoil went from markedly depleted winter values (average: −61 ‰) to maximum enrichment in spring-summer with larger spatial variability (average: −44 ‰) (Fig. 10a). These temporal variations were well within that of δ 2 H in precipitation inputs (Fig. 2c). Yet, the increasing spatial variability of topsoil δ 2 H in spring-summer, and the much more pronounced relative seasonality of topsoil lc-excess (Fig. S2a) (compared to that in precipitation; Fig. 2c), indicated a significant influence of evaporation fractionation on isotopic patterns. During the spring-summer period the highest δ 2 H values, and most negative lc-excess values, were found in the organic soils of the valley bottom and on the higher hillslopes. These locations are where soil evaporation was highest (Fig. 9d) or where the soils are thinnest (rankers Regosols; Fig. 2a). The effect of isotopic fractionation crucially depends on relative storage change (Eq. 7); thus, it had large values either because absolute evaporation was high (valley bottom) or because the available storage was limited (thin soils). Conversely, springsummer lc-excess values were near zero (or even slightly positive), and δ 2 H enrichment less pronounced, in most of the topsoil grid cells where the stream is also present, corresponding to the locations where upslope-routed groundwater exfiltrated (Fig. 9c). Finally, positive winter values for lcexcess across the catchment's topsoil hint at a widespread dominance of winter precipitation and mixing processes (via surface connectivity and infiltration; Fig. 9), over fractionating ones.
The isotopic signature (δ 2 H and lc-excess) in water used by plants for transpiration largely displayed a damped reflection of the topsoil patterns (Figs. 10a and S2a). This reflects distributed root uptake across the soil profile, reaching deeper soil compartments where seasonal isotope dynamics were less marked. One consequence is that the model simulated more isotopically depleted plant water use in the thin Regosols of the upper hillslopes compared to the very shallow topsoil layer (northern and western parts of the catchment).
Finally, groundwater δ 2 H patterns were comparatively more uniform across the catchment (σ spatial = 1.9 ‰) and across seasons (Fig. 10c). Most depleted values were found in the podzolic hillslopes and across the valley bottom, a feature more marked in winter and spring. Lc-excess mostly displayed positive values throughout the year, except for some weakly negative autumn values on the higher hillslopes (Fig. S2c). Markedly positive values were generally found in the organic soil of the valley bottom where fluxes converge. Note that positive values were more spatially homogeneous during winter and springtime, highlighting subsurface Geosci. Model Dev., 11, 3045-3069, 2018 www.geosci-model-dev.net/11/3045/2018/ recharge lagging behind the more superficial compartments by a few months.

Water ages
Simulated water ages showed significant variability across locations in the catchment, as well as a marked seasonality at most sites selected for the analysis (Fig. 11). For convenience, the sites chosen for analysis in Fig. 11a-d were the same as those where isotopic model-data evaluation was conducted. In general, modelled water age increased with distance downhill, consistent with freely draining hillslopes sustaining groundwater fluxes into the riparian area. In the soils, water age ranged from a few weeks on the hillslopes to several years in the valley bottom peat where the topsoil is affected by exfiltration of older groundwater from upslope areas. Groundwater age was more homogeneous across the watershed but still showed significant differences, averaging 1 year in the podzol-covered locations, compared to 2 to 3 years in the riparian area. Seasonal variations were most significant on the hillslopes, from week-old waters in winter to water ages of 2 to 6 months during the growing season in the vadose zone. Weaker intrinsic seasonal variability was generally found in groundwater, which is consistent with the very flat simulated isotope dynamics (Fig. 6). The age of water taken up by plants followed the topsoil age patterns in most cases, reflecting the relatively young water ages from shallow rooting depths. One exception is Forest site A, where the contribution of older water from the second soil layer during the growing season had a clear effect on the age of the water used by vegetation. This latter site interestingly dis- played older water ages compared to other hillslopes locations, suggesting slower drainage conditions, likely linked to less marked local topography and receipt of older water from upslope. In addition, the gley site displayed rather dynamic behaviour in the upper soil layers, similar to podzols, while the turnover of groundwater there was the lowest among all locations, suggesting it is confined and disconnected from the soil profile. Finally, at the peat site, younger water ages were found in groundwater compared to upper soil layers. This surprising result was likely linked to permanently saturated soils with limited infiltration (disconnection from the surface and overland flow) and recharge (disconnection from confined groundwater) where lateral soil water movement was not simulated, but lateral transfers and mixing occurred in the underlying groundwater.
Geosci. Model Dev., 11, 3045-3069, 2018 www.geosci-model-dev.net/11/3045/2018/ Spatial variability was also found in stream water age, as shown in previously referred-to sites and arbitrarily defined locations along the channel network (Fig. 11e). In two of the main tributaries of the BB (HW2 and HW3), simulated water ages were significantly younger (∼ 0.5-1.5 years) than along the rest of the channel network (1.4-2.8 years). The older water ages found in HW1 might be linked to the presence of high water storage in drift deposits and an extensive raised peat bog in this portion of the valley bottom (Sprenger et al., 2017b), while the streams in HW2 and HW3 emerge further upslope at the drift-free ranker-podzol transition (Fig. 2a). There is a localised increase in water age when moving downstream towards the peat site, consistent with increased groundwater exfiltration (Fig. 9c) where stream water is a few weeks older than at the catchment outlet. In this lower part of the catchment lower temporal variability was also evident (1.8-2.4 years). Again, this might be derived from groundwater influxes and the extensive presence of saturated peat soils in this part of the catchment, compared to other sections of the stream.

Performance of the tracer-enhanced model
The model-data comparison demonstrated that EcH 2 O-iso captured a significant part of the isotopic behaviour across multiple ecohydrological compartments and landscape positions monitored in the study catchment. Because no calibration was performed on the isotopic components, these results reveal that the water mixing and storage and the water pathways simulated by the hydrologic core of EcH 2 Oiso correctly reflect the dominant hydrologic dynamics of the basin (Kuppel et al., 2018a). Hydrological states and fluxes in the model evolve driven by the celerity of propagation of local energy gradients (e.g. gravity-driven hydraulic gradient) throughout the landscape, with no direct knowledge about which "water parcels" (e.g. old or young, upslope or downslope) have been mobilised during a given hydrological response (Kirchner, 2003). Conversely, correctly capturing isotope dynamics is conditioned to accurately simulating patterns of water particle velocities, i.e. to routing the correct water parcels all the way from precipitation to their fate in the stream or as evaporative outputs (McDonnell and Beven, 2014). Therefore, the general performance achieved in the present celerity-velocity framework gives reasonable confidence in the mechanistic description of energy-water-plant couplings adopted by the EcH 2 O-iso model.
Despite some of the discrepancies presented and discussed below, the overall isotopic model-observation fit is very encouraging because the evaluated ensemble of model configurations was not derived from any tracer-aided calibration but solely used the information content brought by hydrometric and energy balance datasets in an independent cali-bration exercise similar to Kuppel et al. (2018a). Further, the implementation of water isotope and age tracking, consistent with the original structure of EcH 2 O and including evaporation fractionation of isotopes, was straightforward and followed well-established methodologies  without any parameterisation specific to the study site. By keeping both the isotopic module and calibration as minimalistic as possible, our approach avoids adding new, unnecessary degrees of freedom and reduces the risk of overfitting. Specific model performance might thus be lower than what could be achieved using a dual hydrometric-isotope calibration approach van Huijgevoort et al., 2016;Knighton et al., 2017), but because the isotopic signal remains truly independent of the hydrologic calibration our approach allows unique critical analysis and insight into the physical hypotheses underlying simulated flow generation and water mixing.

Insights into critical processes for model future development
The timing of seasonal isotopic dynamics as well as higherfrequency responses was well simulated in the vast majority of cases (summary in Fig. 8), together with value ranges also broadly consistent with observations. Yet, the amplitudes of modelled temporal isotopic responses displayed variable degrees of agreement with that of measured signals. In general, dynamics of deuterium were better reproduced than those of lc-excess, with a trend to underestimate lc-excess in several compartments. One of these model-data mismatches is the overly enriched signal in the topsoil of the riparian sites during the growing season (Figs. 4 and 8). The concomitant underestimated lc-excess hints at an excessive evaporation fractionation signal. As pointed out in Sect. 4.1, these discrepancies can partially derive from the different information represented by model and observations. For instance, the model simulates the composition of the bulk topsoil water, whereas the observations may reflect only the composition of the free draining portion of the soil water. At the long-term riparian locations (peat and gley sites; Fig. 2a), collection by suction lysimeters  was used, sampling water under low tension and less affected by fractionation (Brooks et al., 2010;Sprenger et al., 2017b). In addition, the samples from a synoptic field campaign across the extended riparian area in flowing surface waters and ponds on two dates  were directly compared to simulated topsoil water ( Fig. 4d-g). This was because the current formulation of EcH 2 O-iso routes all surface water to the next downstream cell and thus does not account for free-standing water such as the ponds and zero-order streams typically forming outside summer in the BB, particularly in the northwest of the catchment ). Yet, the sampled surface water in the riparian area has been shown to have spatially varying sources, presenting distinctively enriched or depleted δ 2 H signals depending if the source is soil water or groundwater, respectively . Systematically comparing soil water to sampled surface water might thus explain the overestimation of δ 2 H, especially in the northwest part of the catchment where limited groundwater seepage is modelled (Fig. 9c). Secondly, the riparian topsoil in EcH 2 O-iso might function as an "evaporation hotspot" to a greater extent than has been found in corresponding sampled surface water locations (Sprenger et al., 2017b). Indeed, topsoil water is not laterally connected in the model, so that evaporation fractionation remains local (horizontally) but immediately mixes across the whole layer -as compared to a vertically stratified isotopic profile in poorly mixed ponded areas. In addition, while fractionation is modest compared to other climatic settings, measurements have shown that ponds and zero-order channels that are not fully evaporated connect to the channel network in spring-summer and drive the seasonal isotopic enrichment observed in the stream (Sprenger et al., 2017b). Further support for this hypothesis was found by disabling evaporative fractionation in our simulations: seasonal variability of isotopes in the stream almost completely disappeared, while short-term, event-driven dynamics remained (not shown). Beyond the idiosyncrasies of our study catchment, and the gap between fine-scale wetland heterogeneity and our model resolution (100 × 100 m 2 ), a large body of literature has reported the importance of riparian wetlands as time-varying "chemostats" controlling stream water quality (e.g. Billett and Cresser, 1992;Smart et al., 2001;Spence and Woo, 2003) or "isostats" mixing isotope signals . Since modelled soil water is not laterally routed to the channel during the onset of the growing season, this might explain some underestimation of summer δ 2 H in the stream outlet, as well as the reported lack of seasonal variability for instream lc-excess (Fig. 7). Further developments of the model to include ponding effects and/or a more dynamic channel network (rather than fixed, as currently conceptualised) would help capture these seasonally varying flow paths in the variably saturated valley bottom of low-energy landscapes.
The isotope and age tracking adopts a complete and instantaneous mixing scheme at each sub-time step where water transfers are computed between the spatially distributed compartments of the simulated domain. This working hypothesis was chosen for simplicity, given the wet and cool climate conditions and the relatively long (daily) simulation time steps. The spatiotemporal variability of simulated fluxes and stores somewhat results in a time-variant partial mixing at the catchment scale at the stream outlet (van Huijgevoort et al., 2016). However, we note, for example, that our simulations of groundwater lc-excess showed an underestimated variability and a consistent negative bias towards near-zero values (Fig. 6). It indicates that the simulated recharge signal is very damped throughout the year and slightly biased towards the signature of over-enriched, evaporation-affected recharge. This contrasts with the evidenced dominance of winter recharge given the markedly positive lc-excess values observed at the monitored wells (Scheliga et al., 2017) as well as in other catchments with comparable ecoclimatic settings (O'Driscoll et al., 2005;Yeh et al., 2011;Bertrand et al., 2014). It might point to an exaggerated mixing across the soil profile in our simulations, overly flattening the precipitation signature and overestimating fractionation signal in the water percolating to the water table. Given that groundwater directly sustains 19 (±16) % of annual stream flow in our ensemble simulations (not shown), one can link this lack of variability in groundwater lc-excess to that simulated in stream water (Fig. 7). While such a link between the degree of unsaturated zone mixing and stream isotopes was not evidenced by Knighton et al. (2017), there was a much lower contribution of baseflow to discharge in the intermittent catchment they modelled. More generally, further developments would benefit from incorporating insights from the growing body of literature on the importance of preferential flow in driving catchment dynamics and tracer mixing (Beven and Germann, 2013). This would first involve implementing conceptualisation of microtopographic controls on overland flow (Frei et al., 2010). Secondly, the significance of subsurface dual pore space (matrix-macropore) representations of tracer flow paths and mixing has long been put forward (Beven and Germann, 1982) but modelling efforts relevant to catchment hydrology remain somewhat scarce (Stumpp et al., 2007;Stumpp and Maloszewski, 2010;Vogel et al., 2010;Sprenger et al., 2018;Smith et al., 2018). Bridging these detailed plot-to hillslope-scale descriptions with a physically based ecohydrological model such as EcH 2 Oiso will likely require a simplified, parsimoniously parameterised implementation and calibration with tracer data.
Our modelling experiment also helps to evaluate the conceptualisation of isotopic fractionation in the soil water of wet, energy-limited catchments. The evaporative fractionation is described by the well-established Craig-Gordon model (Craig and Gordon, 1965), supplemented here with a soil-adapted formulation following Mathieu and Bariac (1996) and Good et al. (2014). As reviewed by Horita et al. (2008), the Craig-Gordon model is very sensitive to the isotopic composition of atmospheric moisture (δ a ), the relative humidity of the atmosphere at the surface (h a ), and the kinetic fractionation factor ( k ). We assumed isotopic equilibrium between rainfall and atmospheric moisture (Eq. 10), as is commonly done when no direct measurement of δ a is available (Horita et al., 2008). While this empirical, and here spatially uniform, approach is valid on monthly timescales in temperate climates (Schoch-Fischer et al., 1983;Jacob and Sonntag, 1991), discrepancies can arise on shorter timescales and/or when local evaporation significantly feeds atmospheric moisture (Krabbenhoft et al., 1990). Second, h a estimates can be a large source of error in wet environments where h a > 0.75 (Kumar and Nachiappan, 1999), which is often the case in our catchment (Wang et al., 2017b). Furthermore, we found a marked sensitivity of isotope dynamics Geosci. Model Dev., 11, 3045-3069, 2018 www.geosci-model-dev.net/11/3045/2018/ to the strategy used to calculate k (Eq. 14), consistent with Haese et al. (2013), who found a large impact on simulated soil δ 18 O in northern latitudes. We chose to use a formulation based on isotopic diffusivity ratios; the latter were taken from Vogt (1976) because their experimental protocol covered a comparatively large range of humidity conditions. Yet it seems that very few (if any) experimental studies estimating these ratios spanned the very humid conditions found at the BB, and further empirical data could help reduce the associated uncertainties (Horita et al., 2008). Finally, we showed that our root uptake simulations for heather shrubs broadly matched the measured isotopic signature in plant xylem. Conversely, a systematic, positive model offset was found for both δ 2 H and lc-excess in Scots pines despite the fact that the model correctly captured the temporal dynamics (Fig. 5). Our simulations assumed identical, exponential root profiles for all vegetation types within soil types, e.g. the podzol, where these experimental heather and forest sites are found (Kuppel et al., 2018a); thus, speciesdependent use of soil water from depth-specific isotopic signature cannot be captured. Heather shrubs have, however, a shallow root system (typically < 5 cm; Geris et al., 2017), and thus its source water might be more affected by evaporation than Scots pine (which can be deeper rooted). However, the observed lc-excess values in the soil of Scots pine (−13 to 5.5 ‰; Fig. 3b, Forest site A not shown) were significantly higher than those measured in the pine xylem (−19.6 to −7.6 ‰; Fig. 5c and d). It mostly seems to stem from significant recorded deuterium depletion while 18 O ratios were consistent or slightly depleted as compared to soil samples, and we found larger simulation biases (relative to the mean value) in xylem for deuterium than for 18 O ratios (not shown). Such isotopic departures between soil and xylem water have been reported in a number of experimental studies, although primarily conducted in seasonally drier environments (Lin and da SL Sternberg, 1993;Zhao et al., 2016;Vargas et al., 2017). Several mechanisms have been proposed, including a discrimination of heavier isotopes during water uptake controlled by root aquaporins (Mamonov et al., 2007) or mycorrhizal associations (Berry et al., 2017), phloem-xylem water cycling on several timescales (Hölttä et al., 2006;De Schepper and Steppe, 2010;Pfautsch et al., 2015;Stanfield et al., 2017), and stem water evaporation through the bark (Dawson and Ehleringer, 1993). While exploring the relevance of these mechanisms to the ecosystems here simulated goes far beyond the scope of this study, it is clear that the complexity of isotopic dynamics in plant xylem cannot be fully captured simply based on a root-profile-weighted mixing of soil pools.

Opportunities for characterising water pathways
The development of EcH 2 O-iso is a methodological "middle path" for modelling conservative tracer transport, between detailed plot-scale models across the soil-vegetationatmosphere continuum (e.g. Mathieu and Bariac, 1996;Melayah et al., 1996;Braud et al., 2005;Haverd and Cuntz, 2010), catchment rainfall-runoff models McGuire and McDonnell, 2015;van Huijgevoort et al., 2016;Knighton et al., 2017), and land surface models for Earth system studies (Haese et al., 2013;Risi et al., 2016;Wong et al., 2017). This reflects the reasons why the original EcH 2 O model was developed, namely to provide a physically based, yet computationally efficient representation of energy-water-ecosystem couplings where intra-catchment connectivity (both vertical and lateral) could be explicitly resolved (Maneta and Silverman, 2013). The combination of these features is critical, since explicit lateral connectivity (surface, subsurface, and channel) is typically the missing piece in land surface models (Fan, 2015) and in plot-scale approaches, and the coupling with vegetation processes is typically missing in rainfall-runoff models (van Huijgevoort et al., 2016). The newly developed model provides, for the first time, a transferable, process-based linkage of spatial-temporal patterns of water fluxes ( Fig. 9) with those of isotopic tracers (Figs. 10 and S2) across a headwater catchment.
Here, a major focus has been put on the isotopic analysis to evaluate the consistency of EcH 2 O-iso using the wealth of data available at the study site, and the limitations stemming from the unavoidable technical trade-off we adopted. Yet, principles used for isotope tracking were applied to track water age across the ecohydrological compartments (Fig. 11). This provides a more complete picture of catchment functioning than stream water age, although the latter metric provides an important first-order benchmark for comparison with other modelling approaches. The mean stream water age of ∼ 2.1 years is consistent with isotope-calibrated rainfall-runoff approaches reporting ∼ 1.55 years (van Huijgevoort et al., 2016;Ala-aho et al., 2017) and ∼ 1.8 years . The low temporal variability found here yields higher discrepancies when considering flowweighted median ages: ∼ 2 years against 1.2 years found by Soulsby et al. (2015) and ∼ 1 year reported using transport model driven by StorAge Selection functions (Benettin et al., 2017). We notably find a slower water turnover in the valley bottom soils (∼ 2.8 years) as compared to compared to the spatially distributed approach of van Huijgevoort et al. (2016) (∼ 2 years), and EcH 2 O-iso conversely simulates much younger water ages in the groundwater both on the hillslope and in the valley bottom (∼ 0.9 and ∼ 2.2 years vs. ∼ 2.9 and ∼ 3.4 years, respectively) and on the hillslope soils (∼ 0.2 years vs. ∼ 0.8 years) than van Huijgevoort et al. (2016). Keeping in mind that these discrepancies might arise from differences in modelling and calibration approaches, these mismatches may also confirm a tendency of EcH 2 Oiso to overemphasise the role of the riparian area as a hydrologic buffer and mixing zone, as well as the contribution of groundwater, in damping the stream isotope response which could be addressed by strategies suggested in preceding sections.
The EcH 2 O-iso model presented in this study is, to our knowledge, the first to simulate catchment dynamics of water isotopes ( 2 H and 18 O) and age by combining a physically based description of hydrological stores and fluxes, a spatially distributed simulation domain, a predictive vegetation component, and non-conservative isotopic processes (evaporative fractionation). Evaluated against a multi-site, extensive isotopic dataset encompassing a wide range of ecohydrological compartments (soil moisture, groundwater, plant xylem, and stream water) across hydropedological units, the model has generally shown good performance in reproducing the seasonal and higher-frequency variations of absolute and relative isotopic content (δ 2 H and lc-excess, respectively). Despite some limitations, this isotope-based evaluation suggests a reasonable capture of the velocity fields (i.e. how fast water parcels move) across the catchment, and complements a previous calibration and evaluation mostly using hydrometric observations (water fluxes and storage dynamics) which indicated a good simulation of catchment functioning from a celerity viewpoint (i.e. how fast energy propagates via the hydraulic gradient) (Kuppel et al., 2018a). Satisfying this dual velocity-celerity perspective is key to characterising water pathways and quantifying the associated travel times in different ecohydrological compartments of headwater landscapes. Complementing more conceptual approaches, the physical basis of the EcH 2 O-iso model further provides the potential to extrapolate these insights beyond recorded conditions and scales, and to notably project the reciprocal feedbacks between plant water use, hydrological pathways, and potential environmental changes. The relatively simple conceptualisation of compartment-scale velocities, e.g. assuming complete mixing and without site-specific parameterisation, and the absence of isotopic calibration, already make the current results particularly encouraging. It also provides a useful framework for hierarchising model development and benchmarking needs. For example, some of the model-data discrepancies in our results stress the necessary incorporation of partial mixing hypotheses, likely to be critical in drier and/or flatter landscapes where diffusive water movement prevails. Second, our model-data analysis of isotope dynamics strongly reflects fractionation effects, be it via soil evaporation or species-specific plant water use. Finally, the versatility of climatic settings in which the original EcH 2 O model has already been evaluated facilitates applying the presented methodology beyond the specifies of a high-latitude, low-energy, wet and steep headwater catchment such as the one simulated here. Further, the flexible spatial domain used by the model will help to provide a process-based modelling framework for plot-to catchment-scale hypothesis testing. This is timely for current challenges in critical zone science, such as exploring the occurrence and mechanisms behind the postulated ecohydrological separation of water fluxes (Berry et al., 2017).
Code and data availability. The source code of the EcH 2 O-iso model used in this study corresponds to the branch master_1.0 in Kuppel (2017). The presented datasets are available from the PURE data repository of the University of Aberdeen (Kuppel et al., 2018b).