Storage dynamics in hydropedological units control hillslope connectivity, runoff generation, and the evolution of catchment transit time distributions.

We examined the storage dynamics and isotopic composition of soil water over 12 months in three hydropedological units in order to understand runoff generation in a montane catchment. The units form classic catena sequences from freely draining podzols on steep upper hillslopes through peaty gleys in shallower lower slopes to deeper peats in the riparian zone. The peaty gleys and peats remained saturated throughout the year, while the podzols showed distinct wetting and drying cycles. In this region, most precipitation events are <10 mm in magnitude, and storm runoff is mainly generated from the peats and peaty gleys, with runoff coefficients (RCs) typically <10%. In larger events the podzolic soils become strongly connected to the saturated areas, and RCs can exceed 40%. Isotopic variations in precipitation are significantly damped in the organic-rich soil surface horizons due to mixing with larger volumes of stored water. This damping is accentuated in the deeper soil profile and groundwater. Consequently, the isotopic composition of stream water is also damped, but the dynamics strongly reflect those of the near-surface waters in the riparian peats. "pre-event" water typically accounts for >80% of flow, even in large events, reflecting the displacement of water from the riparian soils that has been stored in the catchment for >2 years. These riparian areas are the key zone where different source waters mix. Our study is novel in showing that they act as "isostats," not only regulating the isotopic composition of stream water, but also integrating the transit time distribution for the catchment. KEY POINTS: Hillslope connectivity is controlled by small storage changes in soil unitsDifferent catchment source waters mix in large riparian wetland storageIsotopes show riparian wetlands set the catchment transit time distribution.

1. Introduction

Quantifying the processes of water and tracer transport through catchments remains a key research frontier in hydrology, where new technologies, theoretical developments, and novel modeling approaches provide fresh insights into the dynamic controls on streamflow generation and solute fluxes [e.g., ; ; ; ]. Much of this interest has been stimulated by tracer-based insights into the so-called “old water paradox” whereby short-term rainfall-runoff dynamics (over minutes and hours) controlled by the celerity of hillslope responses mobilize water that has usually been stored in the catchment for much longer periods (months to years) but constrained by low pore velocities [; ]. Input-output studies of conservative tracers have provided invaluable insights into the emergent properties of such catchment scale responses [e.g., ]; they have also been useful for catchment intercomparison [e.g., ; ] and convenient for hypothesizing dominant flow paths and mixing processes []. Nested tracer studies integrating soil profile-hillslope-catchment scales over prolonged periods and combining empirical observations with quantitative modeling are a particularly efficient route to improved understanding of the spatial and temporal dynamics of how water is partitioned, stored, and discharged [; ; ]. However, such studies are still quite rare.

Exactly how the spatial distribution of hydrological processes in different hydropedological units integrates to control runoff generation and tracer transport at the catchment scale is still poorly understood [; ] and has proved to be highly variable between different geographical regions []. Recent work [e.g., ] has built on earlier studies [e.g., ; ; ] in emphasizing the ecohydrological importance of soils and their spatial distribution in controlling the catchment hydrological response. The emerging field of hydropedology [] underlines the often crucial role of pedology as a first-order control on catchment runoff generation and solute transfer [e.g., ]. Coupled with this is the role of soil wetness in determining surface and near-surface hydrological connectivity between hydropedological units and the channel network [; ]. However, interrelationships between soils and the deeper subsurface are also important []. There is evidence that even in steep montane catchments, sedimentary [e.g., ], igneous rocks [e.g., ; ; ; ], and superficial drift [] can have significant groundwater storage and affect runoff dynamics.

In northern regions, glacial history and cold, wet hydroclimatic conditions often result in landscapes with distinct hydropedological units and associated ecohydrological communities [; ; ]. Peats are often dominant in low-lying areas of riparian zones and can be differentiated from more freely draining soils in upslope areas [; ]. The absence of trees in such environments often means that the connections and relationships between different landscape units are superficially clear. This provides opportunities for detailed multiproxy investigations of hillslope hydrological connectivity that can test hypotheses regarding hydrological processes [e.g., ; ; ]. It is critical that the fundamental hydrological processes in northern catchments are understood in order to protect ecosystem services in an era of unprecedented environmental change []. Rates of climatic warming are accelerating, with many regions already experiencing warming, with more rapid snowmelt [] and wetter winters and drier summers [; ]. This climatic amelioration is expected to increase development pressures in terms of land use change and exploitation of renewable energy sources (e.g., wind power and forest biofuels) [; ]. Disturbance of near-surface hydrologic processes from such developments may compromise the provision of downstream water resources that provide potable supplies, sustain ecosystems, and support hydropower [].

Nested tracer studies integrating the soil profile-hillslope-catchment scales allow the filtering effect of catchments on tracer input signals to be investigated directly. The passage of conservative, environmental tracers through and between different hydropedological units can be monitored to test hypotheses about how input signals are damped and lagged by internal mixing processes and connections between different spatial units, including deeper groundwaters []. Key to this is understanding the interrelationships between storage and discharge in different landscape units, and how these integrate to control the storage-discharge relationships of the catchment []. Environmental tracers play a critical role in such studies, providing insight into the transit time distributions (TTDs) of catchment and a measure of the internal storage in terms of the volume needed to damp tracer inputs [; ].

Building from this context, we present the results of a 1 year intensive study in a catchment in the Scottish Highlands where peat-dominated soils exert a strong influence on runoff generation. Previous work has identified the main runoff generating areas in valley bottom wetlands from mapping hydropedological units and associated tracer studies []. This has been used as a basis for tracer-aided modeling that has helped us to identify the dynamics of the dominant processes [,2011b; ] and the relationships between storage and discharge [] at the catchment scale. Here we examine the internal dynamics of the major hydropedological units to understand the patterns of water storage, connectivity, and flux and how these affect the filtering of tracer input signals. Using hydrometric and isotopic data, the specific objectives are to understand: (1) where water is stored in different hydropedological units and the associated temporal dynamics; (2) how such storage dynamics affect hydrological connectivity between different units and how this governs the integrated streamflow response; and (3) how connectivity and mixing in different parts of the catchment integrate in riparian wetlands to control the transit time distribution of stream water.

2. Study Site

