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 subcatchments of the Krycklan catchment (0.12–68 km2). 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 (MTTgeo) 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 MTTgeo ranged from 1.2 to 7.7 years, while spring MTTgeo varied from 0.5 to 1.9 years. The modelled variation in annual and seasonal MTTgeo and the fraction of young water (<3 months) was supported by extensive observations of both δ18O 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 MTTgeo, 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 MTTgeo 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.

• Soil type was the most important factor regulating the variation in mean travel times between different sub-catchments • The areal coverage of mires increased the young water fraction, especially during the spring.
• The greatest seasonality in mean travel times was found in catchments dominated by silty soils because of the long travel times in winter and relatively short travel times in spring.

40
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 have on hydrological and biogeochemical processes 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 mobilisation of organic carbon and associated solutes (Tiwari et al., 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 four to five 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 two to three months, can, however, still be quantifiable 70 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 take from the source to the recipient stream .
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, needas far as possibleproper 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). CollectionThe collection of such field data is, however, costly and time-consuming. Therefore, data for calibration and validation is often limited, and the minimum length and types in data-sparse catchments is currently a topic of increasing interest 90 (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 longlasting snow cover that often melts rapidly in the spring create both opportunities and challenges for determining the pathways and 95 travel time of water discharging to streams. The long, snow-rich winters do 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), 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 100 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.

105
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 110 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.

115
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 ten-year observational records from Krycklan, including average seasonal 120 changes of 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 subcatchments, 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).

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 featuresvegetation of this boreal landscape areis 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.
Krycklan has a landscape distinctively formed by the last ice age (Ivarsson and Johnsson, 1988;Lidman et al., 2016). At the higher 140 elevations to the northwest, located above the highest postglacial coastline, the till-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. In this study, we refer to soil as all unconsolidated material aboveThe deeper soils consist of basal till which was deposited and compacted under the bedrock.
Themoving ice. In contrast, the shallower till layers consists 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 145 characteristic for glacial tills in northern Sweden (Bishop et al., 2011;Nyberg, 1995;Seibert et al., 2009) with conductivities close to 5*10 -5 m/s at the ground surface and exponentially decreasing downwards (Nyberg, 1995). 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.

