How catchment characteristics influence hydrological pathways and travel times in a boreal landscape

Understanding travel times and hydrological pathways of rain and snowmelt water transported through the landscape to recipient surface waters is critical in many hydrological and biogeochemical investigations. In this study, a particle-tracking model approach in Mike SHE was used to investigate the pathway and its associated travel time of water in 14 partly nested, long-term monitored boreal sub-catchments of the Krycklan catchment (0.12–68 km²). This region is characterized by long and snow-rich winters with little groundwater recharge and highly dynamic runoff during spring snowmelt. The geometric mean of the annual travel time distribution (MTT_geo) for the studied sub-catchments varied from 0.8 to 2.7 years. The variations were related to the different landscape types and their varying hydrological responses during different seasons. Winter MTT_geo ranged from 1.2 to 7.7 years, while spring MTT_geo varied from 0.5 to 1.9 years. The modelled variation in annual and seasonal MTT_geo and the fraction of young water (<3 months) was supported by extensive observations of both δ¹⁸O and base cation concentrations in the different streams. The travel time of water to streams was positively correlated with the area coverage of low-conductive silty sediments (r=0.90, P<0.0001). Catchments with mixed soil–landscape settings typically displayed larger variability in seasonal MTT_geo, as contrasting hydrological responses between different soil types (e.g. peat in mires, till and silty sediments) are integrated. The areal coverage of mires was especially important for the young water contribution in spring (r=0.96, P<0.0001). The main factor for this was attributed to extensive soil frost in mires, causing considerable overland flow during the snowmelt period. However, this lower groundwater recharge during snowmelt caused mire-dominated catchments to have longer stream runoff MTT_geo than comparable forest catchments in winter. Boreal landscapes are sensitive to climate change, and our results suggest that changes in seasonality are likely to cause contrasting responses in different catchments depending on the dominating landscape type.

1 Introduction

The pathways and associated travel times of water through the terrestrial landscape to stream networks is a widely discussed topic in contemporary hydrology. This interest has emerged because of the significant role travel time and routing of water through various subsurface environments play in hydrological and biogeochemical processes (McDonnell et al., 2010; Sprenger et al., 2018). This includes fundamental implications for weathering rates (Burns et al., 2003), transport and dispersal of contaminants (Bosson et al., 2013; Kralik, 2015), and accumulation and mobilization of organic carbon and associated solutes (Tiwari et al., 2017). The travel time, from precipitation input to the outflow into streams, provides valuable information about catchment sensitivity to changes in land use and climate and the fate of long-range transport of contaminants and nutrients deposited with precipitation (van der Velde et al., 2012). The travel time distribution can vary substantially in time and space, depending on catchment characteristics and hydrological conditions, including, for example, slope, catchment size, soil heterogeneity, and seasonality (Botter et al., 2010; Lin, 2010; Heidbüchel et al., 2012; Hrachowitz et al., 2013). Therefore, estimating travel times for contrasting landscape elements is challenging, but when successful, it will enhance our ability to understand and predict catchment functioning more adequately.

Stream water consists of a blend of overland flow and groundwater of different ages. The mean travel time (MTT) to streams is calculated as the average age of this mix (McGuire and McDonnell, 2006). The baseflow is the part of stream groundwater contribution that generally has travelled the furthest and is the oldest (Klaus et al., 2013; Hrachowitz et al., 2016). In contrast, young stream water is typically connected to overland flow or shallow subsurface pathways, which mainly can be seen at times with large rain or snowmelt inputs (Peters et al., 2014; Hrachowitz et al., 2016). The variability of water sources makes the travel time distribution difficult to quantify, especially on intra-annual timescales, as they vary in time and space depending on numerous scale-dependent and scale-independent processes (Botter et al., 2010). A better understanding of the seasonal variability in the fraction of young and old waters can help provide insights into the fundamental role catchment characteristics play in regulating the hydrology and biogeochemistry of streams and rivers.

Stable water isotopes and biogeochemical tracers are common tools applied in field investigations to locate water sources and follow their pathways through the landscape (Maulé and Stein, 1990; Rodhe et al., 1996; Goller et al., 2005; Tetzlaff and Soulsby, 2008). Isotopic tracer signal dampening can provide an estimate of MTT (Uhlenbrook et al., 2002; McGuire et al., 2005; Peralta-Tapia et al., 2016), and more elaborate time-series analysis can offer quantitative assessments of travel times (Harman, 2015; Danesh-Yazdi et al., 2016). However, the isotope amplitude signal used to estimate MTT in many transfer functions is lost after approximately 4 to 5 years because of effective mixing (Kirchner, 2016), limiting the use of isotopes in catchments with long travel times. The young water fraction, often defined as water younger than 2 to 3 months, can, however, still be quantifiable in such catchments (von Freyberg et al., 2018; Lutz et al., 2018; Stockinger et al., 2019). The main advantage of water isotopes is that they are relatively conservative and fractionate primarily because of evaporation. Hence, once in the subsurface environment, the signal is only affected by mixing different water sources. In contrast, many biogeochemical tracers react and transform on their route to streams (Lidman et al., 2017; Ledesma et al., 2018). Such transformation and reactions depend on the specific solute and soil environment that water encounters and, therefore, give qualitative information about groundwater flow pathways (Wolock et al., 1997; Frisbee et al., 2011; Zimmer et al., 2012). Combined information from conservative and reactive tracers can hence provide an enhanced understanding of hydrological processes as their concentrations and dynamics can tell complementary stories about the specific pathways water takes from the source to the recipient stream (Laudon et al., 2011).

A complementary approach to field experiments is numerical modelling, which can help achieve a more complete system understanding. Lumped hydrological models often describe catchments as single integrated entities. In contrast, distributed numerical models can include spatial heterogeneity in input parameters and therefore have the potential to represent catchment processes more mechanistically. In turn, this can lead to a more process-based understanding of hydrology and biogeochemistry at the catchment scale (Brirhet and Benaabidate, 2016; Soltani, 2017). Common methods to calculate travel times using numerical methods include models using solute transport routines and particle tracking (Hrachowitz et al., 2013; Ameli et al., 2016; Kaandorp et al., 2018; Remondi et al., 2018; Yang et al., 2018; Heidbüchel et al., 2020). Models, however, need – as far as possible – proper tests against empirical observations to build confidence in their output. Stream discharge, groundwater levels, and tracer data are examples of such validation data that can provide vital information (McGuire et al., 2007; Hrachowitz et al., 2015; Wang et al., 2017). The collection of such field data is, however, costly and time-consuming. Therefore, data for calibration and validation are often limited, and the minimum length and types in data-sparse catchments are currently a topic of increasing interest (Bjerklie et al., 2003; Jian et al., 2017; Li et al., 2018).

