Hydroclimatic non-stationarity drives stream hydrochemistry dynamics through controls on catchment connectivity and water ages

We used data from 17 years of routine stream flow sampling to characterise hydrochemical dynamics in the Girnock Burn, an internationally important long-term catchment monitoring site. We combined hydrochemical time series analysis with hydrological modelling to understand short-and long-term dynamics in relation to dominant runoff processes. Isotopic tracer data employed within the model enabled flux tracking of water as it transited the catchment, facilitating water age estimates. The stream drains an upland 31 km 2 catchment with complex topography and variable geology, which influences soil distributions; peaty soils dominate valley bottoms which receive drainage from more permeable podzolic soils on steeper slopes. Short-term hydrochemical dynamics largely reflect the relative dominance of hydrological sources: summer low flows are dominated by alkaline groundwater high in weathering products (Ca 2 + , Mg 2 + , SiO 2 ) whilst high flows generated by overland flow from peaty soils are more acidic, DOC enriched, and dominated by atmospheric solutes (Na + , Cl - , SO 42). Baseflows are dominated by older groundwater ( > ~3years) in contrast to storm runoff dominated by younger ( < ~0.5 years), near-surface, waters. Usual seasonal distributions of stream flows, with high flows in winter and summer low flows, result in strong seasonality in most solutes. However, coupled effects of hydroclimatic variability and catchment heterogeneity dictates marked scatter in flow-concentration relationships for most determinands. A long-term decreasing trend was evident for SO 42-concentrations, reflecting continuing effects of reductions in atmospheric deposition of pollutants from coal burning. Reductions in large storm event frequency post-2016 and protracted periods of low flows in the droughts of 2018 and 2022 have resulted in the catchment becoming dominated by low flows, derived from older groundwater with more alkaline chemistry for longer periods. Changing climatic conditions and stream flow responses seem likely to continue to cause changes in stream water quality and associated ecosystem services of the Girnock Burn and similar upland catchments.


Introduction
Upland streams and rivers provide crucial, diverse ecosystem services ranging from the provision of drinkable water to supporting habitat for ecologically and economically important species (Ferretto et al., 2021;Bernthal et al., 2022).Multiple factors determine the capacity of these streams to provide such ecosystem services, including water management and local land use; however, a key component is the chemical composition of streams and their source waters (Capell et al., 2011).Complex interactions of multiple catchment characteristics may create strong spatial and temporal variation in the hydrochemistry of upland streams (cf.Jarvie et al., 2005;Cooper et al., 2020).The spatial distribution of a catchment's land use and soil types will affect nearsurface flow paths and solute transport during storm events (Vystavna et al., 2023).On the other hand, the catchment's underlying geology will deliver relatively consistent hydrochemical signals at baseflows derived from more deeply stored water (cf.Shand et al., 2005a,b).Such distinctions will further be reflected in water ages, given baseflows are dominated by older water whereas younger waters are more likely to be evident when near-surface runoff is present during storm events (Stevenson et al., 2022).Importantly, the prominence of any spatial variation in water quality drivers can be directly affected by time-variant landscape-stream connectivity that impacts the relative contribution of contrasting waters to the stream (Soulsby et al., 2015).
In northern temperate or boreal uplands, the time-variable connectivity of peaty, acidic, soils to the stream channel is dictated by antecedent precipitation (P) and subsurface storage dynamics.Variation in the contribution from these shallow source waters provides a contrasting alkalinity signature to well-buffered alkaline groundwaters in contact with the underlying carbonate bedrock (c.f.Smart et al., 2001;Stutter et al., 2012;Oni et al., 2013;Isokangas et al., 2019).Such distinction has enabled hydrograph separation to quantify relative contributions of dominant source waters through the explicit consideration of alkalinity dynamics at catchment outflows (Tetzlaff et al., 2007a).Clearly, such dynamics can be highly variable in the short-term, for example at the storm event scale or with seasonal alterations in precipitation, evapotranspiration and catchment saturation distributions (Soulsby et al., 2001).However, longer-term changes in key hydroclimatic drivers of stream water chemical composition are becoming increasingly evident due to climate change, i.e. the frequency, magnitude, timing and distribution of precipitation events (Kendon et al., 2020).These changes may disrupt the seasonality of catchment stream flow generation, and the temporal distribution and dynamics of the chemical characteristics of stream water (c.f.Larson et al., 2023).
Given the tight links between hydrochemistry and the general health and functioning of a stream ecosystem, climate change is likely to impact the ecosystem services of upland streams (c.f.Orr et al., 2008;Tiwari et al., 2022).Therefore, it is important to be able to assess how current catchment stream flow dynamics, water ages and associated hydrochemical characteristics relate to patterns in the recent past.Such information not only provides an understanding of how current conditions may be changing but also how resilient the integrity of ecosystem services may be in the context of a changing climate.Given the increased evidence and pace of climate change, the potential impacts on hydrochemistry have received greater research focus (e.g.Bussi et al., 2018;Okoniewska and Szumińska, 2020;Pastor et al., 2021;Ponnou-Delaffon et al., 2020;Wu et al., 2021;Liu et al., 2023).
Understanding the controls on upland stream hydrochemistry is often location dependent (geographically and land use type) as well as determinand specific (e.g. the hydrochemical parameter(s) considered).In an upland temperate/boreal context, Shah et al. (2021) considered the River Halladales peaty headwaters in Northern Scotland, observing that two winter hurricanes caused a marked reduction in the acidity and rises in base cations, likely due to climate perturbations causing increased weathering and consequential local bank collapses.Furthermore, Tetzlaff et al., (2007b) also found climatic variability caused marked changes in stream water chemistry in central Scotland.Similarly, Sutton and Price (2022) observed that water quality was sensitive to climatic variability with large temporal changes in Na + concentrations in a mosaic of upland peatlands dominated by fens in Alberta, Canada.Modelling in the UK uplands by Lee et al. (2020) predicted DOC surges in Autumn following drought periods.More widely, Marx et al. (2023) analysed 66 years of stream water quality data from the Panke catchment in Berlin, finding that hydroclimatic changes, particularly in drought years, can affect water quality and alter chemostatic/ chemodynamic catchment solute export patterns.Conversely Bozau et al. (2021) observed that extreme hydrological events in other areas of Germany did not shift hydrochemical parameters outside their seasonal bounds, demonstrating the site and system specificity.
Our goal was to seek further evidence of the influence of climate change on hydrochemistry in upland streams by utilising a long-term hydrochemistry dataset from the Girnock Burn (GB) in the Cairngorms National Park, Scotland.The GB is an internationally important longterm monitoring site which is unusually data rich following a 60-year long programme that has sought to understand the catchment's hydrology and Atlantic Salmon population (Glover et al., 2019;Soulsby et al., 2024).Consequently, the catchment hydrology is well understood, especially following two decades of process-based research that have resulted in one of the most intensively studied sites internationally in terms of hydrological function and water ages (Birkel et al., 2011(Birkel et al., , 2014a(Birkel et al., , 2015;;Soulsby et al., 2015).Furthermore, water quality data have been collected regularly over the past 17 years, which we analyse in their entirety for the first time.Specifically, we explored spatial and temporal patterns of hydrochemical variability in the GB to address the following research objectives: 1.To characterise the long-term temporal variability of the Girnock Burn's stream hydrochemistry.2. To assess the influence of stream flow dynamics and water ages on hydrochemical patterns.3. To use insights from these efforts to identify the longer-and shorterterm drivers of hydrochemical variability.