150
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 155 differences in topography (Table 1). 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, 2014, 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.).

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. We defined the winter as late November to late March because it (Table 2, Table 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 zero °C, until the air temperature started to rise above freezing temperatures again, causing snowmelt. This season is characterised 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 characterised by rain, high ET, and relatively little runoff. March, June and, October, and November were excluded because, hydrologically, they are typically transition months between the three distinct seasons. This is because snowmelt can still influenceinfluences runoff in March and June,. October and November are transitional months between summer and winter 180 conditions (with irregularly occurring snowfall, and soil frost, etc.) can sometimes begin to establish in already October events.
In this study, stable water isotopes (δ 18 O) were used to track pathways of precipitation inputs to stream networks (see Appendix for the δ 18 O definition). Ten years of δ 18 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 185 was corrected by accounting for the percentage of lakes in each sub-catchment (Table 2) The comparison of the modelling results to observations of δ 18 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 190 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 accumulates as snow. Hence, we assume that the stream isotopic signature originates from groundwater only (Laudon et al., 2007;Peralta 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 195 to 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 wellmixed 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 18 O) (Laudon et al., 2007;Tetzlaff et al., 2015).

200
We refer tocalculated the difference between the average winter and average spring signature for all years. The mean difference we hereafter refer to as the Δδ 18 Ospring, which we assumed to be negatively correlated to the young water fraction (Fig. 2a).
Similarly, we refer to the mean difference between the averageannual averages of winter and average summer signature as Δδ 18 Osummer, which similarly should be related to the young water fraction during the summer. However, in summer, precipitation is heavier (more enriched in 18 O) than in winter, which hence should give the young water a heavier signal. Therefore, we assumed 205 a positive relationship between the young water fraction and the Δδ 18 Osummer (Fig. 2a).
Another indicator of travel times to stream that we used was the sum of base cations (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 210 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 215 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).    (2016), where sampling and analyses are 235 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 1160approximately 1000 precipitation measurements of δ 18 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 240 collection methodology is reported in Ledesma et al. (2013).  sublimation, and soil evaporation from the unsaturated zone based on a methodology developed by Kristensen and Jensen (1975).

265
For the Krycklan model, the horizontal grid was set to 50*50 m. Vertically, the model is divided into ten calculation layers (CL) 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 m, 3 m, 4 m, 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 discretisation 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.
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.

275
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, Johansson et al. 2015and 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 280 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 . Based on a methodology presented in Johansson et al. (2015), soil freeze and thawing processes were described using time-varying K and infiltration capacity. 290 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 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 -8 m/s to 1*10 -7 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 according to all available data.with the hydrological flow model (Jutebring Sterte et al., 2021). Fine till 1*10 -6 1*10 -7 Bedrock 1*10 -9 1*10 -9 Peat 5 Peat 1*10 -5 5*10 -5 7 Clay 1*10 -9 1*10 -9 To bedrock Fine till 1*10 -6 1*10 -7 Bedrock 1*10 -9 1*10 -9 Silty sediments 3 Silt/clay 1*10 -7 1*10 -7 To bedrock Fine till 1*10 -6 1*10 -7 Bedrock 1*10 -9 1*10 -9 Sandy Sediments 4 Silt/Sand 2*10 -5 2*10 -5 0.9*max depth Sand 3*10 -4 3*10 -5 To bedrock Gravel 1*10 -4 1*10 -4 Bedrock 1*10 -9 1*10 -9 Drain constant Peat 1*10 -6 Till 4*10 -7 Silty sediments 1*10 -7 a 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.

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 310 into the streams, as described in detail in Bosson et al. (2010Bosson et al. ( , 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(Jutebring et al., , 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 advectiondispersion equation governs the transportation of particles for a porous medium. The Darcy velocity is divided by the porosity to 315 calculate the groundwater velocity. Therefore, the only complementary input data needed to run the particle was porosity values (Table 4).
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 from Jutebring Sterte et al. (2018) 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).   Kaandorp et al., 2018;Massoudieh et al., 2012Massoudieh et al., , 2017Unlu et al., 2004), 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 δ 18 O and BC concentrations (Table 2), the SD is much smaller than half of the average. Therefore, the arithmetic mean was used 340 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 345 (Taagepera, 2008;Unlu et al., 2004;Zhang et al., 1996). However, we provide all metrics, including the arithmetic mean, the geometric mean, median, and SD, in the Appendix, Table A1A2.
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.

350
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 A1A2). The young water, 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 three month old. In our case, this includes all water reaching streams as overland flow and as young groundwater (<three 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).
• Particles were removed when they reached a sink, which in most cases, was the stream. However, sinks could also be the unsaturated zone, dried out cells (due to evapotranspiration), or if a particle left the model through the model boundaries.

Travel time results
The particle tracking results were used to establish travel time distributions and MTT of water to the streams of the 14 subcatchments in Krycklan. Since the travel time distributions were significantly skewed, we assumed that the geometric mean of the 370 travel time distributions provided the best representation of MTT (Table 5, Fig. 5). However, all metrics are stated in the Appendix, Table A1A2. The annual MTTgeo 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 one 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 A1. The displayed sub-catchments were: C2 (small till and 375 forest dominated catchment), C4 (small mire dominated catchment), C20 (small silt dominated catchment), and C16 (the full-scale Krycklan 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 385 streams and the MTTgeo displayed strong seasonal trends. The longest seasonal MTTgeo, 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 MTTgeo, 0.5-1.9 years, were connected to events of larger groundwater recharge, such as the spring snowmelt and heavy summer rains.

390
In spring, mire sub-catchments had the shortest MTTgeo. 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 sub-catchments (like C2) in spring and longer MTTgeo than C2 in winter (Table 5). In C4, the MTTgeo decreased from 1.5 years to 0.7 years from winter to spring, while the corresponding change in C2 was 1.2 years to 0.7 years. The seasonality of MTTgeo was 395 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 MTTgeo that decreased from 7.7 years to 1.9 years from winter to spring (Table 5).

405
The geometric mean of the travel time distribution (MTTgeo) is adjusted for the overland flow. The young water fraction (YWF) includes overland flow and groundwater younger than three months (%). An extended version of the results, including arithmetic mean, median, and SD, is included in the Appendix (Table A1A2).

Testing model results to stream isotopic composition and chemistry
In addition to investigating the annual MTTgeo, 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 base cation (BC) data was available for all sites. In winter, the modelled MTTgeo was correlated to the isotopic 415 composition (r=-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 Δδ 18 Ospring and the young water fraction was also significant (r=-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)

Model results compared to catchment characteristics
The main catchment characteristics correlated to MTTgeo 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 Spring MTT geo areal coverage of silt (r=0.90, P<0.0001) (Fig. 8). A positive correlation between catchment size and MTTgeo 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 to 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).

440
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 MTTgeo and the young water fraction was found for the annual and spring seasonal results but were lost for the summer and winter.  Table 6: Correlation matrixyoung water fraction (YWF), geometric mean travel time (MTTgeo), and catchment characteristics. The catchment characteristics include the log catchment size (log A), the areal coverage of mires (Mire), and the 450 areal coverage of silt (Silt). The table includes annual (grey),, winter (blue),, spring (green),, and summer (orange) results. Darker colours showFor |r| >|> 0.5 with, the connected p-value is shown according to a p<0.05 and b p>0.05. Particle tracking in the Mike SHE model provided valuable insights into the annual and seasonal mean travel times (MTTgeo) across the 14 Krycklan sub-catchments. The modelled MTTgeo and the young water fractions were strongly correlated to observed stream δ 18 Owinter signatures, seasonal variation in δ 18 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 (> four years). Overall, we found that soil type was the most 460 important variable explaining MTTgeo and that mires are an important landscape feature regulating the young water fraction in spring (Fig.8).

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- 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 475 a mix of groundwater and overland flow. Therefore, to allow for actual MTTgeo estimates, we corrected the results by reducing the estimated MTTgeo using the overland flow from the flow model as a scaling factor (Appendix , Table A1A2).
This uncertainty primarily affects the mire dominated sub-catchments that have a large fraction of overland flow, especially during the spring.

480
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 485 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 490 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

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 505 into the streams (Fig. 6). In agreement with our conceptual model (Fig. 2), a strong negative correlation between winter MTTgeo and the isotopic stream signatures during winter baseflow was observed (Fig. 7a). At an average travel time older than four 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 MTTgeo older than four years for some sub-catchments, since their stream isotopic signatures were close to the long-term 510 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 of 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, 515 including overland flow. It is well established that the difference in stream isotopic signature between winter baseflow and spring peak flow at snowmelt (Δδ 18 Ospring) 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 Δδ 18 Ospring 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(Laudon et al., , 2007.

520
Similar to the conditions in spring, the conceptual model predicted that the difference in stream isotopic signature between winter baseflow and summer flow, Δδ 18 Osummer, should be correlated to 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.

525
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 Δδ 18 Ospring, there was still a significant correlation between the average Δδ 18 Osummer and the modelled young water fraction (Fig. 7e).

Controls of travel times on base cation concentrations
The annual and seasonal average BC concentrations were positively correlated to the MTTgeo (Fig. 7b, 7d, and 7f). Since the 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 535 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 540 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, rather enhance it.
Despite arguments that can be made against the use of BCs as tracers, they still offer a complementary possibility to test the model 545 performance (Abbott et al., 2016). As McDonnell and Beven (2014) emphasised, the inclusion of tracers in hydrological models are 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 550 chemistry provides additional support that our modelling of these processesand thus the entire systemwas reasonable and consistent with the empirical data.

Mean travel times, young water fractions, and catchment characteristics
All sub-catchments showed similar seasonal patterns in MTTgeo and young water fraction, manifested as water with shorter travel 555 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 MTTgeo 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 560 size. However, it seems that thisa correlation of travel times to 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 565 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 MTTgeo. In contrast, MTTgeo 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 emphasise that one cannot assume that the travel time would increase with catchment size unless the distribution of different soils is comparable throughout the landscape.

575
A larger catchment size does not necessarily imply that considerably more amounts of old groundwater is 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 MTTgeo is between one 580 and almost eight years. The change in seasonal MTTgeo from winter to spring was also largest for the silt dominated catchments.
For example, there was a six-year difference for C20 compared to the two-year difference for the similar-sized till dominated subcatchment 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 that occurred in spring and summer. The contrasting hydrological response of mire and silt areas, respectively, caused greater annual MTTgeo variation for 585 sub-catchments with both features. For example, the MTTgeo for C4, dominated by mires, decreased from 1.5 years to 0.7 years from winter to spring. In contrast, winter MTTgeo 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 MTTgeo, although C2 (dominated by forest and till) and C4 (dominated by mires) had almost an equal annual MTTgeo (Table 5). In spring, the MTTgeo was shorter in C4 than in C2 due to the surface runoff from the frozen mire.

590
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;. Looking only at the part of runoff originating from groundwater (Appendix , Table A1A2), the MTTgeo for C4 decreases from 1.5 years 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 effects the renewal of the groundwater in mire dominated catchments have long-term consequences on the travel times.

595
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;Kaandrop 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 600 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 that deeper groundwater water transport in till soils in Krycklan is slow. Not far below the groundwater 605 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 (Fig 6, Fig. 8 and Appendix Fig. A1).
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 what 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 615 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 620 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.  (Laudon et al., 2013). Several individuals have also helped with the creation of this work. Special acknowledgement goes to Patrik Vidstrand (SKB) for bedrock properties consultation, Jan-Olof Selroos (SKB) for constructive comments and criticism, Hanna Corell (DHI) for initial particle release consultation, and Anders Lindblom (SKB) with the design of Fig. 1 InThe appendix includes the Appendix,classification of the seasons for each year (Table A2). It also includes the statistical measurements for particle tracking and isotope data are further explained, used in, e.g., tTable 5 and table A1Table A2. The statistics used in this study include arithmetic mean (Eq. A1), geometric mean (Eq. A2), standard deviation (Eq. A3), standard error of the mean (Eq. A4). The isotopic signature of δ 18 O has been calculated as Eq. A5. Eq. A6 shows the equation used for adjusting the MTT from the particle tracking (MTTparticle tracking) to the overland flow MTT shown in table 5. Eq. A7 shows the 920 lake regression from Peralta- Tapia  The start and end dates used in this study for the evaluation of isotopic signatures and BC are shown in Table A1 (Krycklan

930
Database, 2021). 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 to September.   The table includes the travel time results based on particle tracking before taking the overland flow into account. The results