Snow-dominated landscapes have received increasing attention in the last decades due to their importance as water resources (Barnett et al., 2005) and their vulnerability to climate change (Tremblay et al., 2011; Aubin et al., 2018). Landscapes with long-lasting snow cover that often melts rapidly in the spring create both opportunities and challenges for determining the pathways and travel time of water discharging to streams. The long, snow-rich winters not only cause protracted periods of winter baseflow with little or no recharge (Spence et al., 2011; Spence and Phillips, 2015; Lyon et al., 2018), but they also cause considerable amounts of water during the often short and intensive snowmelt in the spring. Although attempts to assess travel times generally have provided useful results using, for example, models to reconstruct isotope signal dampening in snow-dominated catchments, the winter season has proven to be especially challenging, suggesting that other methods to assess travel times may be required (Heidbüchel et al., 2012; Peralta-Tapia et al., 2016). The boreal region also consists of numerous patches of lakes and mires, interspersed in a landscape dominated by coniferous forests on different soil types, which makes this task even more challenging. Hence, accounting for the unique circumstances of both baseflow with long travel times and those of the intensive spring snowmelt with potential large overland flow components in heterogeneous landscapes requires models that can handle the complexity and separation of various flow components across scales, soil types, and landscape patches.

To overcome previous limitations, this study used particle tracking in the physically based distributed numerical model, Mike SHE (Graham and Butts, 2005), to enhance our understanding of stream water contribution in boreal landscapes across seasons and landscape configurations. The water movement model in Mike SHE calculates saturated (3D) groundwater flow and unsaturated (1D) flow and is fully integrated with the surface water and evapotranspiration. The water flow model setup and results previously presented by Jutebring Sterte et al. (2018) were used as the study platform for this work. The model has been calibrated and validated to 14 sub-catchments using daily stream-discharge observations and periodical measured groundwater levels in 15 wells throughout the Krycklan catchment in the boreal region of northern Sweden (Laudon et al., 2013; Jutebring Sterte et al., 2018). The model complexity allows for an in-depth investigation of advective travel times by non-reactive particle-tracking simulations in a transient flow field.

The main objective of this study was to quantify annual and seasonal (winter, spring, and summer) travel time distributions and calculate MTT of water runoff to streams of the Krycklan sub-catchments to disentangle how these are related to physical landscape characteristics and variation in groundwater recharge. Firstly, the credibility of the model results was tested by comparing calculated travel times for the 14 sub-catchments to 10-year observational records from Krycklan, including average seasonal changes in stream isotope signatures and base cation concentrations. The usefulness of stream isotopic composition and chemistry record has previously been demonstrated for understanding the connection of hydrological flow pathways and travel times for this site (Laudon et al., 2007; Peralta-Tapia et al., 2015) but with the limitation of studies on only short periods or single catchments. Secondly, the purpose was to go beyond what was previously done by identifying the connection between travel times and different catchment characteristics and test how this varies depending on the hydrological conditions. This was accomplished by capturing contrasting seasons such as the low-flow conditions in winter with limited input of new precipitation, high flow in spring when the system is still partly frozen, and summer when evapotranspiration (ET) becomes a significant process. We focused especially on the catchment characteristics that have been suggested to be important factors for regulating stream chemistry of the Krycklan sub-catchments, including the areal coverage of mires, catchment size, soil properties, and seasonal changes in groundwater recharge (Karlsen et al., 2016; Klaminder et al., 2011; Laudon et al., 2007; Peralta-Tapia et al., 2015; Tiwari et al., 2017).

2 Method

2.1 Site description

The Krycklan study catchment, located in the boreal region at the transition of the temperate/subarctic climate zone of northern Sweden, spans elevations from 114 to 405 m a.s.l. (Fig. 1, Table 1). The characteristic vegetation of this boreal landscape is the dominance of Scots pine (Pinus sylvestris) and Norway spruce (Picea abies), covering most of the catchment (Laudon et al., 2013). In this study, we refer to soil as all unconsolidated material above the bedrock.

Figure 1 The Krycklan catchment. (a) Locations of sub-catchments and their outlets. The areas are colour-coded based on their stream network connections; e.g. all sub-catchments of one colour connect before reaching the white area. For further details of the catchment characteristics, see Table 1. (b) The figure shows the soil map used in the Mike SHE flow model, which is based on data from the Swedish Geological Survey soil map (1:100 000) and field investigations. (c) Soil depth to bedrock map taken from the Swedish Geological Survey (2016) and is shown in metres below the ground surface (m b.g.s.). (d) Catchment topography, shown as metres above sea level (m a.s.l.).

Table 1 Sub-catchment characteristics. The list includes all 14 monitored sub-catchments in Krycklan, called C1 to C20, including the entire Krycklan catchment, C16. Different branches of the stream network are gathered in the table and illustrated in distinct colours in Fig. 1. The table includes the sub-catchment area, average elevation, and average slope. Further descriptions of these characteristics can be found in Karlsen et al. (2016). The table also includes soil proportion based on the soil map (1:100 000) from the Swedish Geological Survey (2016).

Krycklan has a landscape distinctively formed by the last ice age (Ivarsson and Johnsson, 1988; Lidman et al., 2016). At the higher elevations to the north-west, located above the highest postglacial coastline, the soils can reach up to 15–20 m in thickness. Here, the soil primarily consists of glacial till, and the landscape is intertwined with lakes and peatlands. The deeper soils consist of basal till which was deposited and compacted under the moving ice. In contrast, the shallower till layers consist primarily of ablation till, which is less compact since it mainly has been compacted by its own weight (Goldthwait, 1971). This causes a decreasing hydraulic conductivity with depth, which is characteristic for glacial till in northern Sweden (Bishop et al., 2011; Nyberg, 1995; Seibert et al., 2009). At lower elevations, the soils consist of fluvial and glaciofluvial deposits of primarily sandy and silty sediments. Compared to the soil at higher elevations in the catchment, these deposits can reach thicknesses up to approximately 40 to 50 m and have a hydrological conductivity that is more constant with depth because these soil types have mainly been compacted only by their own weight.

For more than 30 years, multi-disciplinary biogeochemical and hydrological studies have been conducted in Krycklan (e.g. Laudon and Sponseller, 2018). Streamflow is monitored in 14 nested sub-catchments, called C1 to C20, with the longest continuously monitored time series starting at the beginning of the 1980s. Connected by a network of streams, the different sub-catchments allow an evaluation of the effects of catchment characteristics on hydrologic transport, including soil type, vegetation, and differences in topography (Table 1).

2.2 Linking seasonal base cation concentration and isotopic signature to travel times to stream water

This study was focused on three seasons in Krycklan: winter, spring, and summer (Tables 2, A1, Appendix). For evaluation of stream chemistry, we defined the winter as late early December to late February, from the time air temperatures were below 0 ^∘C, until the air temperature started to rise above freezing temperatures again, causing snowmelt. This season is characterized by an extensive and permanent snow cover with little or no groundwater recharge. We assumed that the winter stream composition reflects the chemistry of deeper groundwater (Fig. 2). Similarly, we defined spring as the hydrological period directly influenced by the snowmelt. The main part of the snowmelt and spring flood occurs in April–May. During snowmelt, ca. 50 % of the annual precipitation leaves the system in a short period of time, diluting baseflow with new input of water. Finally, we defined the summer season as the period between July and September when the hydrology is characterized by rain, high ET, and relatively little runoff. March, June, October, and November were excluded because, hydrologically, they are typically transition months between the three distinct seasons. This is because snowmelt influences runoff in March and June. October and November are transitional months between summer and winter conditions, with irregularly occurring snowfall and soil frost events.