This study was based in the Bruntland Burn (BB; 3.2 km2), a tributary of the Girnock Burn, a long-term monitoring site in the NE of Scotland (Figure 1). The catchment is granite-dominated with associated metamorphic rocks and has been glaciated, with a wide flat valley bottom at around 250 m above sea level (asl), increasing in gradient to steeper slopes up to around 550 m asl (Figure 1a). The landscape divides into three main hydropedological units following classic soil catenas (Figure 1b). On steeper slopes podzols cover ∼55% (1.76 km2) of the area, but these become thinner and grade into regosols on the interfluves covering ∼10% (0.32 km2) of the catchment. As the gradient declines on lower slopes, these grade into peaty gleys which cover ∼25% (0.8 km2) of the catchment. In the flatter riparian zones, peats (up to 4 m deep) cover 10% (0.32 km2) of the area, overlaying drift deposits that geophysical resistivity surveys show are permeable, stratified, and ∼30 m deep. In these valley bottom areas, the water table is close to, or above, the ground surface. On the steeper hillslopes, the drifts are still 5–7 m deep at around 350 m asl and largely absent above 400 m asl (J. Bradford, Boise State University, personal communication). On these steeper slopes, the water table is usually >1 m deep.

Figure 1

Bruntland Burn catchment showing equipment locations and (a) topography and (b) main hydropedological units; soil stations 1–3 (1, riparian peat soil; 2, peaty gley soil in transitional area between steeper hillslope and valley bottom; 3, freely draining podzol on steeper hillslope). Groundwater level loggers were installed at each soil station and the other sites shown.

The soils have distinct hydrological characteristics reflecting their pedogenesis and support different vegetation communities. The podzols have a freely draining profile: a 15–20 cm O/A horizon mantles a 15–20 cm E horizon, with an illuvial Bs horizon beneath. These horizons thin at higher altitudes grading to shallow rankers above 400 m asl. The vegetation is dominated by heather (Calluna sp. and Erica sp.) moorland. The peaty gleys have a ∼20 cm deep O horizon which overlies gleyed Eg and Bg horizons. Molinia sp. dominates the vegetation. Fringing the stream, as the gradient decreases, peat soils form, characterized by a 20 cm deep less humified O1 horizon that acts as the Spagnum-dominated acrotelm overlying more humified deeper horizons forming the catotelm which can be up to 4 m deep.

Mean annual precipitation (P) is approximately 1100 mm, and evapotranspiration (ET) is relatively low (∼400 mm; Figures 2a and 2b). ET (estimated using Penman-Monteith) is highest in summer. Snow occurs but usually comprises <10% of the annual precipitation. P is evenly distributed with limited seasonality, and most P falls in low-intensity frontal events, with 50% falling in events <10 mm and 75% in events of <20 mm. Most events instigate a streamflow response, though this increases nonlinearly in wetter periods as the saturated zone in the valley bottom expands and connects runoff generating areas to the channel network [].

Figure 2

Spatial and temporal variability in water balance components in (a) precipitation; (b) discharge (black line) and evapotranspiration estimated by Penman-Monteith equation (red line); and soil moisture in the (c) freely draining podzol on steeper hillslope; (d) peaty gley soil in transitional area between steeper hillslope and valley bottom; and (e) riparian peat soil.

3. Data and Methods

Monitoring occurred between 1 June 2011 and 31 May 2012. While the annual precipitation was close to average, the distribution was unusual with a wet summer in 2011, a relatively dry winter 2011/2012, with a warm, dry March 2012, followed by a wetter late spring. Despite this distribution the catchment transitioned between periods of wetness and dryness, which are typical of the annual ranges observed for the 45 year monitoring period for the Girnock catchment. Monitoring in the BB commenced in 2007; precipitation and streamflow were measured, and two Campbell Scientific automatic weather stations (AWS) in the Girnock are located within 2 km of the BB catchment []. Precipitation and stream water have also been sampled for isotope analysis. Sampling takes place on at least a weekly basis, though for this study and the hydrological year 2008/2009 sampling was daily. The samples are collected by automatic samplers (ISCO 3700) which accumulates daily precipitation at 00:00 A.M. bulking four individual stream samples collected at 6 h intervals. The samples were preserved by paraffin, collected weekly for analysis, and analyzed with a Los Gatos DLT-100 laser isotope analyzer following a standard measurement protocol (precision of ±0.4‰ for deuterium; ±0.1‰ for oxygen-18). Isotope signatures are reported in the δ notion (‰) after calibration using Vienna Standard Mean Ocean Water standards. Due to greater precision deuterium was mainly used in the analysis that follows.

Extensive mapping of soil and vegetation [cf. ] and spatially distributed surveys of soil moisture and various tracers [e.g., ] informed the identification of a hillslope transect for detailed investigation of water and tracer movement through and between the main hydropedological units (Figure 1). While it is recognized that one hillslope will capture only a portion of the soil heterogeneity at the catchment scale, the main landscape differences are incorporated in the transect. The transect comprises three stations where soil moisture, groundwater levels, and isotopes were monitored. Volumetric soil moisture content (vsmc) was measured with Campbell time-domain reflectometer (TDR) probes inserted at depths of 10, 30, and 50 cm within the profiles; corresponding to major horizons. At each site duplicate probes were inserted ∼2 m apart and logged at 15 min intervals by a common Campbell data logger. The TDR probes were calibrated against volumetric water content of field samples of each horizon analyzed in the laboratory. Soil water samples were collected from 10, 30, and 50 cm depths using Rhizon soil moisture samplers (Rhizosphere Research Products, Wageningen, Holland). Disposable 60 mL syringes were attached to samplers via a Luer-lock to collect soil pore water under vacuum. After 2 h duplicate samples were collected for isotope analysis. Samples were collected weekly, though frosts occasionally precluded sampling due to soil water being frozen or sample water freezing in the tubes.

Two groundwater wells were installed at each station; water levels were monitored at 15 min intervals using Odyssey capacitance probes at 1.5 m depths at the peat and peaty gley sites. The stony nature of the substrate in the C horizon of the podzol precluded installation below 0.5 m. A further seven wells were distributed around the catchment, particularly in the peat and peaty gley soil units to assess the effects of spatial variability on shallow groundwater dynamics (Figure 1a). Because of the limited depth attainable with hand augured wells (typically <1.5 m), two groundwater springs were sampled for isotope analysis (Figure 1). These are perennial springs, and previous work has shown that they drain deeper groundwater (>5 m) stored in the glacial drift deposit on the catchment hillslopes [].