Study site description
The 31 km 2 Girnock Burn catchment is located in the Cairngorms National Park, Scotland, and is a tributary of the River Dee.The River Dee supports an economically important Atlantic Salmon fishery (Youngson et al., 2002;Scottish Government Marine Directorate, 2023) and provides potable water to > 300,000 people (Stevenson et al., 2022).Detailed accounts of the catchment can be found elsewhere (e.g.Soulsby et al., 2007;Tetzlaff et al., 2007a;Stevenson et al., 2021), but the key information is provided below.
Annual precipitation in the catchment is ~ 1000 mm, where < 10 % of this falls as snow.Annual runoff is ~ 600 mm with average flows of ~ 1.5 mm/d (2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019)(2020)(2021)(2022).Due to the soil types (see below) the GB has a responsive, flashy, hydrological regime with storm runoff dominated by saturation overland flow (SOF) and baseflows sustained from deeper drift deposits.Air temperatures, based on the average annual maximum and minimum of the study period, range from ~>20 • C in summer to ~ -10 • C in winter, when parts of the upper catchment can freeze for prolonged periods (absolute temperature range of the study period being − 15 • C to 29 • C).The catchment ranges in altitude from 230 m to 862 m (average ~ 400 m), with steep hillslopes (up to ~ 61 • ).Generally, it is dominated by heather moorland (~90 %), with areas of conifer forestry and natural pine forest, as well as rough grassland present in the lower catchment which is intermittently grazed by sheep.Natural tree regeneration in the catchment is limited by an overpopulation of Red Deer, Cervus elaphus, with densities approaching 10 deer per km 2 .A small number of residential houses are present in the lower catchment.
The underlying geology is complex with higher elevation areas associated with granite and granite-like rocks (e.g.granodiorite and diorite), with surrounding schists (some calcareous) and other metamorphic bedrock, most notably serpentinite in the East (Fig. 1).Most of the bedrock is covered by glacial and periglacial drift deposits which tend to be poorly draining, particularly in the valley bottoms.Peaty gleys and peaty podzols are the most common soils (~60 % of the catchment).The former tend to dominate in the flatter valley bottom areas and the latter are found on more subdued interfluves and steeper hillslopes.Blanket peat is found in smaller amounts, often adjacent to the stream channels.On the steepest hillslopes in the lower part of the catchment, freely draining humus iron podzols, brown soils, and subalpine podzols are dominant (covering ~ 24 % of the catchment), having developed on a range of bedrock/ drift types.