Table 2 Seasonal stream chemistry. The table includes average winter signatures (‰) and the average difference between annual winter–spring and winter–summer signatures (Δδ¹⁸O). The table also includes average winter, spring, and summer BC concentrations.

^a δ¹⁸O signature (2008–2018); data have been adjusted according to the lake proportion according to Eq. (A7), Appendix.
^b Base cation concentration (2008–2016); data have been adjusted according to the mire proportion.
^c SD: standard deviation; SEM: standard error of the mean.
^d Measured precipitation average for isotopes (2007–2016). The precipitation average is close to equal to isotope measurements of groundwater below ca. 10 m.

Figure 2 Conceptual figure of travel time to stream vs. stream isotopic signature (a, b) and stream base cation concentration (c). (a) The connection between δ¹⁸O and travel time to stream, where the sine curve shows the annual variations of δ¹⁸O in precipitation, and approximate seasonal winter, spring, and summer stream compositions are marked and exemplified by the average annual changes in C4. In winter, the travel times are related to the average deviation in the isotopic signature between the winter baseflow and the long-term precipitation. In spring, the fraction of young water is correlated with the difference between the average spring stream signature and the average winter baseflow. In summer, the fraction of young water is correlated with the difference between the average summer stream signature and average winter baseflow. (b) Seasonal δ¹⁸O averages for three example streams: C2, C4 and C16. (c) The connection between base cation (BC) concentration and soil contact time. The longer the water spent in the mineral soil, the higher the stream concentrations of BCs will be due to soil weathering.

In this study, stable water isotopes (δ¹⁸O) were used to track pathways of precipitation inputs to stream networks (see the Appendix for the δ¹⁸O definition). Ten years of δ¹⁸O measurements for 13 of the 14 sub-catchments were used. Isotopic fractionation caused by lake surface evaporation affects the isotopic signal of some of the sub-catchments (Leach and Laudon, 2019). This fractionation was corrected by accounting for the percentage of lakes in each sub-catchment (Table 2), using the same principle as Peralta-Tapia et al. (2015) but adjusted to newly acquired δ¹⁸O observations (Eq. A7, Appendix).

The comparison of the modelling results to observations of δ¹⁸O was based on a conceptual model of the seasonal variability and differences between precipitation and runoff (Fig. 2a). The precipitation signal varies on a seasonal basis, creating an amplitude difference (Fig 2b). This amplitude is reduced due to groundwater mixing until complete mixing is reached and the groundwater receives the same signal as the long-term precipitation average. There is no or little groundwater recharge during winter because almost all precipitation inputs arrive and accumulate as snow. Hence, we assume that the stream isotopic signature originates from groundwater only (Laudon et al., 2007; Peralta-Tapia et al., 2015). Consequently, the closer the stream signature comes to the long-term precipitation average, the more the groundwater has been mixed. The groundwater isotopic signature, in turn, should be correlated with the travel time to stream until full mixing of the precipitation signal is reached. The closer the signature is to the long-term precipitation average (which is equal to the deep groundwater measurements in Krycklan, Laudon et al., 2007), the more well-mixed and, consequently, the longer travel times will be found. We used the average annual winter signature for the evaluation. In spring, previous studies have shown that the young water fraction can be distinguished by comparing the change in the isotopic signature to the preceding winter because the snow is much lighter (depleted in ¹⁸O) (Laudon et al., 2007; Tetzlaff et al., 2015). We calculated the difference between the average winter and average spring signature for all years. The mean difference we hereafter refer to as the Δδ¹⁸O_spring, which we assumed to be negatively correlated with the young water fraction (Fig. 2a). Similarly, we refer to the mean difference between annual averages of winter and summer signature as Δδ¹⁸O_summer, which similarly should be related to the young water fraction during the summer. However, in summer, precipitation is heavier (more enriched in ¹⁸O) than in winter, which hence should give the young water a heavier signal. Therefore, we assumed a positive relationship between the young water fraction and the Δδ¹⁸O_summer (Fig. 2a).

Another indicator of travel times to stream that we used was the sum of base cation (BC) concentration (Fig. 2b) (Abbott et al., 2016). Previous attempts to follow the chemical development of groundwater in the Krycklan catchment and other streams have shown that the BC concentration increases along the groundwater flow pathway (Klaminder et al., 2011). Therefore, a correlation between the stream concentration of BCs on the one hand and modelled soil contact time on the other were assumed in this study. The BCs are mainly derived from the weathering of local soils in the Krycklan catchment, with only a minor contribution from atmospheric deposition (Lidman et al., 2014). Our assumption is further based on modelling studies of weathering rates in a soil transect in the Krycklan catchment, which indicates that there is a kinetic control of the release of BCs in the soils (Erlandsson et al., 2016). Since all BCs behave relatively conservatively in these environments (Ledesma et al., 2013; Lidman et al., 2014), we used their combined concentration as a proxy for soil contact time. However, the assumption is only valid when the water is in contact with mineral soils, not with peat in mires, which are abundant in some of the investigated sub-catchments. There are little minerals present in the peat and, therefore, the BC concentration cannot be expected to increase during the time the water spends there. Therefore, the BC concentrations were adjusted for the influence of mire, using the sub-catchment mire proportion as a scaling factor to allow a fair comparison to water–soil contact time (Lidman et al., 2014) (Table 2).

All stream chemistry data come from the online open Krycklan database (Table 2) (Krycklan Database, 2013). The isotopic signatures contain approximately 10 years of field observations (2008 to mid-2018), approximately 25 samples per year for each site. Parts of the dataset have been published by Peralta-Tapia et al. (2016), where sampling and analyses are described in detail. It has since been expanded using the same methodology. We used the average winter isotope signatures from these years as a representation of baseflow. These averages were also compared to the volume-weighted average of the long-term precipitation, calculated using approximately 1000 precipitation measurements of δ¹⁸O between 2007 and 2016. The precipitation was measured throughout the year, both as rain and as snow. The long-term precipitation average is −13.4 ‰, which is close to equal to observations of the isotopic signature at the deep groundwater wells of Krycklan (ca. 10 m depth). The BC data collection methodology is reported in Ledesma et al. (2013).

2.3 Water flow model setup

We applied the Mike SHE/Mike-11 hydrological modelling tools to quantify travel times in a pre-calculated 3D transient flow field. The simulated terrestrial hydrological system for the Krycklan catchment includes the saturated and unsaturated flow, ET, snowmelt, overland flow, and streamflow processes. The fully distributed 3D modelling tool uses topography, soil properties, and time-varying climate inputs to calculate the water fluxes throughout a catchment (Rahim et al., 2012; Sishodia et al., 2017; Wang et al., 2012; Wijesekara et al., 2014). The ET processes include canopy interception, open surface evaporation, root uptake, sublimation, and soil evaporation from the unsaturated zone based on a methodology developed by Kristensen and Jensen (1975). Flow in the saturated zone (SZ) is calculated in 3D by the Darcy equation. The flow in the unsaturated zone (UZ) is calculated in vertical 1D using the Richards equation, and overland flow (OL) is calculated using a horizontal 2D diffusive wave approximation in the Saint-Venant equations (Fig. 3). Streams are modelled in 1D using a high-order dynamic wave formulation of the Saint-Venant equations. The river model (Mike 11) is not restricted to the grid size of Mike SHE and allows for a more precise calculation of stream water levels and flow rates. The different model compartments OL, UZ, SZ, and rivers are fully integrated, and water fluxes between and within the compartments are calculated in each time step of the simulation. More in-depth documentation and manuals of Mike SHE and Mike 11 are provided by DHI (DHI, 2021).