Time series were corrected removing spurious outliers, and where duplicate measurements exist (groundwater levels, soil isotopes, and soil moisture in particular horizons) means were used (for replicate soil water isotope variability averaged <1‰, soil moisture variability averaged <0.03 vol/vol at 10 cm and <0.05 vol/vol at 30 and 50 cm, and groundwater level variability averaged <2 cm). Precipitation and discharge were integrated to mean daily values. From the daily catchment precipitation we calculated a 7 day antecedent precipitation index (API7) using a negative exponential function representing antecedent wetness. Previous work demonstrated that a calibrated 7 day recursive algorithm best represents the expansion and contraction of the saturation area linking discharge dynamics to hydrometeorological conditions []. This was used to relate hydrological signatures such as runoff coefficients (RCs) and event water contributions to antecedent wetness. We calculated daily RCs as discharge-precipitation ratios selected only for rainfall events with a clear runoff response indicated by a rising hydrograph. Furthermore, a simple two-component hydrograph separation estimated daily event and pre-event water contributions to flow: where Q is the pre-event contribution to total discharge Q, δ (‰) is the event signature in precipitation, and δ (‰) the observed isotope signature in streamflow, and the pre-event signature in streamflow δ (‰) was defined as the previous day streamflow signature.

We applied lumped convolution integral models to estimate the likely form of TTDs for each isotope record in the soils, groundwater, and stream water in relation to the input variability of precipitation isotope signatures: where δ (t) is the stream tracer composition in ‰, δ (t−τ) is the input composition (‰) at time t−τ, and g(τ) is the transfer function. We used the two-parameter gamma distribution (α as shape and β as scale) as the transfer function [see for details] in combination with daily and weekly input functions according to sample frequency:

The mean transit time (MTT) is calculated as the product of α and β (days). Volume-weighted precipitation inputs were corrected for evapotranspiration which was estimated using the Penman-Monteith method using data from an AWS nearby.

The TTD models were first optimized using a differential evolution generic algorithm applying wide initial parameter intervals allowing the search algorithm to converge selecting only the best mutation of parameter set generations []. Model likelihood was assessed minimizing the root-mean-square error (RMSE). Optimal solutions from this initial optimization were then used to inform prior parameter distributions for a Markov chain Monte Carlo (MCMC) approach. An upper temporal limit for the scale parameter β was imposed at 5 years due to the inability of isotope data in detecting water ages greater than this []. MCMC is based on a random walk through the parameter space, where each sampled parameter is associated to the total likelihood. The MCMC was chosen to identify parameter combinations resulting in high posterior probabilities of the model and use this to assess associated parameter uncertainty of these simple TTD models reported as the 5th and 95th percentiles around the mean of the posterior parameter distributions. The posterior parameter distributions were obtained after 500,000 to 2 million iterations of the Markov Chain from which 10,000 parameter sets were retained. The convergence of the chain (full exploration of the parameter space) was visually assessed using parameter traces similar to and a negligible autocorrelation of the retained parameter sets. This analysis was conducted using the R packages FME and Coda []. The 12 month time series were looped for 5 years using the mean output isotope signature as a constant input to achieve approximate mass balance prior to the transit time assessment []. While the relatively short run of data precluded a full assessment of the time variance of the TTD, the robust statistical framework quantified the inherent uncertainty. For the daily precipitation-streamflow data, we examined 6 month “moving window” time slices to explicitly gain insight into the potential time variance [cf. ]. We also carried out an identical analysis of the daily data collected from 2008 to 2009 (see above) for interannual comparison; this year had 10% more precipitation which was more evenly distributed [].

4. Results

4.1. Soil Water Dynamics

The three units showed distinct soil moisture dynamics in response to climatic variability (Figure 2 and Table1). The freely draining podzol on the steeper hillslope (Figure 2c) exhibited clear wetting and drying phases in response to events. In the organic-rich surface horizon the vsmc reached 45% in wet periods falling rapidly to ∼30% during recessions after which it changed little until subsequent large events. The mineral subsoil vsmc showed similar dynamics, with the E and B horizons varying between 35% and 40% in wet periods to ∼25% in dry periods. Wetting fronts moved quickly through the profile consistent with active preferential flow via macropores in the root zone, with the 50 cm sensor always responding within 60 min of the 10 cm sensor. The soils were in more dynamic wetting/drying phases for around 30% of the study period, and moisture content was low and declined slowly in the remainder. Extrapolating from the vsmc measurements with idealized soil horizon depths shows that the water storage in the upper 60 cm varied between 168 mm in dry periods and 238 mm when the profile was at its wettest (Table2).

Table 1

Summary Statistics of Hydrometric Time Series for the Study Period From 1/6/2011 to 31/5/2012ab

Data	Unit	n	Mean	Min	Max	St. Dev.
Discharge (Q)	mm d−1	366	1.52	0.31	11.2	1.5
Precipitation (P)	mm d−1	366	2.41	0	32.5	4.6
Evapotranspiration (ET)	mm d−1	366	1.26	0	6.2	1.2
Groundwater depth 1	cm	35,086	−3.8	−9.8	6.3	1.2
Groundwater depth 2	cm	32,932	−18.5	−42.9	−8.2	5.7
Groundwater depth 3c	cm	13,716	−22.9	−35.0	−1.85	9.5
Soil water 1 O (0–20 cm)	%	33,895	0.81	0.78	0.84	0.02
Soil water 2 O (0–20 cm)	%	34,738	0.74	0.71	0.76	0.01
Soil water 2 E (20–40 cm)	%	34,947	0.35	0.33	0.37	0.008
Soil water 2 Bg (40–60 cm)	%	34,947	0.35	0.34	0.35	0.002
Soil water 3 O/A (0–20 cm)	%	34,943	0.33	0.3	0.45	0.024
Soil water 3 A/E (20–40 cm)	%	34,943	0.29	0.25	0.39	0.028
Soil water 3 Bs (40–60 cm)	%	34,943	0.27	0.24	0.35	0.031