Hydrochemical, hydrometric, and hydroclimatic data
Between October 2005 and December 2022, the outflow of the GB (Fig. 1) was sampled at least monthly, although more commonly fortnightly, as part of ongoing monitoring by the Freshwater Environment Group of the Scottish Government Marine Directorate.The sampling aimed to characterise water quality variability to improve understanding of Atlantic Salmon productivity within and between catchments of the upper Dee, with only two individual months across the entire time series not sampled.Samples were predominantly taken at low and more average flows (Fig. 2), given such flows dominate, with few samples at the higher end of the hydrograph.Importantly, however, the distribution of stream flows from sampling days generally matched that of stream flow distributions across the entire flow time series (Fig. 2) and thus reflected well the natural distribution of flow variability.Samples were analysed for a wide range of determinands to characterise general water quality parameters as well as major cation, major anion, nutrients, and trace metals (see Table 1 for a full list).Analysis was undertaken at the Marine Directorate Freshwater Fisheries Laboratory (FFL), Pitlochry, Scotland, using methods specifically for low ionic strength waters and following strict protocols for precise hydrochemical monitoring in upland waters (for more information follow references in Gardner, 2008;Monteith et al., 2014).Additional information on the spatial variability of hydrochemistry was provided by a survey of 36 sites on 29th May 2015 covering the main stem of the GB (Fig. 1), ranging from the headwaters to the outflow, with a particular focus on tributary confluences.This high intensity campaign was linked to historical electrofishing sites and aimed to provide a snapshot of spatial variability in water quality.Samples were analysed as per the protocol described above, though with a constrained range of determinands that included Alkalinity, DOC, SiO 2 , Na + , K + , Ca 2+ , Mg 2+ , Cl -, SO 4 2-, Total Soluble N, Total Soluble PO 4 3-and Total Monomeric Aluminium.
Hydroclimatic measurements of precipitation, air temperature, wind speed, relative humidity and net radiation were sourced from within the catchment using an automated weather station recording at 15-minute intervals.Any short period of missing data (e.g.due to sporadic battery or instrument failure) was sourced from the nearby Balmoral weather station, operated by the UK Met Office, and adjusted according to a regression equation (R 2 > 0.9) derived from periods where both GB and Balmoral data were available.Stream flow (Q) was obtained from stage height measurements in a rated section of the channel (via > 150 discrete manual stream flow gauging's in order to create a suitable rating curve) at the catchment outlet.Any missing data was infilled using a rainfall-runoff model previously applied in the catchment and yielding Nash-Sutcliffe Efficiency (NSE) and log NSE performance statistics of > 0.5 (see Birkel et al., 2010).
Finally, P and stream water stable isotopes were collected on a predominantly daily timestep within a sub-catchment of the GB, the 3.2 km 2 Bruntland Burn (BB), between 2007-2009 and 2011 onwards.Previous work over a three-year period showed that the data in the BB were broadly similar to the GB itself (Speed et al., 2011).Samples were collected using ISCO 3700 autosamplers with paraffin added to storage bottles to prevent evaporative fractionation (at 14:00hrs for stream water and cumulatively over a 24-hour period for P).All samples prior to 2018 were analysed for deuterium (δ 2 H) and oxygen-18 (δ 18 O) isotope signatures at the University of Aberdeen using a Los Gatos DLT-100 laser isotope analyser (precision of ± 0.4 ‰ for δ2 H and ± 0.1 ‰ for δ18O).
After 2018, samples were analysed with a same precision isotope analyser at IGB Leibniz-Institute of Freshwater Ecology and Inland Fisheries.Within the BB, from 01/10/2011, we also monitored the average volumetric soil moisture of a 50 cm soil column at the hillslope/ saturation area interface (Fig. 1), which could then be used as a generalised proxy for catchment storage dynamics.

Time series decomposition for trend analysis
To separate longer-term trends from systematic seasonal and finer scale temporal fluctuations we applied a form of time series decomposition to the water quality, hydrological and meteorological (potential driving variables) time series, using locally estimated scatterplot smoothing (LOESS); employing the stlplus function in R (Hafen, 2016).Using this technique, it is possible to set nested smoothing windows at various lengths to extract key temporal components of the data.We trialled multiple combinations of smoothing windows (e.g.smoothing data over various temporal lengths) to characterise the main systematic features of our data, with 2 and 8 year windows found to be optimal in visually elucidating the shorter-and longer-term trends, respectively, of the data from finer scale fluctuations.Moreover, these smoothing periods separated the effect of more hydrologically extreme shorter-term events, associated with inter-annual climate fluctuations, from longer-  term trends.These windows were then fixed across our analyses of different datasets to ensure outcomes were comparable.In short, monthly median data were calculated for each determinand to provide a time series with uniform time steps for the analysis.The (annual) periodicity parameter was set to 12. Underlying assumptions and requirements of the analysis (see Hafen, 2016) required the data be truncated to cover exactly 17 years (i.e. via the removal of the last 3 months of data).Finally, we explored for changes in the stream flow-concentration relationships during the monitoring period by comparing Q and chemical concentration variability over time.We based this analysis on the ratio between the coefficient of variation of each determinand (CV C ) and that of Q (CV Q ); using the formula CV C / CV Q , as per Pohle et al. (2021).We calculated these statistics on two-year moving window using data only from days where coupled Q and determinand samples were available.A threshold value as per Pohle et al. (2021) of > 0.5 was set to indicate chemodynamic behaviour of a determinand (biogeochemically controlled) and ≤ 0.5 indicating chemostatic behaviour (hydrologically controlled).

Water age estimations
The high-resolution isotopic dataset was coupled with a semidistributed, dynamic saturation area flow-tracer model (D-Sat) to provide a first approximation of water ages and thereby improve understanding of hydrological controls on hydrochemistry.Supplementary Material Fig. 1 displays the basic structure of the model, which conceptualizes the catchment as three interacting storage compartments, with hillslope, dynamic saturation area and groundwater reservoirs.Each has an associated dynamic storage, with transfer between these controlled by calibrated linear rate parameters.Stream flow is generated via groundwater using a linear rate parameter, and saturation overland flow, which is controlled by a rate parameter and a non-linearity parameter.Each compartment is also characterized by a calibrated mixing volume parameter to represent the storage needed to dampen the isotope outputs in relation to precipitation inputs, but these volumes do not affect the dynamic water storage fluxes.A key aspect of the model structure is the non-linear expansion and contraction of the saturation area, which is calculated at each timestep using an antecedent precipitation index-type algorithm.Coupling the isotopic data (e.g. the 'fingerprint' of the precipitation and stream flow) with the model structure facilitated the tracking of water as it transited from input (precipitation) to output (stream flow).This time-stamping of inputs and outputs was then able to be used to estimate the age of water egressing the catchment (e.g. the time a particular 'packet' of precipitation spent within the model domain) using the equations detailed in full by Birkel et al. (2015) and Soulsby et al. (2015).For simulations of all fluxes and ages across the 17-year time period, we employed the 500 parameter sets shown (after extensive calibration) by Stevenson et al. (2021) to have the best performance and process representation for the BB catchment (the authors were considering the period 2011 -2018).Model output results were assessed via the Kling-Gupta efficiency statistic (KGE; Kling et al., 2012) for dual, equally weighted, calibration on measured daily stream flow and isotope composition.
Importantly, the D-Sat model was developed specifically for use in the GB (see Birkel et al., 2011;2015) and has proved effective in modelling the similar runoff dynamics, isotopic composition, dominant hydrological sources (e.g.groundwater or saturation overland flow) and water ages of both the GB and BB and elucidating internal catchment functioning and the consequences to stream flow volume and water ages.For instance, Birkel et al., (2011 demonstrated how runoff generation from the high water storage saturation area of the GB and BB occurs predominantly during wetter periods due to enhanced stream network connectivity, resulting in a reduction in water age estimates.The authors also found that the systems water age non-stationarity was likely to be sensitive to more extreme hydroclimatic variability.Furthermore, modelling by Birkel et al., (2014a) and Soulsby et al. (2015) revealed the importance of antecedent conditions to catchment connectivity, with drier period stream flows being sustained from older groundwaters as SOF reduce or entirely cease; causing a disconnect between the saturation area and stream channel.Moreover, when conditions cause the saturation area to reconnect, the age of the consequential SOF may be elevated in comparison to SOF occurring after longer-term wet conditions given the contributing volumes time within the (previously disconnected) saturation area store.
Given the extensive isotope dataset available for the BB, and therefore the process-based insights on catchment functioning across the study period, the model was applied solely to this sub-catchment.Importantly however, runoff generation and stream water age distributions of both the GB and BB have been shown to be very similar (see Birkel et al., 2015;Soulsby et al., 2015), allowing derived age simulations from the BB to be used as a first approximation to assess against hydrochemical dynamics in the GB.Furthermore, given well known complexities and time variable uncertainties in determining precise water ages (especially for older waters), we normalised the generated water age data between 0 and 1, suitable to provide transferable dynamic patterns and time series to GB hydrochemical data analysis.We also grouped the water age data from sampling days into four percentile groupings (0-5 %, >5%-33 %, >33-66 %, >66 %).This provided broad age categories to assess against determinand concentrations, whilst inclusion of the constrained lowest category was designed to ensure that the information content from the relatively few higher flows sampled (and when it is likely the water will be youngest) was independently grouped and captured.

Hydroclimatic time series trends
The study period (2005 -2022) encompassed a wide range of climatic conditions, with some extreme events, especially in relation to the magnitude and distribution of P and air temperature (Fig. 3).The first two years of the study period were typical of the longer-term average conditions with wet winters, characterised by medium sized, evenly distributed, rainfall and runoff events, with drier summers.However, the summer of 2007 was notably wet, whilst the winters of 2009/10 and 2010/11 were cooler than average, with conditions after this having more typical seasonality until a warmer, wetter, winter in 2013/14 (Fig. 3).The winter of 2015/16 was notable for a large storm (Frank) where heaviest rain occurred on 30th December but continued to bring wet weather in the first week of January.The rainfall peak in December saw the most extensive flood event in the Dee catchment for over 200 years (Soulsby et al., 2017).Subsequently, conditions turned generally drier than average (notably in the winters) and variable, especially following a drought in 2017/18, when negative precipitation anomalies became more common; most severely during the summers of 2018 and 2022.These years were, however, interspersed by a wetter 2020/21 (Fig. 3), whilst the very end of the time period saw a series of relatively large events occurring in autumn 2022.
Q largely mirrored P dynamics, especially during the wetter periods of higher catchment connectivity.During the more extreme dry periods (e.g.summers of 2018 and 2022) baseflows, sustained by groundwater, dominated for extended periods (Fig. 3).Air temperature exhibited relatively consistent inter-annual seasonality (Fig. 3), though in 2009, 2010 and post 2018 winters were notably cold (to − 15 • C) which caused parts of the catchment to freeze.Contrastingly, summer highs remained fairly consistent across the time period.

Hydrochemical time series trends
The stream was typically mildly alkaline with a mean pH of 7.19 (though notably mildly acidic during higher, near-surface dominated, flows), whilst also being of low ionic strength and with relatively high DOC concentrations (Table 1).Nutrient levels were generally low, though intermittent enrichment was evident given the maximum levels of PO 4 3-and Total Soluble Organic N (Table 1).More generally, the relative contributions of major ions to the hydrochemistry of the stream water was highly variable, as evident in a Piper plot (Fig. 4).Many major ions showed clear inter-and intra-annual dynamics (Figs. 5 and 6, respectively), which can be related to the influence of antecedent wetness (as inferred from soil water content; Supplementary Material Fig. 2) on catchment functioning.For instance, during drier periods and lower flows (approximately May -Oct), stream water showed higher concentrations of groundwater influenced solutes such as alkalinity, Ca 2+ , Na + and Mg 2+ (Fig. 6).Peaks in these solutes were also observed during, and preceding, drought periods when catchment storage was being reduced and consequently reached time series lows (Supplementary Material Fig. 2).Contrastingly, these hydrochemical trends generally reversed during higher flow winter periods when SOF increased.Across the time series as a whole, and especially post 2018, SO 4 2decreased notably, as did Cl -, though to a lesser extent (Fig. 5).In contrast, SiO 2 concentrations increased, most markedly in the early part of the record pre-2011, whilst DOC showed a general trend towards lower summer concentrations, most notably between 2011 and 2018, though higher peak values were also present as the time series progressed post 2017 (Fig. 5).
The drought of 2018 had a marked short-term impact on a range of hydrochemical characteristics (Fig. 5).For instance, Alkalinity, Na + , K + , Mg 2+ , Cl -and SO 4 2-all experienced substantial increases ~ 2018 in context of their monthly median average and shorter-term trends.Conversely, after this period concentrations declined to their lowest values of the time series, whilst DOC, Total Soluble N and Total Soluble PO 4 3-experienced relative peaks, with similar patterns evident during and after the next drought period in summer 2022 (Fig. 5).

Flow -Concentration relationships
Most chemical determinands exhibited relationships with stream flow, with all but four exhibiting strong linear relationships (threshold of p < 0.05; Table 3 and Fig. 7).However, the majority of these relationships also exhibited considerable scatter (e.g.Figs. 4 and 7).This was likely related to the hydroclimatic variability experienced throughout the study period, particularly in the last few years, and the impact on antecedent conditions and more extreme catchment storage dynamics (Supplementary Material Fig. 2).Overall, alkalinity, conductivity, Ca 2+ , Mg 2+ and pH had the strongest relationships with Q (R 2 > 0.50), given the first four of these are predominantly associated with groundwater, which dominates low flows.Conversely, pH is most influenced by overland flows which have come into contact with acidic soils, these being most prevalent during higher Q periods.These soil contacting SOF flows also caused DOC to evidence a relatively strong relationship with Q (Fig. 7).
More moderate, though inverse, relationships were present for SiO and Total Monomeric Aluminium, the former evidencing flow dilution.Conversely, all three PO determinands exhibited no clear relationships with stream flow (Table 3).During the very highest sampled flows, concentrations of all determinands were diluted.Visualisation of monthly flow-concentration relationships (selective solutes shown in Fig. 7) further illustrated seasonal effects, with samples from June and July frequently grouping together (though less so for Ca 2+ ), as did the winter months of November and December.Conversely, samples from January/ February and especially August showed considerable scatter, coinciding with periods when the catchment would more commonly be partially frozen or relatively dry and potentially subject to convective rainfall, respectively (Fig. 3).
Eight of the 23 determinands demonstrated CV C /CV Q ratios of less than 0.5 throughout the time series, with these being predominantly weathering derived solutes (e.g.Ca 2+ , Mg 2+ and Na + ) or those heavily influenced by the peaty soils (e.g.pH).During the two-year windows between 2013 and 2017 all determinands were chemostatic, coinciding with a period of relatively increased P (Fig. 3).Conversely, apart from three solutes (SiO 2, NH 4 + and PO 4 3-), all chemical determinands exhibited their greatest CV C / CV Q ratio value after the 2018 drought.Notably, post 2018 DOC switched from being < 0.5 to consistently being > 0.5 (Table 2).

. Model application and simulated dynamics
Overall, the tracer-aided model was able to capture the inter-and intra-annual variability in both Q and isotope data reasonably well with KGE's of 0.64 and 0.77 for median simulation values, respectively; demonstrating the goodness of fit over the time period.Moreover, confidence bounds were relatively well-constrained (Fig. 8), even during the drought of 2018 when catchment functioning was substantially altered.However, more extreme isotopic signals and flows derived from more unusual, often individual, high flow events were captured less well (Fig. 8).Nonetheless, the model performance gave confidence that the derived water age dynamics were representative (Fig. 8).The simulated SOF dynamics from model outputs (Fig. 8) provided a proxy for connectivity between the landscape and the stream channel, reflecting the importance of this component water source in the catchment outflow.This is particularly so under higher flows when SOF was simulated as contributing considerable volumes to stream flow (see Section 2 and Supplementary Material Fig. 3).Consequently, during late autumn and winter months when higher flows were more commonplace, SOF flow was also more dominant (Fig. 8 and Supplementary Material Fig. 3).
Water was generally younger during higher flows (averaging a few months) and periods when such flows were most common (e.g.winter; Fig. 8 and Supplementary Material Fig. 3), whilst the inverse was true in summer when groundwater contributions proportionally increased (to ages > 3 years).However, there were also instances when higher flows were simulated as being comprised of relatively older water  (Supplementary Material Fig. 3), representing periods when the saturation area store was reconnected to the stream via SOF following drought (and associated aging of water within the saturation store).There was also a clear trend in the smoothed data of simulated Q ages increasing from the start of the time series to ~ 2016/17, when consecutive drier winters substantially elevated water ages until the 2018 summer drought ended.This was followed by a reduction in water ages as the catchment rewetted, though water remained older in respect to ages before the drought.Moreover, the highest median (calendar) monthly age estimate was experienced in 2020, with peaks also observed during the summer drought of 2022.Following this a series of relatively large autumn storms caused a sudden reduction in water ages, to the lowest of the entire study period.Collectively, these events demonstrated more extreme variability post 2018 in water age dynamics.
Simulated SOF closely mirrored Q variability (Fig. 8 and Supplementary Material Fig. 3) remaining relatively consistent until (Fig. 8), after which volumetric contributions exhibited a downward trend reaching their minimum during the drought of summer 2018.During this period, the saturation area in the model was essentially disconnected from the stream for ~ 2 months.Contributions increased again during 2020/2021, though remained below the timeseries previous norms given relatively low P (Fig. 3).During the summer of SOF simulations once again ceased before increasing again in the autumn of 2022 following relatively large and repeated P events (Fig. 3).

Water age -Concentration relationships
Most solutes demonstrated a statistically significant relationship with water age (P < 0.05, Table 3), whilst grouping water ages into four categories and comparing to select determinands also demonstrated these trends (Fig. 9).However, as with Q -Concentration analysis, there was considerable scatter in these relationships.The highest coefficients of determination (e.g.those exhibiting the strongest statistical relationships) were for predominantly groundwater and weathering derived determinands such as alkalinity, Na + , Ca 2+ and Mg 2+ .Conversely, Total Soluble PO 4 3-showed little relationship with water age categories, whilst there was considerable variability in K + concentrations, with the 6 to 33 % age bracket having outliers spanning the full range of sample concentration (Fig. 9).Given the relationships between Q and determinand concentrations, it was unsurprising similar relationships between water ages and determinands were also observed (Table 3).However, the results were not directly comparable, with Absorbance(s), DOC and Total Soluble N (Table 3) exhibiting much stronger relationships with water age than Q, reflecting potential variation in source water contributions under similar stream flows, as discussed in Section 5.1.

Spatial characteristics of stream hydrochemistry
The spatial sampling undertaken in May 2015 (Section 3.1) highlighted the influence of the catchment's geology and soils (Fig. 1) on hydrochemistry (Fig. 10).The upper part of the catchment was generally characterised by lowest concentrations of most solutes, reflecting the predominance of granite in the underlying geology.Notably, however, Total Soluble N and PO 4 3-and Total Monomeric Al were higher in the upper catchment.The more weatherable metamorphic geology, and its effect on drift deposits and soil mineralogy, influenced hydrochemistry in the mid-catchment, with increasing concentrations of SiO 2 , Na + , K + , Ca 2+ and particularly Mg 2+ which increased in response to the Serpentinite outcrop in the east (Figs. 1 and 10).The influence of soil characteristics was reflected throughout the catchment, for example DOC concentrations were greatest just downstream of the mid-section where peaty soils dominate (Fig. 1).Conversely, the higher easterly tributaries contributed relatively minimal DOC concentrations due to differences in soil types (Section 2 and Fig. 1) and reduced soil depths on steeper slopes.Finally, at the catchment outflow, there was a notably

Hydrochemistry and the impact of flow, water age and spatial characteristics
Monitoring over 17 years captured a wide range of climatic conditions and revealed a complex and changing hydrochemical regime in the GB, with recent reductions in P similar to those being observed and predicted elsewhere (Konopala et al., 2020;Smith et al., 2021).In general, the stream was mildly alkaline with low ionic strength, high DOC and low nutrient concentrations (Table 1).However, the hydrochemical composition was also highly temporally variable (Figs. 5 & 6), influenced by the relative contribution of different source waters.SOF was derived from peaty, acidic, DOC rich, soils which caused a hydrochemical signal dominated by ions associated with stronger acids.This contrasted with the groundwaters which sustained lower baseflows and were typified by greater alkalinity and weaker acids given the prevalence of weathering products; characteristics typical of boreal regions (Burd et al., 2018;Scheliga et al., 2019).For most of the monitoring period, i.e. up to ~ 2016, this distinction in water sources underpinned relatively consistent seasonal cycles in predominantly chemostatic determinands, as has been observed elsewhere (c.f.Minaudo et al., 2019).DOC, alkalinity, Ca 2+ and Mg 2+ varied strongly with P as groundwater dominated drier summer months and SOF had greater influence during the wetter winter months when catchment connectivity was higher.Post 2016, P inputs became more variable and extreme (Fig. 3), altering the spatial extent of the catchment saturation area, with consequences for rate and timing of SOF thereby influencing hydrochemical signatures (discussed below).
Given the importance of water source area to stream hydrochemistry   (Tetzlaff et al., 2007b;Jarvie et al., 2008), Q was a statistically significant predictor of stream concentration for the majority of determinands (Table 3 and Figs. 4 and 7) and a stronger predictor for weathering derived, groundwater influenced solutes.However, scatter in these relationships also reflected the complexity of stream flow generation, including potential hysteresis effects masked by the sampling pattern.
Complexity was also more likely due to the extensive monitoring period and substantial variability in antecedent conditions prior to sampling, both in terms of precipitation dynamics and catchment storage (Fig. 3 and Supplementary Material Fig. 2).For example, during drier periods small P events may have connected individual saturation area pools to the stream leading to short-term increases in DOC concentration that were disproportionally high relative to the increase in Qthough further high temporal/ spatial resolution data would be required to corroborate this hypothesis.Furthermore, the fact that water age was often a better predictor of determinand concentration than stream flow reflected potential variation in source water contributions under similar stream flows.Such variation would result from changes in antecedent conditions, storage dynamics (Supplementary Material Fig. 2) and SOF; factors directly effecting both determinand concentrations and water ages.These findings demonstrated the value of considering water age estimations in the analysis of stream water quality and highlight how process-based understanding derived from modelling helped elucidate the effects of internal catchment functioning.Furthermore, water ages proved informative in supporting linkages between the predominant source waters and determinand concentration, with older waters dominated by groundwater derived solutes, as observed by others (c.f.Tunaley et al., 2016;Jepson et al., 2019).Moreover, changes in the hydrochemical composition of the stream in the later monitoring period when large storms followed drought periods were clearly reflected in modelled water ages (see Section 5.2).The importance of atmospheric deposition and pollutant control measures on stream hydrochemistry was also evident, with an overall decline in SO 4 2-in response to declining sulphur emissions from the burning of coal, a key driver of acid rain (Monteith et al., 2014;Lu et al., 2021).Reductions in SO 4 2-concentrations decelerated towards the end of the time series, with the remaining concentrations reflecting previous soil retention of SO 4 2-and its continued gradual release (Palmor et al., 2013).Cl -also saw a notable decrease, most clearly post 2018, which possibly reflected reductions in storm prevalence and lower inputs of marine-derived solutes, as well as reductions in coal burning (Evans et al., 2011).DOC showed a general trend of reduced summer concentrations towards the end of the timeseries, indicating reduced connectivity between the peat soils and stream during these months.Conversely, over the monitoring period SiO 2 concentrations increased substantially, despite being already elevated due to weathering influences (Fig. 5).The drivers of this behaviour are complex, though may in part reflect wider biogeochemical changes over time which would, for example, impact diatom uptake of SiO 2 and water column concentrations (see Conley et al., 1993).
The spatial variability of stream water composition in the catchment had a clear impact on the determinand distribution and concentration, as has been observed elsewhere (c.f.Shand et al., 2005a,b;Kiewiet et al., 2019;Lupon et al., 2023).In the middle of the catchment, increasing concentrations of weathering derived solutes were observed in response to increasing influence of metamorphic geologies, in contrast to the granites that dominate the headwaters, the latter likely having lower weathering rates in this environment (though this may be less so elsewhere, c.f. Hayes et al., 2020) and less extensive fracture networks for groundwater flow (Soulsby et al., 1998).Moreover, the hydrochemical importance of geological constituents was exemplified by Mg 2+ , derived predominantly from a Serpentinite outcrop, which contributed solutes to the main stem via the mid-easterly tributary (Fig. 10).The prevalence of peaty soils and their proximity to the main stem also increased other determinands such as DOC in central catchment (Birkel et al., 2014b;Dick et al., 2015).A local source of Total Soluble N near the outflow was likely derived from sheep grazing on the lower catchment's rough grassland and drainage from forested areas contiguous to the stream.Muller et al., 2015).
Finally, it should be noted that the routine sampling occurred predominantly during lower and average flows, given such flows dominate the hydrograph.As a result, the higher end of the hydrograph and the relationship with determinands are not fully characterised in this dataset, which could have promoted stronger statistical relationships between groundwater-associated determinands and flow, given lower flows will coincide with higher groundwater contributions.However, as shown in Fig. 2, the distribution of hydrochemical sampling generally covered flow distributions more widely.Therefore, sampling encompassed well the general variability of stream flows, including time periods when SOF would have been present, diluting groundwater contributions (Fig. 2, Fig. 8 and Supplementary Material Fig. 3).Moreover, those higher flows not captured by hydrochemical sampling were relatively rare (Fig. 2), meaning even if inclusion within statistical analysis had been possible their influence is likely to have been relatively minor.

Impact of negative rainfall anomalies
The summer drought of 2018 in Northern Europe (Kleine et al., 2020) was especially acute given it followed two anomalously dry, consecutive, winters which caused depleted soil moisture stores within the catchment (Soulsby et al., 2021).Consequently, stream flow was dramatically reduced during the summer of 2018 and derived from groundwater fed baseflows for extended periods (Figs. 3 and 8).This was reflected in our modelling, which simulated the cessation of SOF during these periods (Fig. 8), implying that the peaty saturation area effectively became hydrologically disconnected from the stream, as was observed in field observations (Soulsby et al., 2021).Furthermore, modelled water ages increased rapidly during the drought (Fig. 8).
The hydrochemical impact of such changes with respect to historic norms was substantive, with the majority of determinands showing large changes in concentration dynamics, most notably during the summer of 2018 (Fig. 5), with ratios in variability generally dropping prior to this period and spiking after (Table 2).Given the proportional predominance of groundwater at these times, weathering derived determinands were elevated and recorded their highest concentrations (e.g.alkalinity, Mg 2+ , Ca 2+ etc.), with a corresponding reduction in those derived from saturated peaty soils (e.g.DOC; Fig. 5).Following the 2018 drought, there was a partial, relatively short-lived, move back to towards longerterm hydroclimatic patterns (as observed elsewhere in Europe, c.f. Smith et al., 2020).However, the consequential catchment storage deficit (Supplementary Material Fig. 2) and insufficient P inputs (Fig. 3) meant runoff responses were substantially damped, even during the winter.This shift caused further anomalies in determinand concentrations, with time series peaks of DOC, Total Soluble N and PO 4 3-; dynamics inconsistent with prevailing flowdeterminand relationships.Furthermore, the moving window CV C / CV Q ratio showed flowconcentration relationships changes for determinands after 2018, with 19 of the 23 determinands characterised by their highest CV C / CV Q value during this time (Table 2).Though this does not mean that all solutes became chemodynamic, it did hint towards a substantive change in catchment functioning and a more complex relationship between determinand concentration, flow rates and water sources during the post-drought rewetting period.Such dynamics demonstrate how catchments depleted in long-term storage after long droughts are less likely to show time-invariant outputs in comparison to extensive well mixed reservoirs present in more normal hydrological periods (Soulsby et al., 2021).For example, drying of the saturation area during 2018 would facilitate increased carbon decomposition rates (in ecosystems valued for their low decomposition rates and carbon storage) thereby flushing greater concentrations of DOC when connection to the stream was reestablished (c.f.Tunaley et al., 2017), even with relatively small increases in stream flow.Such findings are consistent with other studies where changes in climate have been hypothesised to be driving changes in stream hydrochemical dynamics (e.g.Ponnou-Delaffon et al., 2020), though the extent of this is dependent on the sensitivity of the landscape in question (c.f.Ulloa-Cedamanos et al., 2022).
The summer of 2022 was also anomalously dry (Fig. 3) and characterised by reductions in DOC and relative increases in weathering derived products (Fig. 5), similar to those seen in 2018.However, in contrast to 2018 a series of relatively large P events in the early autumn caused a rapid rewetting of the catchment.Consequently, intracatchment connectivity was extensively and quickly re-established, with SOF spiking (Fig. 8), along with determinands associated with short residence water sources (Fig. 5).Taken together the anomalous and extreme effects of the 2018 drought and 2022 summer low flows clearly demonstrate how climate change induced negative rainfall anomalies and extremes can exacerbate the chemostatic and hydrochemical dynamics of the system (c.f.Tiwari et al., 2022;Zhang et al., 2022).

Uncertainties in ecosystem service provision in a changing hydroclimatic environment
Catchments such as the GB provide multiple ecosystem servicesfor instance, flood protection, nutrient cycling and support of the wider riverine ecosystemthough most importantly the stream is a tributary of the Dee which provides potable water to > 300,000 people and is economically important for Atlantic Salmon fishing (Soulsby et al., 2007;Marine Scotland, 2017;Keele et al., 2019;Stevenson et al., 2022).Therefore, changes in catchment functioning that will affect these services need to be understood so that their impact can, where possible, be managed and mitigated (Curtis et al., 2014).The data presented here have shown that changes in stream flow, water ages and water quality are principally driven by hydroclimatic variability (Fig. 3) and catchment storage (Supplementary Material Fig. 2).Therefore, mitigation of any negative hydrochemical change would need to be mediated indirectly, e.g. through changes in catchment management.However, land use changes would also need to be carefully considered to avoid synergistic, negative, impacts to stream hydrochemistry in the context of a changing hydroclimatic environment.
For example, afforestation is relatively commonplace in Scotland due to its economic importance (Scottish Government, 2023) and potential for flood alleviation (Black et al., 2021), though the efficacy of this in large floods is questionable (Soulsby et al., 2021).Importantly, increased afforestation of the uplands is expected to elevate catchment water losses via higher transpiration and interception evaporation, thereby reducing soil storage, groundwater recharge, catchment connectivity and potentially summer low flows (Neill et al., 2021).In periods of negative rainfall anomalies it is, therefore, likely that catchment scale re-forestation could exacerbate summer low flows and thus affect the hydrochemical composition of streams more extensively and for a longer period of time.Such changes would likely propagate to impacts on ecosystem services such as water extraction and treatment (e.g.removal of DOC colouration), which are already being pressured by current changes in climate.However, our understanding of such mechanisms is predicated on historical patterns, with future climate changes likely to become more variable, and so may be the impacts.Consequently, future work needs to both continue monitoring hydroclimatic and hydrochemical trends, especially in longer-term experimental catchments and those associated with declines in key ecosystem service provisions (Tetzlaff et al., 2017).

Conclusion
Over a 17-year period (2005 -2022), we monitored the hydrochemistry of the GB, an internationally important monitoring site with a rich, diverse and extensive hydrological dataset, including stable water isotope data which facilitated the application of a tracer-aided model to derive water age estimates and understand internal catchment functioning.Using these data we characterised the long-term hydrochemistry of the GB and assessed its relation with stream flow, water ages, catchment connectivity and spatial heterogeneity in geology, soils and land use.
Our results showed that short-term hydrochemical dynamics largely reflect the relative dominance of hydrological sources, with summer low flows dominated by alkaline groundwater high in weathering products (Ca 2+ , Mg 2+ , SiO 2 etc.) whilst high flows generated by SOF from peaty soils were acidic, enriched in DOC and dominated by atmospherically derived solutes (Na + , Cl -and SO 4

2-
). Meaningful relationships were demonstrated between most of the 23 determinands and stream flow or water age, given lower stream flows are dominated by baseflows from older groundwaters.However, there was considerable scatter in the majority of these relationships, related to the hydroclimatic nonstationarity experienced throughout the study, which impacted, for example, catchment antecedent storage (Supplementary Material Fig. 2) and connectivity dynamics.
The evident non-stationarity in catchment functioning became increasingly variable and extreme towards the end of the monitoring period due to climate change induced alterations in P rates and timings, in turn impacting hydrochemical dynamics of the stream.For example, extended periods of groundwater sustained summer low flows occurring during the droughts of 2018 and 2022 as the peaty saturation area became hydrologically disconnected from the stream, with associated increases in weathering derived solutes.Conversely, post drought peaks in determinands such as DOC occurred, though these did not match historic relationships with flow rates, possibly reflecting increased decomposition rates during droughts and relatively rapid catchment reconnectivity and consequent flushing.
Such changes in catchment functioning and the resultant impact on hydrochemical dynamics is likely to impact key ecosystem services such as water extraction, water quality and treatment, especially if coupled with other factors such as land use changes.Future work should continue to monitor hydroclimatic and hydrochemical dynamics in longer-term experimental catchments such as the GB, and especially

Fig. 1 .
Fig. 1.General map of (A) bedrock geology and (B) soil types of the Girnock Burn, indicating flow monitoring, routine sampling location at the out flow, synoptic sampling locations and location of the Bruntland Burn sub-catchment.(C) shows detailed soil map and sampling locations of the Bruntland Burn.Top left map indicates geographical reference to United Kingdom.

Fig. 2 .
Fig. 2. Histograms of stream flow volumes across entire daily timeseries (bottom panel) and stream flow volumes on only days when a hydrochemical sample was taken (top panel).Note differences of y axis scale.

Fig. 3 .
Fig. 3. Daily (left column) and smoothed, using seasonal time series decomposition Locally Estimated Scatterplot Smoothing regression analysis (right column) time series for precipitation (top row), stream flow (middle row) and air temperature (bottom row).Red vertical shading indicates general timings of 2018 and summer droughts.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 4 .
Fig. 4. Piper plot of all stream flow determinand samples, coloured according to flow percentile.

Fig. 5 .
Fig. 5. Smoothed timeseries of selected determinands using seasonal time series decomposition Locally Estimated Scatterplot Smoothing (LOESS) regression analysis.Monthly values in grey indicate median observed value of calendar month.Red vertical shading indicates general timings of 2018 and 2022 summer droughts.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 6 .
Fig. 6.Distribution of monthly mean values, calculated from the full timeseries, for select determinands.Outer box indicates first to third quartile of the distribution, inner line the median and outer vertical lines the variability outside these first and third quartiles.Circles represent outlying values.

Fig. 7 .
Fig. 7. Selected determinands plotted against the natural log of stream flow.Colours indicate the month of the year, with the statistics derived from linear regression analysis across the entire dataset.
.Stevenson et al.   strong pulse of Total Soluble N.

Fig. 8 .
Fig. 8. A) Timeseries of monitored and simulated deuterium at the catchment outflow, B) stream flow, C) saturation overland flow and D) normalised stream flow water age.In panels A and B black points represent observed values, blue lines median simulations and orange shading 05th to 95th percentile confidence intervals.In panels C and D, red line represents seasonal time series decomposition Locally Estimated Scatterplot Smoothing (LOESS) regression analysis employing a 2 year window of smoothing.Grey lines indicate monthly sum and median for panels C and D, respectively.Red vertical shading indicates general timings of 2018 and 2022 summer droughts.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 9 .
Fig. 9. Box plotting of stream flow water age brackets against selected determinands.For reference, the youngest water bin was <~0.5 years, with the oldest >~3 years.Outer box indicates first to third quartile of the distribution, inner line the median and outer vertical lines the variability outside these first and third quartiles.Circles represent outlying values.

Fig. 10 .
Fig. 10.Determinand concentration maps from the high spatial resolution sampling undertaken on 29th May 2015.

Table 3
Determinand relationships with natural logarithm of stream flow rate and water age.Calculated via linear regression analysis using all continuous data.