Figure 3 Schematic of a general Mike SHE model setup. Precipitation falls on the ground as rain or snow. Evapotranspiration (ET) processes include canopy interception, open surface evaporation, root uptake, and soil evaporation from the unsaturated zone (UZ). The overland flow (OL), saturated zone (SZ), and UZ interact depending on the saturation level. The SZ is divided into 10 calculation layers (CLs), while the UZ has a much finer description. Streamflow is modelled through Mike 11 and is not restricted to the Mike SHE resolution. The figure is used courtesy of SKB. Figure illustrator: LAJ.

For the Krycklan model, the horizontal grid was set to 50 × 50 m. Vertically, the model is divided into 10 calculation layers (CLs) and extends to a depth of 100 m below ground. The SZ-CLs vary with depth and are thinner closer to the soil surface; the first CLs extend to 2.5, 3, 4, and 5 m, respectively, below the ground surface, with the soil properties and depth extension following the stratigraphy (Table 3). The UZ and SZ interact throughout the soil. If the soil is unsaturated, the UZ discretization and equations are used. The influence from ET and UZ processes on the SZ is only fully active to the depth of the uppermost SZ-CL. Here, the ET and UZ are calculated at a finer resolution, leading to a detailed calculation of the groundwater table level. The first SZ-CL depth was set to 2.5 m and was calibrated using the influence of the CL thickness on groundwater table level, UZ, and ET dynamics.

Table 3 Flow model setup. Flow model setup from the calibrated and validated Mike SHE model presented in Jutebring Sterte et al. (2018). The “soil-type surface” corresponds to the soil type shown in Fig. 1b. A drain constant was used to account for coarser material of the upper half metre of the soil.

^* The table shows the depth down to which the same description extends. For example, the first description of peat extends down to 5 m, while the first calculation layer is 2.5 m.

Following the thickness of the SZ-CL in the Krycklan model, all soils above 2.5 m depth are prescribed as one soil type, with hydraulic properties being an average of all the soil types throughout the vertical profile from the ground surface to 2.5 m depth. In Mike SHE, horizontal hydraulic conductivity (Kh) is averaged using the thickness of each soil layer. Vertical flows are more dependent on the lowest vertical hydraulic conductivity (Kv). Therefore, the harmonic weighted mean value is used to calculate the new Kv instead (Table 3). A drain function was used in this model and several previous studies (Bosson et al., 2012, 2013; Johansson et al., 2015; Jutebring et al., 2018) to account for the higher hydraulic conductivity in the uppermost part of the first CL. In the Krycklan model, the function was activated whenever the groundwater reached 0.5 m below the ground surface, above which higher K values have been observed (Table 3) (Bishop et al., 2011; Nyberg, 1995; Seibert et al., 2009). The model also accounted for soil freezing processes, which in Krycklan have been shown to have a strong influence on the water turnover in mires (Laudon et al., 2011). Based on a methodology presented in Johansson et al. (2015), soil freeze and thawing processes were described using time-varying K and infiltration capacity.

The Krycklan flow model was able to reproduce daily accumulated stream discharge, groundwater levels, and timing of precipitation events (Jutebring Sterte et al., 2018). This includes daily discharge observations (14 streams) and weekly to monthly observed groundwater levels (15 wells) for 2009–2014. The accumulated error in stream discharge was on average 11 % and highest for sub-catchments with few observation points (<25 %). For this study, a few changes were made to the Krycklan model (Jutebring Sterte et al., 2018). Most importantly, new field data from the Krycklan database gave a more precise location of the observation station at C5 (red circle in Fig. 1). The Kh of silt was also increased from 1 × 10⁻⁸ m/s to 1 × 10⁻⁷ m/s due to new soil property samples, which gave a slightly better flow representation of the sites' affected silty sediments. However, the corrections and additions did not influence the model results in any substantial way. The improvements were small but were made to better represent the site with the hydrological flow model (Jutebring Sterte et al., 2020).

2.4 Establishing travel times – particle tracking

Particle tracking in Mike SHE enables investigations of groundwater travel time from the recharge to the SZ until the discharge into the streams, as described in detail in Bosson et al. (2010, 2013). The model calculates the location and age of separate particles added with infiltrating water along their flow lines. The particles move by advection governed by the pre-calculated groundwater flow field from the Mike SHE model (Jutebring et al., 2018, 2021). This method allows for long-term transport calculations where particle tracking can be run for several annual cycles based on the same transient or steady-state flow field. The advection–dispersion equation governs the transportation of particles for a porous medium. The Darcy velocity is divided by the porosity to calculate the groundwater velocity. Therefore, the only complementary input data needed to run the particle were porosity values (Table 4).

Table 4 Porosity values for different soil types used in the Mike SHE model.

^a Average of Morris and Johnson (1967). ^b Joyce et al. (2010). ^c Average value between sand and clay. ^d Average value between silt and clay.

Particle tracking was used to assess groundwater travel times from groundwater recharge to stream runoff for each sub-catchment. The model was run for 1000 years to capture the travel times of all discharging groundwater for each sub-catchment. One year of simulated flow results was cycled 1000 times to extend the particle-tracking simulation. The year 2010 was selected, as the water balance was close to the long-term annual averages observed for the Krycklan catchment. All particles were released at the top of the transient groundwater table the first year. Numerical constraints restricted the number of particles released to 0.5 particles/10 mm modelled groundwater recharge per grid cell, which corresponds to a total of approximately 0.6 million particles for the entire modelled area in the first year. This number of particles was assumed to be enough to capture the timing of recharge patterns (Fig. 4).

Figure 4 Particle model setup. (a) Steps of particle tracking. (b) Average depth to the groundwater table. The main part of the model area has a calculated depth to the groundwater table between 0 and 3 m and varied daily. (c) Schematic illustration of particle-tracking setup. Particles were added to groundwater recharge at the transient groundwater table. The age of these particles was zero at the time of recharge. Thereafter, they followed the groundwater flow, increasing in age until reaching a stream or lake.

2.5 Analysis of modelled travel times, relationship with stream chemistry and landscape characteristics