n, number of observations; Mean; Min, minimum, Max, maximum, and St. Dev., standard deviation.

Numbers 1–3 refer to sampling stations: 1, riparian peat soil; 2, peaty gley soil in transitional area between steeper hillslope and valley bottom; 3, freely draining podzol on steeper hillslope.

Data series with data gaps as level did not measure when water table dropped below soil profile.

Table 2

Estimated Water Storage in Soil Profiles (Based on Soil Water Ranges in Table1 and Integrated Over Profile (for Station 1) or Horizon Depths)a

Data	Min (mm)	Max (mm)
Station 1 O (0–100 cm)	820	840
Station 2 O (0–20 cm)	142	152
Station 2 E (20–40 cm)	66	74
Station 2 Bg (40–60 cm)	68	70
Station 3 O/A (0–20 cm)	60	90
Station 3 A/E (20–40 cm)	50	78
Station 3 Bs (40–60 cm)	48	70

Numbers 1–3 refer to sampling stations: 1, riparian peat soil; 2, peaty gley soil in transitional area between steeper hillslope and valley bottom; 3, freely draining podzol on steeper hillslope.

The peaty gley in the transition zone between the steeper slope and valley bottom (Figure 2d) showed little variability in vsmc and remained at or close to saturation. The vsmc was highest in the organic surface horizon varying between 71% and 76%. This decreased in the minerogenic, stony subsoil to ∼35% with little variability. This state of saturation in the soil matrix indicates little additional storage capacity, with the inference that precipitation events will rapidly initiate lateral flows in the larger pores in the organic surface horizons. Indeed, a first approximation of the range in profile storage indicated variation between 278 and 296 mm—a maximum storage change of 18 mm—between the driest and wettest periods (Table2).

The riparian peats also exhibited little variability with the soil matrix remaining saturated all year-round. At a depth of 10 cm the vsmc varied between 78% and 84% again (Table1) indicating only 12 mm dynamic storage change in the upper 20 cm of the profile and consistent with rapid runoff response to precipitation events from saturation overland flow. The total matrix storage in the upper 1 m of the peat varied between 828 and 840 mm (Table2).

4.2. Water Table Dynamics and Streamflow Responses

Water table variability in the three units reflects the dynamics of the soil moisture (Figure 3). On the steeper slopes, wettest periods resulted in the podzol profile becoming saturated to within 10 cm of the surface (Figure 3b). However, the water table fell below the soil profile between events. The water table was present in the B horizon for around 30% of the study period, while it extended into the O horizon <5% of the time. This indicates slower vertical drainage in the lower profile (which tracer tests indicate have Ks values of <0.5 m d−1) but rapid transient lateral flow in the O horizon (Ks > 10 m d−1) in larger events.

Figure 3

Spatial and temporal variability in groundwater levels: (a) discharge, (b) steep podzolic slope, (c) transitional lower slope, and (d) riparian peatland. Groundwater levels are measured as depth beneath surface in centimeters.

On the transitional slope, the water table was generally within 20 cm of the soil surface (Figure 3c). The only exception was the drier period in March when it dropped to 43 cm (Table1). The water table response was extremely dynamic, and most precipitation events resulted in a 5–10 cm rise which took it into the organic O horizon initiating lateral flow. Indeed, lateral flow from soil and groundwater seepage upslope as well as direct precipitation contributes to this response. Recovery after March was immediate.

In the riparian peatland, the water table was always within 10 cm of the soil surface, reflecting saturated conditions due to the retentive nature of the peat, the low topographic gradient, and seepage from upslope (Figure 3d). Empirical evidence has shown that during wetter periods the water table is lying up to 10 cm above the ground surface with ubiquitous surface saturation. This implies more-or-less continual seepage from the peat toward the river in the acrotelm and development of increasingly deep water “ponding” on the profile during precipitation events when upslope drainage accumulates as it migrates laterally toward the channel as saturation overland flow. This is evident in wetter periods in Figure 3d in August, October, December, and April, though the microtopography of the riparian peatland masks this effect.

The stream response to precipitation events was consistent with that of soil water and groundwater (Figures 2 and 3). Precipitation events as small as 3 mm initiated responses in flow, with simultaneous water table response in the transitional lower slope consistent with displacement of water from the permanently saturated riparian wetlands. In most small to moderate events (<5 mm d−1) RCs were generally 5–30% depending on antecedent conditions and precipitation intensity (Figure 4). In larger events, when the steeper podzols are fully wetted and the saturation area expanded, the connectivity appears to be more efficient, with RCs exceeding 40%. The highest values occurred when there were snowmelt events, usually associated with rain-on-snow.

Figure 4

Selected runoff coefficients plotted as time series against discharge. The inset box shows the relationship of runoff coefficients with the 7 day antecedent precipitation index.

This was also consistent with the lag times between the groundwater and streamflow responses exhibiting spatial variability when all 13 groundwater wells in the catchment were considered (not all data shown). Depending upon antecedent conditions, the riparian wells close to the channel exhibited peak water table levels 1–5 h before peak discharge, while those on the upper hillslopes peaked 1–3 h after the highest flow (Figure 5). This also tended to result in positive clockwise but narrow hysteresis curves in the relationships between water table levels in riparian wells and stream discharge (Figure 5a). This was usually reversed deeper in the podzol, indicating drainage contributions on the stream recession (Figure 5b). The nonlinear response of the steeper hillslopes was also evident, where for the larger events the water table was close to the ground surface, but in moderate events only reached into the B horizon. Consequently, the magnitude of water table variations tends to increase with distance from the channel. This resulted in a positive relationship between altitude and mean water table depth and the degree of water table fluctuations.

Figure 5

Depth to groundwater versus streamflow for two of the wells for the entire study period. Markers were plotted at 15 min intervals: (a) an example for a riparian well and (b) an example for a well at steeper hillslope. Arrows show general direction of hysteresis loops.

4.3. Isotope Dynamics