The time it took for particles to reach a stream or lake via groundwater (hereafter called “travel time”) was calculated for each sub-catchment. The calculated travel time distributions were based on all particles arriving in a stream within a certain period of time, either annually or for a specific season, for the entire modelling period. The distributions were analysed using four statistical measurement tools, the arithmetic mean, the geometric mean, the median, and the standard deviation (SD). The arithmetic mean, the geometric mean, and the median are common choices to describe the central tendency of a distribution (Destouni et al., 2001; Kaandorp et al., 2018; Massoudieh et al., 2012, 2017; Unlu et al., 2001), which all have their strengths and weaknesses. If the distribution is not significantly skewed, the SD is smaller than half of the average (Taagepera, 2008). In the case of the observed δ¹⁸O and BC concentrations (Table 2), the SD is much smaller than half of the average. Therefore, the arithmetic mean was used to describe the central tendency of the data set. However, if the travel time distribution becomes skewed, the arithmetic mean becomes highly sensitive to the tail of the distribution and produces considerable uncertainty. In these cases, the median and the geometric mean are often better as a measure of the central tendency of mean travel time (MTT) than the average. However, to compare the MTT of discharged water of different streams, we still wanted the metric to account for the length of the tail. Therefore, we used the geometric mean because the median only states the middle value of a distribution regardless of the tail length (Taagepera, 2008; Unlu et al., 2004; Zhang and Zhang, 1996). However, we provide all metrics, including the arithmetic mean, geometric mean, median, and SD, in the Appendix, Table A2.

The MTT was compared to stream chemistry, which is a mix of both groundwater and surface water. In winter, all streamflow contributions originate from groundwater. Here the results from the particle tracking reflect the actual travel time to the streams. However, in summer and especially in spring, some water will reach the streams via overland flow (OL), which has not spent any time in the ground. Since the particle tracking does not take surface flow into account, OL was accounted for by reducing the MTT by using the OL fraction as a scaling factor (Appendix, Table A2). The young water fraction was also used as an evaluation criterion. Like previous studies (Kirchner, 2016; von Freyberg et al., 2018; Lutz et al., 2018; Stockinger et al., 2019), we assumed young water fraction to be the sum of all water less than 3 months old. In our case, this includes all water reaching streams as overland flow and as young groundwater (< 3 months). The modelled MTT and young water fraction were also used to identify the main factors determining the travel times to stream. The catchment characteristics tested included important terrain factors such as catchment size, slope, and main soil types (Table 1).

3 Results

3.1 Travel time results

The particle-tracking results were used to establish travel time distributions and MTT of water to the streams of the 14 sub-catchments in Krycklan. Since the travel time distributions were significantly skewed, we assumed that the geometric mean of the travel time distributions provided the best representation of MTT (Table 5, Fig. 5). However, all metrics are stated in the Appendix, Table A2. The annual MTT_geo for all sub-catchments ranged from 0.8 to 3.1 years (Table 5). Most groundwater discharging to a stream had a travel time of less than 1 year in all sub-catchments (34 % to 54 %). The longest stream MTTs were connected to the larger catchments, such as C16, and the silt-dominated catchments such as C20. We used some sub-catchments for result representation, but all results are provided in Table 5 and Appendix Table A2. The displayed sub-catchments were C2 (small till- and forest-dominated catchment), C4 (small mire-dominated catchment), C20 (small silt-dominated catchment), and C16 (the full-scale Krycklan catchment).

Table 5 Annual and seasonal (winter, spring, and summer) travel times. The geometric mean of the travel time distribution (MTT_geo) is adjusted for the overland flow. The young water fraction (YWF) includes overland flow and groundwater younger than 3 months (%). An extended version of the results, including arithmetic mean, median, and SD, is included in the Appendix (Table A2).

Figure 5 Examples of particle-tracking results. The figure shows the distribution of all particles reaching the different streams for the entire modelling period. The solid line shows the statistics for C16, including the 25th percentile, the median, the geometric mean, the arithmetic mean, and the 75th percentile (Appendix A). Moreover, the figure shows three other example distributions, including C2 (small forest- and till-dominated catchment), C4 (small mire-dominated catchment), and C20 (small silt-dominated catchment).

On an annual basis, a fraction of water reached the streams as overland flow. A major part of the overland flow occurred during the snowmelt in spring, especially in sub-catchments with mires such as C4 (Fig. 6). Both the fraction of young water reaching the streams and the MTT_geo displayed strong seasonal trends. The longest seasonal MTT_geo, 1.2–7.7 years, and the smallest young water fraction were found during the winter season. In winter, the fraction of older water successively increased until the spring snowmelt began in early April. Conversely, the smallest fraction of old discharging water and short MTT_geo, 0.5–1.9 years, were connected to events of larger groundwater recharge, such as the spring snowmelt and heavy summer rains.

Figure 6 Seasonal fraction of discharge to streams. The figure shows the proportion of annual stream discharge arriving as groundwater and overland flow. Four sub-catchments are exemplified, including (a) the small till- and forest-dominated C2, (b) the small mire-dominated C4, (c) the silt-dominated C20, and (d) the full-scale Krycklan catchment C16 with mixed mires and forests (extended version in Appendix Fig. A1). The figure showcases the travel time fraction of water discharging to the streams. The fractions are both shown as part of the total annual discharge as well as the water composition. The bands below the months highlight the three investigated seasons, spring, summer, and winter.

In spring, mire sub-catchments had the shortest MTT_geo. However, as exemplified by the similar-sized C2 and C4 sub-catchments, groundwater was not renewed to the same extent in mire-dominated systems due to a larger fraction of surface runoff (Fig. 6). Mire-dominated sub-catchments (like C4) displayed stronger seasonal variations in MTT_geo, with shorter MTT_geo than till-dominated sub-catchments (like C2) in spring and longer MTT_geo than C2 in winter (Table 5). In C4, the MTT_geo decreased from 1.5 to 0.7 years from winter to spring, while the corresponding change in C2 was 1.2 to 0.7 years. The seasonality of MTT_geo was even more pronounced for catchments with a larger areal coverage of mires combined with a larger areal coverage of silt. For example, C20 had an MTT_geo that decreased from 7.7 to 1.9 years from winter to spring (Table 5).

3.2 Testing model results on stream isotopic composition and chemistry

In addition to investigating the annual MTT_geo, three distinct seasons were evaluated regarding the stream chemistry: winter, spring, and summer. The isotopic composition was available for 13 out of 14 sub-catchments (C20 excluded because of short time series), while the BC data were available for all sites. In winter, the modelled MTT_geo was correlated with the isotopic composition ($r=-\mathrm{0.80}$, P<0.01), with older stream water being closer to the long-term precipitation average (Fig. 7a). Some of the sub-catchments had an isotopic signature close to the precipitation average, suggesting almost complete mixing (e.g. C16). The negative correlation between the Δδ¹⁸O_spring and the young water fraction was also significant ($r=-\mathrm{0.90}$, P<0.0001, Fig. 7c), following the conceptual model (Fig. 2a). The same was also true for the summer season but with a weaker positive correlation compared to the spring (r=0.80, P<0.001, Fig. 7e), again agreeing with the conceptual model (Fig. 2b). The opposite sign of the slope was due to the heavier summer precipitation compared to the winter baseflow. The correlation between the BC concentration and MTT_geo was strongest in winter (r=0.90, P<0.0001) and weakest in summer (r=0.79, P<0.001). The sub-catchments with the longest travel times to streams and highest BC concentration included the sub-catchments with larger areal coverage of silt, for example, C16 and C20. The shortest travel times and lowest BC concentrations were connected to smaller sub-catchments in till-dominated areas, such as C2 and C4.

Figure 7 Relationships of seasonal MTT_geo and young water fractions (YWFs) with seasonal stream isotopic composition and base cation concentration. Note that δ¹⁸O results are for 13 sites, while the BC record comprises all 14. The sub-plots (a) to (f) show the δ¹⁸O (winter) or Δδ¹⁸O_{spring/summer} and BC concentrations as a function of the MTT_geo in winter, spring, and summer, respectively. The standard error of the mean (SEM) shown as whiskers denotes variations in field observations.

3.3 Model results compared to catchment characteristics

The main catchment characteristics correlated with MTT_geo and young water fraction were catchment size, the areal coverage of low-conductive silty sediments, and the areal coverage of mires. The strongest positive correlation was found between the young water fraction and the areal coverage of mires (r=0.96, P<0.0001). There was also a strong positive correlation between MTT_geo and the areal coverage of silt (r=0.90, P<0.0001) (Fig. 8). A positive correlation between catchment size and MTT_geo was also found, albeit weak due to one catchment, C20, yet significant (r=0.63, P<0.05) (Fig. 8). However, the catchment size was also correlated with the areal coverage of silt, which may be the underlying reason for this correlation (Table 6), as C20 is the only relatively small monitored sub-catchment located in the area with sorted sediments. The annual and seasonal patterns were similar (Table 6). However, the positive correlation between mires and the young water fraction was lost in winter due to a lack of new precipitation input into the system. A weak negative correlation between MTT_geo and the young water fraction was found for the annual and spring seasonal results but was lost for the summer and winter.

Figure 8 Catchment characteristics are important for travel times. The figure shows the annual averages: (a) the areal coverage of mires and the young water fraction (YWF), (b) areal coverage of silt and MTT_geo, and (c) catchment size and MTT_geo.

Table 6 Correlation matrix – young water fraction (YWF), geometric mean travel time (MTT_geo), and catchment characteristics. The catchment characteristics include the log catchment size (log A), the areal coverage of mires (Mire), and the areal coverage of silt (Silt). The table includes annual, winter, spring, and summer results.

For $\left|r\right|\mathit{>}\mathrm{0.5}$, the p value is shown according to ^a p<0.05 and ^b p>0.05.

4 Discussion

Particle tracking in the Mike SHE model provided valuable insights into the annual and seasonal mean travel times (MTT_geo) across the 14 Krycklan sub-catchments. The modelled MTT_geo and the young water fractions were strongly correlated with observed stream δ¹⁸O_winter signatures, seasonal variation in δ¹⁸O, and base cation (BC) concentrations. This model validation suggests that particle tracking is a useful complementary tool to tracer-based travel time studies, at least in snow-dominated catchments, areas with pronounced seasonality, and streams dominated by older groundwater (> 4 years). Overall, we found that soil type was the most important variable explaining MTT_geo and that mires are an important landscape feature regulating the young water fraction in spring (Fig. 8).

4.1 Model assumptions and limitations of estimated travel times

Comparing the results from this modelling study to previous Krycklan investigations of MTT conducted in the C7 sub-catchment demonstrates that different model approaches have provided similar results. While our study suggested a MTT_geo of 1.1 years and a median of 0.8 years (Appendix, Table A2), Peralta-Tapia et al. (2016) calculated a MTT of 1.8 (minimum 0.8 and maximum 3.3) years by applying a mathematical method for isotopic dampening to fit a model to the observed stream isotopic response. In another recent study using the Spatially distributed Tracer-Aided Rainfall-Runoff (STARR) model for the same stream, the median age was estimated to 0.9 years (Ala-aho et al., 2017). The close agreement between the different studies strengthens the overall reliability of the results. However, like all modelling techniques, particle tracking in Mike SHE is associated with some uncertainties and limitations.

In contrast to the Mike SHE flow model, which estimates groundwater and overland flow pathways, the particle-tracking model is restricted to the subsurface hydrological component. This is a limitation in the modelling approach as water reaches the streams as a mix of groundwater and overland flow. Therefore, to allow for actual MTT_geo estimates, we corrected the results by reducing the estimated MTT_geo using the overland flow from the flow model as a scaling factor (Appendix, Table A2). This uncertainty primarily affects the mire-dominated sub-catchments that have a large fraction of overland flow, especially during the spring.

Another uncertainty related to the particle-tracking model in Mike SHE is related to the travel time from the point of infiltration through the unsaturated soil horizons to the saturated groundwater. Due to technical limitations, this travel time cannot be accounted for in the particle-tracking calculations. Particles are placed at the groundwater table proportionally to the groundwater recharge (Fig. 4). Therefore, the main fraction of particles introduced to the model occurs at high infiltration rates when the groundwater level is close to the soil surface. Under these conditions, the water has, in most cases, spent a relatively short time in the unsaturated zone. However, some particles are also introduced when the groundwater level is lower, such as early snowmelt or following extended dry periods. Under such conditions, the model uncertainty increases. In this context, the smallest potential uncertainty occurs in mires, where the groundwater table always is close to the ground surface. The uncertainty becomes somewhat larger in the till areas where the unsaturated zone on average is above 1 m but can extend down to 3 m below the ground during low flow. C14 and the lower part of C16 are exceptions to these relatively shallow saturated conditions as a deep esker traverses the sub-catchments resulting in a groundwater level up to 10 m below the soil surface (Fig. 1). Accounting for the travel time from infiltration to recharge could impact the results and provide, especially for C14 and C16, longer MTT than if the groundwater level were at the same level throughout the whole catchment. This limitation primarily affects catchments with the longest MTTs and, therefore, does not seriously question the general pattern observed. The distance from the ground surface to the groundwater table is, for most model cells, much shorter than the distance to the nearest stream, so most of the transit time should be related to the groundwater flow rather than to percolation. Although water, especially during dry conditions, no doubt can spend considerable time in the unsaturated zone, it must also be acknowledged that this water volume is small compared to the groundwater inventory in the saturated zone. Therefore, its impact on the average MTTs should be relatively small.

4.2 Seasonality of isotopic composition

Following the conceptual model (Fig. 2), patterns in stream isotopic signatures can be explained by seasonal changes in travel times. The modelling results show that all sub-catchments discharged water with the longest travel times in winter, somewhat shorter travel times in summer, and water with the shortest travel times in spring. When winter arrived, the main precipitation was snow, resulting in a cessation of the groundwater recharge. This caused an increasing proportion of old groundwater discharging into the streams (Fig. 6). In agreement with our conceptual model (Fig. 2), a strong negative correlation between winter MTT_geo and the isotopic stream signatures during winter baseflow was observed (Fig. 7a). At an average travel time older than 4 years, it can be expected that the groundwater has reached full mixing. Hence, older water can no longer be accurately quantified using amplitude dampening of the water isotope signal (Kirchner, 2016). These theoretical considerations strengthen the results of a winter MTT_geo older than 4 years for some sub-catchments since their stream isotopic signatures were close to the long-term precipitation average and, therefore, should have reached complete mixing.