As previously reported for the BB [cf. ], the δ2H variability in precipitation (mean = −56.3‰; range −143.0 to −12.6‰) is damped in stream water (mean = −58.1‰; range −65.8 to −53.6‰; Table3). Precipitation signatures were highly variable, and while seasonality was evident, the atypical nature of the study year is apparent with large fluxes of depleted inputs in the cool summers, as well as more modest depleted inputs in the winter (Figure 6a). Despite the damping, stream water did reflect the major departures of isotope inputs. However, use of mixing models with daily precipitation and streamflow data showed that “new” water from each day's incoming precipitation rarely accounted for 20% of runoff peaks and was usually <10% (Figure 7). Some of the highest new water contributions occurred when the API7 was low, though this paradox might be explained by winter events when snowmelt or frozen soils resulted in greater new water contributions. Likewise, the first month after a warm, dry period in April/May 2011 may have resulted in hydrophobicity in the organic horizons of the drier soils and increased new water inputs. During these warmer periods stream water also became enriched consistent with fractionation in the wet soils in the valley bottom (Figure 8).

Table 3

Summary Statistics of Unweighted Isotope Time Series (δ2H) of Different Waters for the Study Period From 1/6/2011 to 31/5/2012a

Data	Unit	n	Mean	Min	Max	St. Dev.
Stream	‰	317	−58.1	−65.8	−53.6	1.8
Rain	‰	192	−56.3	−143	−12.6	23.2
Groundwater springs	‰	42	−61.2	−63.2	−58.3	1.1
Soil water 1 (in saturation area; 10 cm depth)	‰	45	−55.9	−61.6	−50.7	2.9
Soil water 1 (in saturation area; 30 cm depth)	‰	44	−59.2	−61.8	−56.5	1.1
Soil water 2 O (0–20 cm)	‰	46	−57.3	−61.1	−53.6	2.2
Soil water 2 E (20–40 cm)	‰	44	−59.1	−63.6	−55.4	2.0
Soil water 2 Bg (40–60 cm)	‰	46	−60.5	−63.9	−55.6	2.2
Soil water 3 O/A (0–20 cm)	‰	47	−58.8	−82.9	−44.4	9.6
Soil water 3 A/E (20–40 cm)	‰	46	−59.2	−63.8	−53.7	3.1
Soil water 3 Bs (40–60 cm)	‰	43	−59.3	−66.7	−53.9	3.1

n, number of observations; Mean; Min, minimum; Max, maximum; and St. Dev., standard deviation.

Figure 6

Isotopic composition of (a) precipitation (weighted, shown as proportional circles of deuterium flux) and runoff; soil waters in the (b) freely draining podzol on steeper hillslope, (c) peaty gley soil in transitional hillslope area, (d) riparian peat soil, and (e) in two deeper groundwater springs.

Figure 7

(top) Daily deuterium values for precipitation (weighted, shown as proportional circles of deuterium flux) and runoff relative to streamflow. (bottom) Daily proportion of pre-event water and inset relative to antecedent precipitation index.

Figure 8

Isotope (δ2H and δ18O) signatures of rainfall, soil water in the riparian peat (10 and 30 cm depth), groundwater, and stream water. The solid line represents the global meteoric water line (GMWL): δ2H = 8δ18O + 10.

The response of soil waters to isotopic inputs varied between sites (Table3 and Figure 6). In the podzol upper horizons, isotopic variability was greatest and rapidly reflected precipitation inputs and the effects of evaporative fractionation (Figure 6b). Thus, δ2H showed initial enrichment in June 2011, followed by depletion in July and August, then further enrichment in the autumn of 2011 prior to marked depletion in January 2012 and recovery in the spring. At 30 cm, a similar, but more muted dynamic was observed and one that was lagged by around 2–4 weeks. At 50 cm the variability was all but completely damped.

Soil water in the peaty gley showed similar sequential down-profile damping but more restricted initial variability in the upper horizon compared with the podzol (Figure 6c and Table3). This is consistent with the much greater storage of water for the precipitation to mix with, particularly in the organic surface horizon (Table2). In the O horizon at 10 cm, initial summer enrichment was followed by depletion in August–September, with changes lagging precipitation by around 2 weeks. The enrichment in the autumn was followed by moderate depletion in January/February and limited change thereafter. A similar, but damped dynamic occurred at 30 cm, while less marked variation was evident at 50 cm.

In the riparian peat at 10 cm, the isotopic variation was damped but tracked precipitation with depletion occurring in July/August 2011, December 2011/January 2012, and then finally in May 2012 (Figure 6d). However, at 30 cm isotopic concentrations remained relatively constant. This 30 cm response was similar to groundwater in the two springs which exhibited highly damped isotopic dynamics (Figure 6e). Generally, groundwaters were depleted compared to stream and soil waters, and G1 was slightly more depleted than G2 (Figure 6e). This would be consistent with recharge mainly occurring in winter when precipitation is usually highest and most depleted and effects of evapotranspiration are limited.

The isotopic composition of stream water exhibited dynamics that generally lie between the 10 and 30 cm deep waters from the riparian peat soils, with the latter being similar to groundwater (Figure 8). Only in some larger events did stream water values deviate outside these toward that of precipitation. The mixing model for daily stream samples indicated that “old” pre-event water generally accounted for 80–95% of streamflows (Figure 7). This is consistent with the hydrometric and isotopic data indicating the riparian area acts as the main source of streamflow where near-surface and deeper groundwater sources mix. Previous modeling showed that around 25–40% of annual stream discharge is derived from return flows of this hillslope drainage []. With additional inputs from the saturated riparian zone, this provides 52–65% of streamflow over the course of the year, with the remaining 35–40% of accounted for by the direct discharge of groundwater to the channel.

4.4. Transit Time Distributions

The isotope dynamics were characterized quantitatively by fitting TTDs to the time series of stream water, soil water, and groundwater using a gamma function. Given the potential significance of nonstationarity [cf. ] we used daily precipitation data to fit TTDs to the daily stream water time series (Table4 and Figure 9). We used 6 month time slices with a moving window approach [cf. ] through the year to test how the derived TTD changes with changing hydroclimatological conditions. Table4 reports the accepted models, goodness of fit, and associated parameter uncertainty derived from the MCMC approach. For wetter periods the TTDs were all similar with acceptable models having identifiable α mean values in the general range 0.7–0.8 resulting in constrained posterior distributions compared to the priors reported in Table4. However, β was generally not identifiable with large uncertainties as determined by the posterior distributions resulting in large uncertainties in the derived MTT. Measures of the goodness of fit in terms of the coefficient of determination (R2) and the NSE statistic were reasonable (>0.5 and >0.2, respectively) especially considering the effects of fractionation in warmer months for daily samples. The optimized α increased when the drier March period was included, though the model fit was poorer with low NSE and R2 values (Table4). However, the isotopic signal in precipitation at this time was also less variable, so that—as well as the drier conditions per se—contributes to the poorer fit (see also Figure 6).