When snowmelt began in late April or early May, the MTTs consistently decreased in all sub-catchments. The fraction of young groundwater in different sub-catchments was well reflected in the change in the isotope signal (Fig. 7). For snowmelt in spring, the calculated young water fraction was used to evaluate the proportion of water reaching the stream through rapid pathways, including overland flow. It is well established that the difference in stream isotopic signature between winter baseflow and spring peak flow at snowmelt (Δδ¹⁸O_spring) is mechanistically linked to the amount of new water reaching the stream (Tetzlaff et al., 2009). In agreement with this, we found a strong statistical relationship between Δδ¹⁸O_spring and the calculated young water fraction (Fig. 7c). These results are well in line with previous work in Krycklan using end-member mixing of new and old water in the same streams (Laudon et al., 2004, 2007, 2011).

Similarly to the conditions in spring, the conceptual model predicted that the difference in stream isotopic signature between winter baseflow and summer flow, Δδ¹⁸O_summer, should be correlated with the young water fraction in summer but with the opposite sign, due to isotopically heavier summer rains (Fig. 2). A larger inter-annual variation in precipitation and high ET likely caused the relationship to be less evident compared to the spring results as the snowmelt conditions are more consistent from year to year. The groundwater signal reaching the streams during the summer season may also be affected by a lingering signal from the snowmelt. However, although less evident than compared to the Δδ¹⁸O_spring, there was still a significant correlation between the average Δδ¹⁸O_summer and the modelled young water fraction (Fig. 7e).

4.3 Controls of travel times on base cation concentrations

The annual and seasonal average BC concentrations were positively correlated with the MTT_geo (Fig. 7b, d, and f). Since the weathering rates were assumed to be kinetically controlled and hence related to the exposure time of water to minerals, spatial and temporal variability in BCs can be used as a relative indicator for transit time (Erlandsson Lampa et al., 2020). However, reducing weathering to travel times may be an oversimplification as the rate is also affected by differences in mineralogy, particle size distributions and the chemical conditions in the groundwater. However, previous research in the Krycklan catchment has suggested that the chemical composition of the local mineral soils is surprisingly homogeneous, even when comparing till and sorted sediments (Klaminder et al., 2011; Peralta-Tapia et al., 2015; Erlandsson et al., 2016; Lidman et al., 2016). Therefore, we did not expect mineralogical differences between soil types to impact the release of cations significantly. However, one exception is peat deposits, which strongly affect the cation concentrations on the landscape scale. The effect of the peat was accounted for by adjusting the concentrations following Lidman et al. (2014). Differences in particle size distribution may be important because coarser soils will have less surface area per volume unit, allowing for less weathering. However, such soils can also be expected to have higher hydraulic conductivities, leading to higher flow velocities and, consequently, less time available for weathering. Therefore, differences in area–volume ratios between different soil types would not counteract the effect of travel times on the weathering but rather enhance it.

Despite arguments that can be made against the use of BCs as tracers, they still offer a complementary possibility of testing the model performance (Abbott et al., 2016). As McDonnell and Beven (2014) emphasized, the inclusion of tracers in hydrological models is necessary to ensure that a model reproduces the speed of flow, which is an important parameter when assessing travel time distributions. For catchment-scale models, this could be an isotopic tracer or a solute transported with the water (Hooper et al., 1988; Seibert et al., 2003; Fenicia et al., 2010; Hrachowitz et al., 2013). Although neither the travel time distribution nor the kinetics of weathering is fully understood, the strong agreement between the calculated travel times and the observed stream water chemistry provides additional support that our modelling of these processes – and thus the entire system – was reasonable and consistent with the empirical data.

4.4 Mean travel times, young water fractions, and catchment characteristics

All sub-catchments showed similar seasonal patterns in MTT_geo and young water fraction, manifested as water with shorter travel times discharging in spring and water with longer travel times discharging in winter. Some of the catchment characteristics influenced the magnitude of these seasonal patterns across the landscape. On a landscape level, the main causal mechanism determining the annual MTT_geo was the areal coverage of silty sediments (Table 1), which largely overshadowed the importance of other catchment characteristics (Fig. 8, Table 6). This finding partially stands in contrast to earlier studies in Krycklan by Peralta-Tapia et al. (2015) and Tiwari et al. (2017) that suggested that the groundwater travel times are nonlinearly linked to the catchment size. However, a correlation of travel times with catchment size appears to be a spurious relationship since there is a correlation between the catchment size and the areal coverage of silty sediments (Table 6). The reason is that the silty sediments are located in the lower parts of the Krycklan catchment, which implies that all large catchments contain at least some proportion of silty sediments. The importance of the silt, rather than the catchment area, for the groundwater travel time is most clearly illustrated by the small silt-dominated C20 catchment, which was a distinct outlier to such a scale-dependent pattern. This indicates that catchment size may not be the primary factor determining the variability of travel times of different catchments (Fig. 6). Instead, the long travel times in C20 suggest that the groundwater flow velocity is slower in the silt areas than elsewhere in Krycklan, even though the average catchment slope is steeper than comparably sized sub-catchments in till areas (Table 1 and Fig. 1).

Similarly, silt may also explain the long travel times at C14 and C16. Although C14 is smaller than C15, which mostly lacks silt, C14 still has a longer MTT_geo. In contrast, MTT_geo in C15 is much closer to C12 and C13, even though the C15 catchment is twice the size (Table 5). Hence, the results suggest that the critical difference between these sub-catchments and other sub-catchments is related to the soil hydraulic conductivity rather than the catchment size. The results further emphasize that one cannot assume that the travel time would increase with catchment size unless the distribution of different soils is comparable throughout the landscape. A larger catchment size does not necessarily imply that considerably more amounts of old groundwater are discharged directly to their streams. Instead, it may just be mainly a confluence of stream water from smaller sub-catchments that determined variation, and in that case, no relationship between catchment area and groundwater travel times should be expected.

The effect of the areal coverage of silty sediments was especially prominent in winter when the range in MTT_geo is between 1 and almost 8 years. The change in seasonal MTT_geo from winter to spring was also largest for the silt-dominated catchments. For example, there was a 6-year difference for C20 compared to the 2-year difference for the similar-sized till-dominated sub-catchment C6. These intra-annual variations can also be linked to another landscape feature, namely the areal coverage of mires. Mires affected the young water fraction only when new precipitation or snowmelt input into the system occurred in spring and summer. The contrasting hydrological response of mire and silt areas, respectively, caused greater annual MTT_geo variation for sub-catchments with both features. For example, the MTT_geo for C4, dominated by mires, decreased from 1.5 to 0.7 years from winter to spring. In contrast, winter MTT_geo for the C20 catchment dominated by silt was 7.7 years, which decreased to 1.5 years in spring. The results also show that groundwater recharge is affected by the soil frost in mires. For example, C4 showed more variations in its seasonal MTT_geo, although C2 (dominated by forest and till) and C4 (dominated by mires) had almost an equal annual MTT_geo (Table 5). In spring, the MTT_geo was shorter in C4 than in C2 due to the surface runoff from the frozen mire. Besides the slightly higher specific discharge (Karlsen et al., 2016), empirical studies suggest that the soil frost on mires causes a large fraction of overland flow (Laudon et al., 2007, 2011). Looking only at the part of runoff originating from groundwater (Appendix, Table A2), the MTT_geo for C4 decreases from 1.5 to 1.2 years. However, C4 still showed more seasonal variation than C2, due to longer travel times in winter. The results suggest that the lack of recharge during spring affects the renewal of the groundwater in mire-dominated catchments have long-term consequences for travel times.