Table 4

Different Time Windows (1 Year and 6 Months) and Periods (1/6/2011 to 31/5/2012 and 1/10/2008 to 30/9/2009) Applied to Calibrate Daily Gamma TTD Models (Shape Parameter α and Scale Parameter β) Using a MCMC Approach as a Surrogate for Time Variance of TTDsa

Period	Initial [α, β]	Posterior			Best Fit
		α Mean [5th, 95th]	β Mean [5th, 95th]	MTT Mean [5th, 95th]	RMSE (‰)	NSE	R2
1 year: 1/6/11 to 31/5/12	[0–2,1–1832]	0.84 [0.53, 1.24]	1223 [534, 1780]	1027 [283, 2207]	16.2	0.2	0.53
6 months: 1/6/11 to 30/11/11	[0–2,1–1832]	0.74 [0.54, 1.08]	1207 [469, 1778]	893 [253, 1920]	10.3	0.38	0.63
6 months: 1/7/11 to 31/12/11	[0–2,1–1832]	0.74 [0.56, 1.12]	1157 [446, 1763]	856 [250, 1975]	10.7	0.36	0.61
6 months: 1/8/11 to 31/1/12	[0–2,1–1832]	0.71 [0.57, 1.17]	1107 [481, 1780]	786 [274, 2083]	10.2	0.41	0.7
6 months: 1/9/11 to 29/2/12	[0–2,1–1832]	0.83 [0.62, 1.29]	1054 [406, 1663]	875 [252, 2145]	10.8	0.28	0.67
6 months: 1/10/11 to 31/3/12	[0–2,1–1832]	0.84 [0.64, 1.31]	1090 [500, 1773]	916 [320, 2323]	11.2	0.25	0.61
6 months: 1/11/11 to 30/4/12	[0–2,1–1832]	0.86 [0.65, 1.39]	1182 [480, 1764]	1017 [312, 2452]	11.3	0.21	0.54
6 months: 1/12/11 to 31/5/12	[0–2,1–1832]	0.94 [0.67, 1.47]	1179 [510, 1771]	1108 [342, 2603]	10.9	0.13	0.5
1 year: 1/10/08 to 30/9/09	[0–2,1–1832]	0.7 [0.48, 1.22]	1185 [498, 1762]	830 [239, 2150]	19.0	0.41	0.61
6 months: 1/10/08 to 31/3/09	[0–2,1–1832]	0.62 [0.41, 0.96]	998 [273, 1736]	619 [112, 1667]	13.5	0.49	0.73
6 months: 1/11/08 to 30/4/09	[0–2,1–1832]	0.55 [0.39, 0.81]	1171 [418, 1766]	644 [163, 1430]	12.4	0.6	0.8
6 months: 1/12/08 to 31/5/09	[0–2,1–1832]	0.61 [0.4, 1.02]	1229 [509, 1784]	750 [204, 1820]	12.5	0.32	0.57
6 months: 1/1/09 to 30/6/09	[0–2,1–1832]	0.74 [0.49, 1.25]	1177 [456, 1757]	871 [223, 2196]	13.3	0.35	0.58
6 months: 1/2/08 to 31/7/09	[0–2,1–1832]	0.71 [0.5, 1.17]	1204 [495, 1771]	855 [248, 2073]	14.0	0.35	0.56
6 months: 1/3/08 to 31/8/09	[0–2,1–1832]	0.79 [0.52, 1.15]	1206 [503, 1772]	953 [262, 2038]	13.0	0.34	0.56
6 months: 1/4/08 to 30/9/09	[0–2,1–1832]	0.74 [0.52, 1.23]	1187 [491, 1778]	878 [255, 2187]	13.1	0.31	0.52

The models were optimized using a root-mean-square error (RMSE) criterion, but the best fit Nash-Sutcliffe efficiency (NSE) and coefficient of determination (R2) are also given for comparison purposes. Uncertainties are expressed as 5th and 95th percentiles around the mean of the posterior parameter distributions.

Figure 9

Simulated daily stream isotope signatures using the best fit and most extreme (i.e., wettest and driest) MCMC TTDs from monthly time window calibration as a surrogate of time-variable TTDs for the study year 2011/2012 (see Table4 for parameter values and associated uncertainty). The inset box shows the TTDs as cumulative distribution functions over time.

The mean of accepted models for stream water for the 12 month period was 2.8 years, but the 90th percentile range was high (Table4). If the TTD for the wettest 6 month periods (1 August 2011 to 31 January 2012) was used, this reduced the MTT to 2.2 years for the mean of acceptable models. In contrast, parameters fitted to the driest 6 months gave mean acceptable MTT of 3 years, though even the best fit model was poor and such dry conditions over a full year are unrealistic in the context of the Scottish climate (Figure 9). Interyear comparisons using daily data for the wetter year of 2008/2009 showed a fairly consistent TTD, with a mean MTT of acceptable models of 2.3 years varying between 1.7 years for the wettest 6 months and increasing to 2.6 years during the driest (Table4). Figure 9 shows the cumulative distribution function of the best fit and the wettest and driest periods all indicating a substantial proportion (>20%) of tracer recovery still needed after 5 years.

Given the similarities of the daily stream water TTDs averaged over longer-term periods, the limitations of 1 year of data, and the weekly sampling of the soil water and groundwater isotopes, time-invariant TTDs were also fitted to the different source data using the same gamma approach (Table5). Using aggregated weekly input data, model fits for stream water were similar to the daily data, with a MTT of acceptable models around 2.5 years. For best fit models an α of 0.61 was derived with the R2 and NSE statistics improved due to the reduced fractionation.

Table 5

Best Fit RMSE Performance Criterion Applied to Weekly Rainfall Input (P) and Weekly Discharge (Q), Groundwater (GW), and Weekly Soil Water (S1, S2, and S3 at Depth 10, 30, and 50 cm) Output Deuterium Time Seriesa