Earlier studies have demonstrated that fluxes of old groundwater are more stable throughout the year than younger groundwater, which shows a more variable temporal pattern (Rinaldo et al., 2011; van der Velde et al., 2015; Kaandorp et al., 2018). In this system, such a pattern can mechanistically be linked to the till soils that dominate most sub-catchments. The response of groundwater in till soils to precipitation events can be described by transmissivity feedback processes (Bishop, 1991), which are caused by the fact that the hydraulic conductivity increases exponentially towards the soil surface. When water infiltrates the ground, the water table rises and activates more conductive soil layers, resulting in a rapid increase in the lateral flow. This implies that much of the water transport in till soils occurs relatively close to the surface. Simultaneously, the groundwater in deeper layers is more stagnant, further explaining the relatively short and consistent MTTs of till areas. Measurements of chlorofluorocarbons (CFCs) further support deeper groundwater transport in till soils in Krycklan being slow. Not far below the groundwater table, CFCs have indicated that the groundwater can be several decades old, suggesting that most groundwater transport occurs close to the surface (Kolbe et al., 2020). Consistent with this explanation, silt-dominated areas had much longer MTTs than comparatively sized sub-catchments underlain by till soils because till soils have a greater decline in hydrologic conductivity with soil depth (Figs. 6, 8 and A1).

5 Conclusions

The combination of stable water isotopes, stream water chemistry, and particle tracking provided a consistent picture of the hydrological functioning of a boreal catchment and which processes and factors are most important for regulating water pathways and travel times. We identified specific landscape characteristics that impact the seasonal distribution of travel times by combining a distributed hydrological model with empirical observations from 14 nested sub-catchments. In the wake of a changing climate and intensified pressure from forestry and other types of land use, this study provides a useful baseline for assessing the intricate connections and feedbacks between hydrological and biogeochemical processes throughout the boreal landscape. Our results showed that travel times to stream could vary considerably on annual and seasonal scales between different types of catchments. This was mainly related to soil properties, with low-conductivity silty sediments leading to the longest travel times on annual and seasonal timescales. In contrast, mires led to increased fractions of young water and hence shorter travel times, but mainly in spring when the soil was frozen. However, mire-dominated catchments experience longer travel times than similar-sized forested catchments in winter. Generally, for the boreal landscape, a warmer climate is predicted with reduced snow cover and snow duration, accompanied by increases in the frequency of winter thawing episodes and reduction in soil frost (IPCC, 2014; Jungqvist et al., 2014; Brown et al., 2017; Lyon et al., 2018). Our results suggest that these changes would reduce the intra-annual variations of MTT created by the freezing of mires, while the impact on other parts of the landscape would remain relatively low.

Appendix A

The Appendix includes the classification of the seasons for each year (Table A2). It also includes the statistical measurements for particle tracking and isotope data used in e.g. Table 5 and Table A2. The statistics used in this study include arithmetic mean (Eq. A1), geometric mean (Eq. A2), standard deviation (Eq. A3), and standard error of the mean (Eq. A4). The isotopic signature of δ¹⁸O has been calculated as Eq. (A5). Equation (A6) shows the equation used for adjusting the MTT from the particle tracking (MTT_{particle tracking}) to the overland flow MTT shown in Table 5. Equation (A7) shows the lake regression from Peralta-Tapia et al. (2015), adjusted according to new data.

$$\begin{array}{}\text{(A1)}& {\displaystyle }\begin{array}{rl}& \text{Arithmetic mean of the travel time distribution}=\\ & \left({\displaystyle \frac{\mathrm{1}}{n}}\right){\sum }_{i=\mathrm{1}}^n\mathrm{ai}\end{array}\text{(A2)}& {\displaystyle }\begin{array}{rl}& \text{Geometric mean of the travel time distribution}=\\ & {\mathrm{10}}^{\left(\frac{\mathrm{1}}{n}\right){\sum }_{i=\mathrm{1}}^n\log \left(\mathrm{ai}\right)}\end{array}\text{(A3)}& {\displaystyle }\mathrm{SD}=\text{standard deviation}=\sqrt{{\displaystyle \frac{\sum (\mathrm{ai}-\mathrm{aMTT}{)}^{\mathrm{2}}}{N}}}\text{(A4)}& {\displaystyle }\mathrm{SEM}=\text{standard error of the mean}={\displaystyle \frac{SD}{\sqrt{n}}}\text{(A5)}& {\displaystyle }\mathit{\delta }=\left({\displaystyle \frac{R_{\mathrm{sample}}}{R_{\mathrm{standard}}}}-\mathrm{1}\right)\mathrm{‰},R={\displaystyle \frac{{}^{\mathrm{18}}\mathrm{O}}{{}^{\mathrm{16}}\mathrm{O}}}\text{(A6)}& {\displaystyle }\mathrm{MTT}={\mathrm{MTT}}_{\text{Particle tracking}}(\mathrm{1}-\text{fraction OL})\text{(A7)}& {\displaystyle }\begin{array}{rl}{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}& =\mathrm{0.2}\left(\text{lake coverage,\%}\right)\\ & -\mathrm{13.2}(R^{\mathrm{2}}=\mathrm{0.9},p\mathit{<}\mathrm{0.01})\end{array}\end{array}$$

ai: data set values; n: number of values; OL: overland flow.

The start and end dates used in this study for the evaluation of isotopic signatures and BC are shown in Table A1 (Krycklan Database, 2013). The winter season generally occurs from early December to late February, from negative air temperatures until increasing temperatures cause small snowmelt events. Here, the isotopic signal is more sensitive to snowmelt than the BC. Spring is the period of the main part of the spring flood, mainly occurring in April to May. Moreover, the summer season includes all observations occurring between July and September.

Table A1 Dates used for chemistry investigation. The dates are the start and end dates for observations within the seasons classified as winter, spring, and summer, respectively. Note that BC only includes dates from 2008 to 2016.

The Appendix also includes an extended version of Table 5 and Fig. 6, including all sub-catchments (Table A2 and Fig. A1). The figure shows the travel time fraction of water reaching the streams of the sub-catchment in Krycklan annually. The figure shows both the groundwater fraction (travel time fraction calculated using the particle-tracking results) and the simulated direct runoff fraction. The table shows more statistical information regarding the travel time distribution, including the skew, SD, and SEM and the fraction of overland flow on an annual and seasonal timescale.

Figure A1 Travel time fraction of water discharging to the streams of Krycklan, in increasing annual geometric mean travel time, MTT_geo (a–n). The black vertical line shows the 50 % mark for visual aid. All charts begin in spring (late March/early April) and end in winter (early March).

Table A2 Extended version of Table 5 – annual and seasonal (winter, spring, and summer) travel time results. The table includes the travel time results based on particle tracking before taking the overland flow into account. The results include arithmetic mean (A), median (M), geometric mean (Geo), skew, standard deviation (SD), and standard error of the mean (SEM). The table also shows the fraction of overland flow (OL) used as a scaling factor for the calculation of MTT.