δ2H	Initial [α,β]	Posterior			Best Fit			MTT (days)
		α Mean [5th, 95th]	β Mean [5th, 95th]	MTT Mean [5th, 95th]	RMSE (‰)	NSE	R2
P-Q	[0–2; 1–1832]	0.81 [0.47, 1.32]	1099 [357, 1757]	889 [168, 2317]	6.0	0.42	0.72	1183
P-GW	[0–20; 1–1832]	5.5 [1.5, 9.6]	945 [140, 1772]	5194 [210, 16,534]	3.2	−0.1	0.11	1043
P-S1_10	[0–2; 1–1832]	0.84 [0.42, 1.33]	1050 [322, 1778]	882 [135, 2366]	6.4	0.48	0.73	1050
P-S1_30	[0–20; 1–1832]	11 [2.6, 19.2]	623 [70, 987]	6853 [182, 18,949]	3.9	0.05	0.09	1806
P-S2_10	[0–10; 1–1832]	2.2 [1.1, 4.8]	728 [135, 1505]	1603 [147, 7224]	10.8	0.51	0.82	518
P-S2_30	[0–10; 1–1832]	2.0 [1.3, 3.7]	686 [110, 1708]	1372 [143, 6321]	12.4	0.29	0.6	609
P-S2_50	[0–10; 1–1832]	3.0 [1.9, 4.7]	469 [217, 1498]	1407 [413, 7042]	7.5	0.53	0.85	931
P-S3_10	[0–10; 1–1832	1.43 [0.41, 2.54]	329 [39, 1414]	469 [16, 3591]	24.9	0.44	0.77	92
P-S3_30	[0–10; 1–1832]	2.17 [0.56, 4.6]	308 [73, 1701]	665 [41, 7826]	16.0	0.08	0.12	154
P-S3_50	[0–10; 1–1832]	1.81 [1.17, 2.23]	581 [56, 1673]	1050 [66, 3731]	18.0	0.06	0.21	205

The initial parameter ranges as well as the posterior parameter distribution characteristics (5th and 95th percentiles around the mean) are given. Additionally to the RMSE criterion the best fit NSE and coefficient of determination (R2) are given. For comparison purposes units are given in days.

Best model fits in the upper soil horizons (in terms of the NSE statistic) were reasonable (>0.4) but generally declined with depth. On the steeper hillslopes this gave best fit MTTs of approximately 3–4 months in the upper soil profile (O horizon) of the podzols, increasing to around 5–7 months in the lower profile, though the lack of variability in the soil isotope response with depth resulted in a poorer model fit (NSE < 0.1). The best fit α parameter and mean of acceptable models were much higher (>1) than stream water (Table5). On the transitional lower slope, the TTDs were better identified (NSE ∼0.5 and R2 >0.6). Here best fit MTTs increased to ∼1.5 years in the upper horizons with down-profile aging to >3 years in the lower horizons. In the valley bottom peatland, the surface horizon of the peat had a TTD which corresponds to a classic gamma function model with an α parameter of 0.58 and a MTT of 2.5 years, which was strikingly similar to the stream water. The TTD for 30 cm in the peat was poorly identified, given the lack of isotopic variability, but inferred MTTs increased to >5 years. The two groundwater springs were also extremely damped, and fitted models were poorly identified and very similar to subsoil horizons with long MTTs.

The spatial variation of the transit times along the hillslope is summarized conceptually in Figure 10 as probability density functions of the TTDs for the soil, ground, and stream waters using the best fit models reported in Table5. This shows the greatest isotopic damping of incoming precipitation occurred as a result of mixing processes in the upper half meter of the catchment soils on the lower hillslopes. The more rapid movement of water through the drier, more freely draining podzol profile was apparent, as was the greater attenuation in the peaty gley. The effects of both down-profile and downslope advection and dispersion are evident for these two soils in Figure 10. However, the TTD in the upper layers of the riparian peatland was very similar to that of stream water. This strongly suggests that drainage of waters from different sources mixes in the large storage zone in the upper 30 cm of the valley bottom peats, and that this is the critical mixing zone in the catchment. This riparian zone integrates the drainage of shallow soil water and deeper groundwater emerging from upslope, mixing it with resident storage and incident precipitation which largely regulates the isotopic composition of the stream and the catchment TTD.

Figure 10

Conceptual hillslope diagram showing the three monitored soil profiles. Transit time distributions were estimated from the best fit weekly input (rain)-output (stream, groundwater gw and soil profiles at different depths) gamma function convolution integrals using a MCMC approach (see Table5 for parameter values and associated uncertainty).

5. Discussion

5.1. Where Is Water Stored and What Are the Associated Storage Dynamics?

The importance of soils—especially, retentive organic horizons—in filtering input isotope signals in peat-influenced catchments is apparent from this work (Figure 6). The storage of water in the soils (∼300 mm at the catchment scale) is large relative to most precipitation events (50% of which are <10 mm), particularly in the valley bottom where soils are saturated all year-round, and there is rapid mixing of inputs (Table1). In these soils the input tracer signals are substantially damped in the upper horizons. The retentive nature of organic soils is well known, and poor drainage afforded by the low topographic and hydraulic gradients ensures that storage is high all year-round. However, even on steeper slopes, inputs are damped in the upper profile (Figure 6) where evaporation also results in fractionation.

In addition to soil water, storage in drift deposits is substantial. Geophysical surveys showed that the drift depth increases to 25–30 m in the valley bottoms where water tables are close to the soil surface (Bradford, personal communication). On the lower slopes, drifts are ∼5–7 m deep up to an altitude of around 350 m asl and are generally absent >400 m asl. These are variable, mainly comprising moraines and ice marginal deposits. The local dominance of granite means that most drift has a relatively coarse sand matrix with abundant clasts (pebbles, cobbles, and boulders). This is permeable and covers ∼70% of the catchment. Assuming an average depth of 10 m with a typical porosity of 15% for drift deposits, this gives a first approximation of storage ∼1500 mm, equivalent to over a year's precipitation. In addition to the soil storage (∼300 mm) this is the same order of magnitude of storage as that derived from tracer damping (approximately 1100–1600 mm on the basis of 2 to 3 year MTTs) or tracer-based hydrological models [, 2011b; ]. Deeper stores of bedrock groundwater and their flow contributions are unknown [cf. ] but maybe important in fractured granite [].

5.2. How Is Connectivity Between Hydropedological Units Affected by Storage Dynamics?

Despite these large residual stores of water at the catchment scale, the dynamic storage changes that regulate the catchment runoff response are small (typically <10 mm). While the lower slopes receive continuous seepage from the steeper surrounding hillslopes as some groundwater reemerges (Figure 1), this connectivity usually produces relatively slow fluxes. Modeling studies indicate that this exfiltration of shallow soil water and deeper groundwaters accounts for 6–13% and 19–27% of annual streamflows, respectively []. Most streamflow responses to smaller precipitation events result in low runoff coefficients (<30%); both the hydrometric and tracer data imply that the stream response is generated from the valley bottom areas. Much greater fluxes from the upland hillslopes occur when connectivity is greatest due to higher hillslope water tables during wet periods when higher runoff coefficients are evident in larger events []. It seems that the delivery of water at such times is more threshold-like [] with higher water tables initiating transmissivity feedback in the upper podzol profiles []. However, the storage change between the wettest and driest states of the podzols is only 70 mm, and the precipitation inputs that trigger strong connectivity are smaller still.

The much larger residual catchment storage—relative to these small dynamic storage changes which drive connectivity—is also reflected in the isotope dynamics of rainfall-runoff transformations. These show that typically <10% of runoff during small to moderate events is “new” event water consistent with displacement of resident soil water in the riparian wetlands (Figure 8). In larger events this increases to ∼20%, indicating steeper hillslopes are connected but contributions are restricted to displacing riparian storage. Exceptions may be where snowmelt, frozen soils, or hydrophobicity results in high new water contributions.

5.3. How Do Connectivity and Mixing in Different Parts of the Catchment Integrate and Control the TTD of the Stream?

The valley bottom riparian area acts as a critical mixing zone where different catchment source waters mix and generate up to 65% of annual discharge []. While the dynamics of streamflow and stream isotope composition are driven by precipitation inputs, these are filtered through different hydropedological units where—as Figure 10 shows—the modifications to the input isotope signal through advection, dispersion, and mixing are markedly different []. The most distinct and variable soil water stores are those in upper horizons of the podzols which are hydrologically disconnected from the channel network most of time. But even here, damping occurs at depth within the profile (Figure 6). Despite this, during both event (through transmissivity feedback) and drier periods (as return flow of groundwater) water from the steeper slopes drains into the riparian wetlands and can mix with a large store of near-surface water. Incident precipitation also mixes with these resident waters as displacement occurs during storm events, and the saturated area may extend to 40% of the catchment []. While such runoff mechanisms have been observed previously in northern catchments, the effects on the catchment TTD have never been demonstrated. Thus, the isotopic variability of the stream is similar to that of the near-surface horizons of the wetlands and can generally be explained by the mixing of shallow soil water with deeper groundwater which also discharges directly into the stream (Figure 8). This results in the stream having a very similar TTD and MTT to near-surface riparian water as this area integrates different water sources contributing to streamflow.

This role of the riparian wetlands as an “isostat” which largely determines the stream water TTD is consistent with other studies in northern catchments that have stressed the significance of the riparian zone in controlling runoff generation and stream water chemistry [; ; ; ; ; ]. This isostat behavior, together with the wet, cool climate [see ], may result in the less temporal variance in the TTDs of the catchment than has been reported from catchments with drier and/or more seasonal climates [e.g., ; ]. While the daily and weekly isotope data used here will underestimate the significance of event peaks [cf. ], both the low new water contributions and the relatively stable longer-term catchment TTDs indicate that these will be moderated by the riparian zones. However, the MCMC analysis indicated considerable uncertainty in transit times, so absolute values should be treated cautiously.

This role of the riparian zone as a “hot spot” for mixing at the catchment scale has important relevance for hydrological models that attempt to integrate tracers [e.g., ]. The mixing assumptions for different model stores/reservoirs are critical elements of tracer-aided models, and assumptions of complete mixing are often inadequate [] resulting in parameterization for partial and dynamic mixing [e.g., ; ]. The results of this study show that the nature of mixing processes varies spatially and the effects on runoff composition depend on the topology of mixing zones in relation to the channel network. Thus, while mixing processes may be partial and nonlinear in some parts of the catchment, others—like the riparian peatlands in this study—may approach more complete mixing as indicated by the TTDs.

5.4. Wider Implications

This study provided insight into the spatial distribution of storage within a northern headwater, the way in which stores fill, connect, and control stream discharge. The hydrometric and tracer data can be used in hypothesis testing by applying rainfall-runoff models to quantify critical internal state variables that characterize the systems function [; ; ] or to understand mixing processes []. However, integrating empirical and modeling work has revealed many unknowns about the functioning of the Bruntland Burn. A key uncertainty is the nature of deeper storage, its dynamics, and contribution to streamflow generation. It is unclear how much of this deeper groundwater (i.e., >10 m deep) interacts with the soil water to influence the streamflow response. Similarly, the significance of contrasting ecohydrological relationships in the different hydropedological units in terms of evaporation losses has been considered only very crudely here [e.g., ]. Given the importance of near-surface processes in the hydropedological units which dominate northern catchments the sensitivity to environmental change is great and the need for protection of riparian zones in particular is clear [; Tetzlaff et al., 2013b].

6. Conclusion

We examined runoff processes in a peat-influenced catchment in the Scottish Highlands using integrated hydrometric and isotope-based techniques. The hydrometric data showed that small changes in dynamic storage in riparian wetlands generated runoff following most precipitation events which were small. In contrast, larger events connected the steeper catchment hillslopes leading to the expansion of saturated areas and a nonlinear increase in runoff response. The isotope data showed that variability in precipitation inputs was almost completely damped in the upper 0.5 m of the catchment soils, with old water accounting for >90% of streamflows in most events. The isotopes also indicate residual catchment storage volumes (>1000 mm) that are two orders of magnitude greater than the dynamic storage changes (<20 mm) generating most storm events. Investigation of a hillslope transect showed downslope aging of waters, and it was clearly shown that the riparian wetlands acted as “isostats,” where the major mixing of different catchment source waters occurred and the isotopic composition and TTD of stream waters were set.