Groundwater flow paths drive longitudinal patterns of stream dissolved organic carbon (DOC) concentrations in boreal landscapes

Preferential groundwater flow paths can influence dissolved organic carbon (DOC) concentration and export in the fluvial network because they facilitate the inflow of terrestrial DOC from large upslope contributing areas to discrete sections of the stream, referred to as discrete riparian inflow points (DRIPs). However, the mechanisms by which DRIPs influence longitudinal patterns of stream DOC concentrations are still poorly understood. In this study, we ask how DRIPs affect longitudinal patterns of stream DOC concentrations under different hydrologic conditions, as they can simultaneously act as major sources of terrestrial DOC and important locations for in-stream processes. To answer this question, we tested four model structures that account for different representations of hydrology (distributed inflows of DRIPs vs. diffuse groundwater inflow) and in-stream processes (no DOC uptake vs. in-stream DOC uptake downstream of DRIPs) to simulate stream DOC concentrations along a 1.5 km headwater reach for 14 sampling campaigns with flow conditions ranging from droughts to floods. Despite the magnitude and longitudinal patterns of stream DOC concentration varying across campaigns, at least one model structure was able to capture longitudinal trends during each campaign. Specifically, our results showed that during snowmelt periods or high-flow conditions (>50 L s⁻¹), accounting for distributed inputs of DRIPs improved simulations of stream DOC concentrations along the reach, because groundwater inputs from DRIPs diluted the DOC in transport. Moreover, accounting for in-stream DOC uptake immediately downstream of DRIPs improved simulations during five sampling campaigns that were performed during spring and summer, indicating that these locations served as a resource of DOC for aquatic biota. These results show that the role of DRIPs in modulating DOC concentration, cycling, and export varies over time and depends strongly on catchment hydrology. Therefore, accounting for DRIPs can improve stream biogeochemistry frameworks and help inform management of riparian areas under current and future climatic conditions.

1 Introduction

Streams and rivers play a critical role in the global carbon (C) cycle because they transport, store and process large amounts of dissolved organic carbon (DOC) (Ciais et al., 2013). Accounting for spatial patterns of stream DOC concentrations within stream networks is vital for understanding net in-stream C retention along rivers (Alexander et al., 2007; Bernal et al., 2018) and catchment-integrated evasion of C (Wallin et al., 2013) as well as for assessing and managing the brownification of large water bodies and coastal ecosystems (Kritzberg et al., 2020). Yet the main drivers controlling spatial variations in DOC concentrations remain unclear, partly because processes occurring at various scales interact in complex ways to influence the concentration and export of DOC to downstream aquatic ecosystems (Laudon and Sponseller, 2018).

In boreal regions, landscape features such as wetlands, headwater lakes, and riparian zones are major controls of the spatial variability in stream DOC concentrations (Frost et al., 2006; Laudon et al., 2011; Lottig et al., 2013; Kothawala et al., 2015). In peat-rich riparian soils, typical for boreal forest catchments, the combination of wet soil conditions and organic matter accumulation can result in elevated DOC concentrations in subsurface water (Grabs et al., 2012). The organization of groundwater flow paths can also regulate spatial patterns of stream DOC concentration by conveying substantial fluxes of water from large upslope contributing areas through wet corridors to discrete sections of the stream (Jencso et al., 2010; McGlynn and McDonnell, 2003), referred to as discrete riparian inflow points (DRIPs; Ploum et al., 2018). DRIPs have been shown to have high groundwater concentrations of DOC associated with sustained water-saturated conditions, moss-dominated vegetation, and organic matter accumulation (Demars et al., 2020; Ploum et al., 2021). Further, the strong connectivity between DRIPs and adjacent streams makes DRIPs important sources of DOC that can be adsorbed, photodegraded, or mineralized by aquatic microbial communities (Berggren et al., 2009; Mineau et al., 2016). This processing happens quickly and over relatively short distances (Demars, 2019), thereby generating hotspots of in-stream uptake immediately downstream of DRIPs (Lupon et al., 2020). To integrate the role of DRIPs as both suppliers of terrestrial DOC and hotspots of DOC uptake in aquatic ecosystems, we need to combine source-transport hydrochemical and in-stream C cycling frameworks (Li et al., 2020). While there are hydrological frameworks that account for flow path convergence (Jencso et al., 2009; Seibert et al., 2009), these features are often not explicitly considered in biogeochemical studies and monitoring strategies (Briggs and Hare, 2018). As a result, the extent to which DRIPs can affect stream DOC concentrations, cycling, and overall C exports at different spatial scales remains largely unknown.

The relative contribution of DRIPs to shaping downstream DOC concentrations and processing likely depends on catchment hydrology. During baseflow conditions and small rain events, DRIPs are major contributors of water to streamflow (Leach et al., 2017; Ploum et al., 2018), and hence they could drive spatial variation in stream DOC concentrations by acting as sources of DOC along streams, as observed for C gases (Duvert et al., 2018; Lupon et al., 2019). In contrast, the relevance of these flow paths as primary drivers of stream DOC concentrations might be less important during extreme hydrological events (i.e., droughts and floods), when the DRIP-stream hydrological connectivity is low or overwhelmed by either upstream fluxes or diffuse lateral inflows (Leach et al., 2017; Gómez-Gener et al., 2020). Further, catchment hydrology also affects the potential for aquatic biota to act upon the DOC in transport (pulse-shunt concept; Raymond et al., 2016). Large residence times during low and moderate flows can promote in-stream DOC mineralization (Casas-Ruiz et al., 2017), while elevated water velocities might overwhelm in-stream DOC uptake during high flows (Bernal et al., 2019). Since the hydrology of many boreal landscapes is rapidly changing due to global change (Laudon et al., 2021), it is important to understand where and when DRIPs hydrologically connect to headwaters as well as their broader effects on stream C cycling and overall catchment C export under current and future climatic scenarios.

In this study, we assessed the relevance of DRIPs as primary drivers of spatial patterns of stream DOC concentrations along boreal headwater streams. Specifically, we aimed to disentangle the role of DRIPs as terrestrial DOC suppliers vs. hotspots for in-stream DOC uptake during different flow conditions. To do so, we tested four different models to simulate stream DOC concentrations along a 1.5 km headwater reach for 14 campaigns with flow conditions ranging from droughts to floods. Models accounted for two types of transport mechanisms. We assumed either (i) uniform, diffuse inflow of groundwater along the reach or (ii) the existence of DRIPs by weighting groundwater inflow relative to their upslope contributing area (UCA). These two assumptions about groundwater inflow were combined with the assumption that stream biota do not take up the supplied DOC (i.e., pulse-shunt concept) or that in-stream DOC uptake takes place directly downstream of DRIPs (i.e., hotspot concept; Table 1).

Table 1 Overview of the different model assumptions. The first column indicates the model name. The second column indicates whether streamflow is represented as a uniform diffuse rate along the reach or distributed based on upslope contributing area. The third column indicates whether in-stream uptake of dissolved organic carbon (DOC) by biota is included.

2 Methods

2.1 Study area

We conducted our study in the Krycklan catchment in northern Sweden (64^∘14^′ N, 19^∘46^′ E) along a 1.5 km stream reach located between the gauging stations C5 and C6 (Fig. 1) (Laudon et al., 2013). The gauging station C5 is the outlet of Lake Stortjärn (4.2 ha) and has a catchment area of 65 ha. The gauging station C6 is situated 1.5 km downstream of C5 and has a catchment area of 110 ha. The catchment contributing to the C5–C6 reach consists of pine-dominated forest, mostly underlain by post-glacial till soil (72 %). Iron podzols and thin soils can be found in the upland areas, while the shallow subsurface soils of the riparian zone (<1.2 m deep) are dominated by peat. Furthermore, soil wetness and flow accumulation maps based on 2×2 m digital elevation models identify seven wet corridors of 1–10 ha contributing area that extend from upland areas to the stream (Ågren et al., 2014). From those, five are considered DRIPs based on previous field validation using vegetation surveys and thermal and isotopic tracing (Kuglerová et al., 2014; Leach et al., 2017). The five DRIPs collectively account for approximately 60 % of the lateral groundwater inflows along the reach, while the remaining lateral inflow is diffuse (Leach et al., 2017). No tributaries are present along the stream reach, and deep groundwater inflows are minimal in this catchment (Tiwari et al., 2017).

Figure 1 The Stortjärnsbacken catchment in Krycklan, Sweden. The stream reach (blue line) starts at the outlet of Lake Stortjärn (gauging station C5) and ends at the downstream gauging station C6. Stream sampling sites at approximately 50 m increments are indicated with small black squares. Groundwater wells along the reach are indicated with red circles (DRIPs) and orange circles (non-DRIPs). At DRIPs, groundwater flow paths (grey lines) converge in the riparian zone.

The average annual temperature in Krycklan is 1.8 ^∘C (period 1981–2010; Laudon and Ottosson-Löfvenius, 2016), and the average annual precipitation is 614 mm, of which 35 %–50 % falls as snow (Laudon et al., 2013). Approximately 50 % of the annual precipitation translates to streamflow. The hydrological regime at the C5–C6 reach is dominated by the annual snowmelt peak, occurring around May (100–200 L s⁻¹). In summer and fall, low flows (<10 L s⁻¹) alternate with medium to high flows (25–75 L s⁻¹) in response to rain events. During winter, the stream is snow- and ice-covered, with flows < 3 L s⁻¹. At C5, streamflow is mostly driven by lake-level variations. As a result, peak flow events are dampened and recession limbs decline gradually (Leach and Laudon, 2019; Fig. 2). At C6 (1.5 km downstream), streamflow responds much more quickly to hydrological events compared to C5 and is characterized by steep rising limbs (Fig. 2; Ploum et al., 2018).

Figure 2 Hydrographs of the gauging stations C5 (grey) and C6 (blue) during the study period (spring 2017, summer 2017, spring 2018 and spring 2019). The vertical dashed lines and letters correspond to the 14 sampling campaigns. The percentages indicate the net gain in streamflow between the gauging stations C5 and C6.

2.2 Study design, field measurements and laboratory analysis

Field measurements were collected between May 2017 and May 2019. In total, we conducted 14 sampling campaigns with different streamflow conditions, which ranged from drought to peak flow conditions (Fig. 2). Nine sampling campaigns were centered around the snowmelt periods of 2017–2019 and five around a lake damming experiment in summer 2017. In this experiment, the upstream lake was blocked and, after a period of artificial drought, a series of controlled flows were released using a pump. During the course of the artificial drought, the strength of DRIP-stream hydrological connections declined, generating a patchy distribution of lateral DOC inputs similar to those occurring under natural droughts (Gómez-Gener et al., 2020). For each sampling campaign, stream water was collected along the stream reach at approximately 50 m intervals over 1200 m, dividing the stream reach into 25 sections (Lupon et al., 2019). Five of the 25 sections had a DRIP discharging into it, while the other 20 sections were fed by small diffuse groundwater inputs (Figs. 1 and 3). For 10 of those sections, we sampled riparian groundwater inputs from a well network setup, which included five pairs of DRIP and non-DRIP wells located 1–5 m from the stream edge (Ploum et al., 2020). Therefore, we sampled the phreatic groundwater of all DRIPs but not all the diffuse groundwater inputs discharging into 15 reaches. The PVC wells (30 mm diameter) had a mean depth of 95 cm (σ=37 cm) below the soil surface and were fully screened every 5 cm.

Figure 3 Relative upslope contributing area along the stream reach. The solid line represents the UCA model, which assumes that the net gain in streamflow is proportional to the gain in UCA between sampling sites (squares). The dashed line represents the “Diff” model, which assumes uniform, diffuse inflow of groundwater along the entire reach. Grey vertical bars indicate the locations of discrete riparian inflow points (DRIPs) along the stream reach, for which we sampled dissolved organic carbon (DOC) concentrations during the study period and in-stream DOC uptake was considered.

Stream water was collected from the thalweg with acid-washed high-density polyethylene bottles. Groundwater was sampled from PVC wells using suction cup lysimeters and evacuated glass bottles or using a peristaltic pump to fill acid-washed high-density polyethylene bottles. The wells were pre-pumped to ensure we did not sample stagnant water. Bottles for both stream water and groundwater were rinsed three times before filling with minimal headspace. Within 24 h, all samples were filtered (0.45 µm mixed cellulose ester (MCE) syringe filters, Millipore^®) and kept refrigerated at 4 ^∘C until analysis (<7 d after filtering). DOC analysis consisted of acidification of the sample for removing inorganic carbon, followed by combustion using a Shimadzu TOC-V_CPH (analytical error: 2 %; Laudon et al., 2011). The analysis was repeated at least three times per sample, resulting in a DOC concentration (mg L⁻¹) and a percent standard deviation.

2.3 Model framework and data input

We used a mixing model that considered the stream DOC concentration at location i to be a result of upstream DOC flux (location i−1) and the net lateral riparian groundwater flux that is gained along the stream section between locations i−1 and i. In addition, we considered that riparian DOC inputs were subjected to in-stream uptake (Eq. 1).

$$\begin{array}{}\text{(1)}& \begin{array}{rl}& {\mathrm{DOC}}_{\mathrm{sw}i}=\\ & {\displaystyle \frac{\left({\mathrm{DOC}}_{\mathrm{sw},i-\mathrm{1}}\times \left(Q_{i-\mathrm{1}}\right)+{\mathrm{DOC}}_{\mathrm{gw}i}\times Q_i-Q_{i-\mathrm{1}}\right)-{\mathrm{uptake}}_i}{Q_i}},\end{array}\end{array}$$

where DOC_sw,i and ${\mathrm{DOC}}_{\mathrm{sw},i-\mathrm{1}}$ are the stream DOC concentration measured at locations i and i−1, respectively, Q_i and Q_i−1 are the estimated streamflows at locations i and i−1, respectively, the difference (Q_i–Q_i−1) is the net groundwater inflow for that stream section, DOC_gw,i is the estimated groundwater DOC concentration between locations i−1 and i, and uptake_i is the in-stream DOC uptake associated with lateral groundwater labile DOC inputs (see below).

We modified the abovementioned model (Eq. 1) to represent different assumptions about catchment hydrology (diffuse vs. distributed groundwater inputs) and in-stream DOC uptake (no uptake vs. in-stream uptake downstream DRIPs), resulting in four different models (Table 1). The “Diff” model assumed diffuse groundwater inputs and no in-stream DOC uptake. The “Diff-Bio” model also assumed diffuse groundwater inputs but accounted for in-stream DOC uptake downstream from DRIPs. The “UCA” model assumed that groundwater inputs were distributed proportionally to their UCA and no in-stream DOC uptake. Finally, the “UCA-Bio” model assumed groundwater inputs proportional to UCA and accounted for in-stream DOC uptake downstream from DRIPs. Below, we outline the approaches used for estimating riparian groundwater DOC concentrations, groundwater inputs, and in-stream uptake.

2.3.1 Estimates of riparian groundwater DOC concentrations

For each date, we used direct measurements of groundwater DOC concentrations from wells to estimate DOC_gw for 10 sections (i.e., 5 sections with DRIPs and 5 sections without DRIPs). For the remaining 15 sections, we assumed that DOC_gw equaled the average of the non-DRIP wells. For all instances where we did not have direct measurements of DOC_gw and used estimates based on the means of the non-DRIP observations, we also computed the standard deviation of the mean as a measure of uncertainty in the mixing model framework.

2.3.2 Estimates of streamflow and lateral groundwater inputs

Streamflow at each location (Q_i) was represented in the model in two ways. Both approaches assume that all net gain in streamflow between the two hydrological stations C5 and C6 is a result of lateral groundwater input from the riparian zone. One scenario assumed that the local gains in streamflow were driven by diffuse groundwater inflow (hereafter referred to as “Diff”, Fig. 3), where the net gain in streamflow is distributed evenly along the C5–C6 reach:

$$\begin{array}{}\text{(2)}& Q_{\mathrm{diff},i}=\left(Q_{\mathrm{C}\mathrm{6}}-Q_{\mathrm{C}\mathrm{5}}\right)\cdot {\displaystyle \frac{\left(L_i-L_{i-\mathrm{1}}\right)}{L_{\mathrm{total}}}},\end{array}$$

where Q_C5 and Q_C6 are streamflow at gauging stations C5 and C6, respectively (both L s⁻¹), L_i and L_i−1 are the distances between gauging station C5 and sampling locations i and i−1 (both m), respectively, and L_total is the total length of the C5–C6 stream reach (1200 m).

The other scenario (hereafter referred to as “UCA”, Fig. 3) was based on Leach et al. (2017), in which lateral groundwater inputs were distributed proportionally to the gain in upslope contributing area at each stream reach:

$$\begin{array}{}\text{(3)}& Q_{\mathrm{uca},i}=\left(Q_{\mathrm{C}\mathrm{6}}-Q_{\mathrm{C}\mathrm{5}}\right)\cdot {\displaystyle \frac{A_i-A_{i-\mathrm{1}}}{A_{\mathrm{C}\mathrm{6}}-A_{\mathrm{C}\mathrm{5}}}},\end{array}$$

where Q_C5 and Q_C6 are streamflow at gauging stations C5 and C6, respectively (both L s⁻¹), A_i and A_i−1 are the catchment area at locations i and i−1, respectively (both ha), and the difference A_C6–A_C5 is the total gain in catchment area between gauging stations C5 and C6 (55 ha). This approach emphasized the hydrological contributions of DRIPs, because of their large contributing areas relative to the rest of the riparian zone.

2.3.3 Estimates of in-stream DOC uptake

We considered two different scenarios regarding in-stream DOC uptake. One model assumed that all terrestrial DOC inputs were transported to downstream ecosystems (i.e., pulse-shunt concept; no in-stream DOC uptake), while the other model assumed that stream biota rapidly take up the DOC coming from lateral groundwater inputs (i.e., hotspot concept). We did not consider the scenario that DOC coming from the upstream lake was taken up along the stream, as previous studies in the Krycklan catchment have suggested that this DOC is highly recalcitrant and rarely used by stream biota (Tiwari et al., 2014; Kothawala et al., 2015).

At each location (i), in-stream DOC uptake of lateral groundwater DOC inputs (uptake_i, mg C s⁻¹) was estimated as follows:

$$\begin{array}{}\text{(4)}& {\mathrm{uptake}}_i={\mathrm{DOC}}_{\mathrm{gw},i}\cdot V_f/\mathrm{60}\cdot {\mathrm{width}}_i\cdot {\mathrm{length}}_i,\end{array}$$

where DOC_gw,i is the DOC concentration of riparian groundwater (mg L⁻¹), V_f is the DOC uptake velocity (mm min⁻¹) associated with riparian carbon, and width_i and length_i are the mean channel width and the reaction path length of each section, respectively (both m). Based on previous work at this particular study reach, we assumed that in-stream DOC uptake mostly occurred immediately downstream from DRIPs (Lupon et al., 2019). We accounted for this by setting the length of all sections to zero, except for those where a DRIP was located (Fig. 3). At these sections, length_i was the distance between DRIPs and location i instead of the total length between i−1 and i. This prevented overestimations of reaction times and path lengths over which in-stream uptake took place. For in-stream DOC uptake from riparian groundwater, we used $V_f=\mathrm{0.6}\pm \mathrm{0.06}$ mm min⁻¹. This value is the median V_f for DOC reported in the literature and has been shown to realistically simulate in-stream DOC uptake at a whole river network (Mineau et al., 2016). Because V_f depends on temperature, streamflow, DOC composition, and microbial assemblages, we tested values for V_f ranging between 0.25 and 1.11 mm min⁻¹. These values yielded similar model results for the simulations that considered in-stream DOC uptake.

2.4 Model uncertainty and performance criteria

For each model, we accounted for uncertainty in modeled stream DOC concentrations (Eq. 1) by incorporating errors in Q observations, water sample analysis, and estimates of both DOC_gw and DOC V_f. We assumed normally distributed errors for Q_c5 and Q_c6 (±10 %, based on repeat streamflow gauging; Karlsen et al., 2016), DOC_gw (either ±2 % for sites with measurements based on laboratory analytical precision or ±1 SD – standard deviation – of the mean for sections that rely on estimates), and V_f (±10 %, based on Mineau et al., 2016). For each date, each model was run 10 000 times using random values selected from these parameter distributions. Error estimates in DOC_sw,i were tracked downstream since the values become ${\mathrm{DOC}}_{\mathrm{sw},i-\mathrm{1}}$ in the computation for the next stream reach.

We evaluated the simulation of each run using two goodness-of-fit metrics, computed using the “hydroGOF” R package (Zambrano-Bigiarini, 2020). First, we computed the percent bias (PBias, %), which measures the average tendency of the simulated values to be larger (PBias > 0 %) or smaller (PBias < 0 %) than their observed ones. We considered that a model successfully simulated the magnitude of stream DOC concentrations if the median value of all the runs was within −5 % and +5 % bias. Second, we calculated the Spearman correlation (R), which shows whether the longitudinal patterns of simulated DOC concentrations mimicked the observed ones. In this case, we considered that a model was capturing the general direction of stream DOC concentrations if the median R of all runs was higher than 0.70.

Figure 4 Longitudinal patterns of DOC concentrations along the C5–C6 reach. Each panel, indicated by label and date, shows one sampling campaign. The black dots are the observed stream DOC concentrations. The colored bands show the simulations of the four models. The vertical grey lines show the locations of DRIPs with wells (solid) and without wells (dashed). The streamflows (Q) at gauging stations C5 and C6 are shown for each sampling campaign.

3 Results

3.1 Stream hydrology

Across the 14 sampling campaigns, hourly Q ranged from 0 to 116 L s⁻¹ and from 2 to 152 L s⁻¹ at C5 and C6, respectively, which were comparable to the range of flows observed throughout the whole ice-free period (C5: 0–150 L s⁻¹, C6: 2–200 L s⁻¹; Fig. 2). At both gauging stations, maximum Q occurred during the snowmelt, whereas minimum Q was observed during the artificial drought in summer 2017. As a result, the net gain in streamflow along the C5–C6 reach ranged from 8 % (artificial flood, event G) to 90 % (artificial drought, event H) (Fig. 2). During the other sampling campaigns, the net gain in streamflow along the reach was between 20 % and 50 %, with a mean of 37 % (Fig. 2).

3.2 Stream DOC concentrations

During the study period, stream DOC concentration ranged from 15 to 32 mg L⁻¹ and varied over time as well as along the stream reach (Fig. 4). Seasonally, average DOC concentration decreased as the snowmelt period progressed during each year – events A–E (2017), J–L (2018), and N–M (2019). In summer 2017, stream DOC concentrations were relatively constant at C5 (19-20 m L⁻¹), whereas they decreased during the same period at C6 (from 28 to 18 mg L⁻¹) (events F–I). Spatially, stream DOC concentrations generally decreased along the C5–C6 reach (8 out of 14 sampling campaigns) (Fig. 4). However, stream DOC concentrations clearly increased along the reach for the increasing limb of snowmelt 2017 (event B), the summer storm event (event F) and the artificial drought (event G). During the recession limb of the snowmelt peak 2017 (event D) and the lake flooding experiment (events H and I), stream DOC concentrations were relatively constant along the reach.

Abrupt changes in stream DOC concentrations occur in those sections affected by DRIPs (Fig. 4). During most snowmelt campaigns (events A–C, J–N) and summer baseflow conditions (event E), stream DOC concentrations sharply decreased in sections fed by DRIPs. The only exception was the section affected by the last DRIP, from which stream DOC concentrations tended to increase. Peaks in DOC concentrations immediately downstream from DRIPs also occurred during the summer rain event and the experimental drought (events F and G). For the other sampling campaigns (events D, H, and I), both increases and decreases in DOC concentrations occurred at DRIP locations.

Figure 5 Relative bias (PBias) by model and sampling campaign. For each model, boxplots show the median and 25th and 75th percentiles, and whiskers extend to the 10th and 90th percentiles of the 10 000 runs. Values close to 0 indicate that the model successfully simulates the magnitude of stream DOC concentrations. The horizontal line at PBias = 0 is shown as a reference.

3.3 Model simulations

The ability of the models to simulate the magnitude and longitudinal patterns of stream DOC concentrations varied across sampling campaigns (Fig. 4). For most sampling campaigns performed during the snowmelt period, at least one of the models was able to capture either the magnitude or spatial variations in stream DOC concentrations (Fig. 4). For six events (B–D, K, and M–N), all the models captured the magnitude of stream DOC concentrations (median PBias from −5 % to 5 %; Fig. 5), yet only for event K were all of them also able to simulate their longitudinal patterns (median R>0.80; Fig. 6). Indeed, none of the models captured spatial patterns for events B, D, and M (median R<0.50), although the Diff-Bio model tended to perform better than the others (Fig. 6). For event C, spatial patterns were captured by both the Diff-Bio and UCA-Bio models (median R∼0.75; Fig. 6), whereas models UCA and UCA-Bio were able to simulate patterns of stream DOC concentrations for event N (median R∼0.70, Fig. 6). For the three other field campaigns performed during the snowmelt period (events A, J, and L), the magnitude of stream DOC concentrations was only successfully captured by models UCA and/or UCA-Bio (median PBias ∼ 0 %, Fig. 5), despite all models being able to simulate the spatial patterns of DOC concentrations (R>0.70; Fig. 6).

Figure 6 Spearman regression (R) by model and sampling campaign. For each model, boxplots show the median and 25th and 75th percentiles, and whiskers extend to the 10th and 90th percentiles of the 10 000 runs. Values close to 1 indicate that the model successfully simulates the magnitude of stream DOC concentrations. The horizontal line at R=0 is shown as a reference.

For the sampling campaigns performed during summer 2017 (events E–I), there were large inconsistencies across models (Fig. 4). For summer baseflow conditions (event E), the Diff-Bio and UCA-Bio models successfully simulated both the magnitude (PBias < 3 %) and spatial patterns (R>0.90) of stream DOC concentrations (Figs. 5 and 6). For the natural rain event (event F), all the models underestimated stream DOC concentrations (PBias < −5 %) and also failed to predict their spatial variation (median R<0.35). Yet the UCA model performed better than the others (Figs. 5 and 6). None of the models was able either to capture spatial patterns during the lake flooding (event H, R<0.5), even though all the models captured the overall magnitude of DOC concentrations (Figs. 5 and 6). For the experimental drought (event G), the Diff and UCA models successfully simulated spatial patterns of stream DOC concentrations (R>0.70; Fig. 5), yet only the Diff model accurately captured their magnitude (median PBias = 1 %; Fig. 6). Similarly, only the Diff model was able to simulate both the magnitude and spatial pattern of stream DOC concentrations for the post-flooding campaign (event I; Figs. 4–6).

4 Discussion

4.1 Modeling spatial patterns of stream DOC concentration

DRIPs are an important, and often primary, source of water and C to headwaters (Briggs and Hare, 2018; Demars, 2019) and play a major role in regulating spatial variation in stream DOC concentration, processing, and export. Our spatially explicit surveys revealed that longitudinal patterns of stream DOC concentrations varied across flow conditions. In general, DOC concentrations tended to decrease along the C5–C6 reach, indicating that DOC was generally diluted or taken up along the stream segment. However, we observed step changes in stream DOC concentrations at DRIPs, indicating that these locations play a key role in stream C export. Similar patterns have been observed in boreal headwaters (Duvert et al., 2018; Lupon et al., 2019), yet our study is the first to reveal the mechanisms by which DRIPs shape spatial patterns of DOC concentrations and fluxes in these streams.

Our results show that accounting for spatial variability in lateral groundwater inflows in the models (i.e., UCA) improved simulations of stream DOC concentrations for five out of the nine sampling campaigns performed during the snowmelt period (events A, C, J, L, and N), when groundwater inputs were high (>20 L s⁻¹) and/or contributed significantly (>40 %) to streamflow. During these events, sharp decreases in DOC concentrations were observed in sections fed by DRIPs, suggesting that these preferential groundwater flow paths mostly diluted the DOC concentrations in the stream. Further, UCA improved the simulations of stream DOC concentrations during the summer rain (event F), when the net gain in streamflow along the reach was 46 %. In this case, however, DOC concentrations increased in most DRIP locations, indicating that these flow paths were acting as important sources of C to the stream. Previous studies have observed that high groundwater tables associated with rain events often increase groundwater DOC concentrations by activating the dominant source layer (Ledesma et al., 2018), which might explain the observed increase in DOC concentrations along the reach observed for event F. Regardless of the process (i.e., DOC delivery or dilution), these findings corroborate that the spatial variability in groundwater flow paths related to landscape topography has a major influence on stream C patterns when streams are mostly fed by groundwater flow (Covino et al., 2021; Dupas et al., 2021; Rocher-Ros et al., 2019).

Accounting for spatial variability in groundwater inflow was important during some, but not all, events. Representing UCAs did not improve model simulations for the sampling campaigns close to a snowmelt peak (events B, K, and M) despite groundwater inputs being elevated (25–36 L s⁻¹). Our explanation is that, during these events, increases in the groundwater level might homogenize groundwater inflows along the reach, potentially generating overland flow because of soil frost causing impervious conditions (Ploum et al., 2020). Similarly, a potential homogenization of overland flow during the snowmelt can explain the decline in DOC concentrations during events K and M (Ploum et al., 2018; Laudon et al., 2011). In any case, from our work it is evident that model frameworks that integrate the spatial arrangement of groundwater flow paths (i.e., DRIPs) can help represent the variability in the hydrological connectivity along the stream.

Model simulations also revealed that accounting for in-stream uptake downstream of DRIPs improved predictions of longitudinal patterns of stream DOC concentrations for the five sampling campaigns occurring from May to August (events C–E, H, and M), but especially during summer low-flow conditions (event E). It is likely that these results are explained by the seasonal pattern of microbial activity in boreal streams, which often mirror the temporal variation in water temperature (Burrows et al., 2017). However, in-stream uptake did not improve model simulations during those dates in summer characterized by very low (event G) or high (event I) flows. These results concur with other recent studies (Lupon et al., 2019; Seybold and McGlynn, 2018; Demars, 2019) and suggest that aquatic biological activity is enhanced at the transition between low and high flows due to increases in labile DOC supply from terrestrial systems. Conversely, in-stream DOC uptake may be minimal during low flows due to C limitation (Burrows et al., 2017) or overwhelmed by high water velocities during rain events (pulse-shunt concept; Raymond et al., 2016). Most importantly, our findings support the idea that spatial patterns of DOC concentrations along headwaters are associated with aquatic biological activity and stream water permanence driven by terrestrial flow path organization (Hale and Godsey, 2019). Collectively, these results suggest that DRIPs are important hotspots of DOC uptake and, thus, the capacity for processing DOC of boreal headwater streams is closely tied to the spatial arrangement of lateral inputs of DOC from riparian zones.

4.2 Limitations of the model

Our model framework represented the source of lateral DOC inputs based on groundwater samples from a riparian well network that compared DRIP and non-DRIP groundwater chemistry (Ploum et al., 2020). This allowed us to distinguish between the spatial variability in riparian groundwater chemistry associated with different soil wetness regimes (Vidon, 2017). For example, during the experimental drought (event G), the Diff model provided better simulations of both the magnitude and spatial patterns of stream DOC concentrations compared to the assumption of uniform inputs along the reach, suggesting that the representation of spatial variability in groundwater DOC concentrations was more important than hydrology or in-stream uptake. Hence, under these conditions, stream DOC patterns might not be directly related to groundwater fluxes but rather to the thermal and chemical conditions that groundwater discharge creates at the local level (Briggs and Hare, 2018). However, there are also limitations in our groundwater sampling approach. For example, our groundwater sampling was not able to represent temporal DOC dynamics associated with variability in groundwater travel times (Heidbüchel et al., 2020), event-scale variability in riparian DOC mobilization (Werner et al., 2019), or the activation of DOC from different soil layers (Ledesma et al., 2018).

Apart from the limitations of our groundwater sampling, we identified some limitations in the hydrological and biogeochemical components of the model as well. The sampling campaign of the summer rain event (event F) is a clear example of mismatch between our simulations and the observations. During this event, stream DOC concentrations increased along the reach, but most of our models simulated a decreasing pattern. Similarly, none of the models properly simulated the longitudinal patterns for two sampling campaigns performed around the snowmelt peak (events B and M). These examples suggest that the model framework has some limitations in representing the complex hydrologic and biogeochemical dynamics occurring in headwater streams (Ambroise, 2004; Klaus and Jackson, 2018). For instance, our models do not account for local conditions affecting snowmelt rates on hillslopes (i.e., shading, sun exposure) or local variations in precipitation, interception, or infiltration that are relevant during rain events (Laudon et al., 2004; Lyon et al., 2010). For the biogeochemical component of our model, we did not take into account processes that produce (i.e., resuspension) or remove (i.e., photodegradation, sorption, flocculation) DOC from the water column (Droppo et al., 1998; Kaiser and Kalbitz, 2012). Further, in-stream DOC uptake was assumed to occur only downstream of DRIPs and at a uniform rate across flow conditions. Previous studies have shown that uptake rates can vary over time as a function of temperature, DOC composition, and microbial assemblages (Berggren et al., 2009; Mineau et al., 2016). While the use of other values for V_f generally resulted in similar model output (Fig. 5 and 6), we cannot rule out the idea that V_f varied among DRIPs and/or over time due to changes in groundwater DOC composition and temperature. To better understand the role of DRIPs in stream hydrology and biogeochemistry, future empirical studies testing how DRIPs affect specific processes are needed. Nevertheless, our study, even with its limitations, demonstrated that both lateral discrete and diffused inputs as well as biological activity are essential components of the DOC patterns in boreal streams. These findings shed new light on the understanding of C dynamics across the boreal aquatic–terrestrial interface.

Another major limitation of our models is their large uncertainty, especially during events with large groundwater contributions such as event B. For six events (A–C, G, J, and L), all the models showed large inconsistencies among runs, resulting in simulated DOC concentrations at C6 that vary over 10 mg L⁻¹. Moreover, the uncertainty in groundwater DOC concentrations was large, because not all stream sections were sampled and groundwater inputs of DOC had to be estimated based on means of the available DOC concentrations from non-DRIP wells. For future studies, we have identified two more directions that can be useful for improving the simulations of stream DOC dynamics along boreal headwaters. For the representation of the spatial heterogeneity in riparian hydrochemistry, the hydrological representation of lateral groundwater inputs through the distinction of DRIP and non-DRIP riparian zones can be further developed. For this matter, integrative hydrochemical frameworks that represent fluxes from various soil layers would be useful to include, especially at non-DRIPs, because here groundwater levels are more dynamic compared to DRIPs (Seibert et al., 2009; Ploum et al., 2020). Furthermore, it can be of interest to downscale the number of riparian groundwater chemistry samples to understand which minimum set of samples is required to represent the spatial heterogeneity in sources of lateral DOC inputs from riparian zones to streams. A preliminary analysis indicated that the most optimal strategy for reducing model uncertainty was to monitor DRIPs individually while averaging DOC concentrations at non-DRIPs (Ploum, 2021). However, given that non-DRIP groundwater chemistry changes with groundwater table fluctuations (Ledesma et al., 2015; Ploum et al., 2020), it is likely that optimizing groundwater sampling campaigns requires careful consideration of the antecedent groundwater conditions.

5 Conclusions

This study provides new insight into the role of DRIPs in stream DOC concentrations in boreal headwater catchments. We showed that DRIPs influence longitudinal patterns of stream DOC concentrations at small spatial scales (a few meters) by controlling both the hydrology and biogeochemistry of the streams they feed. However, our study also shows that the role of DRIPs can change over time depending on hydrologic conditions. During high flows, DRIPs control DOC concentrations by diluting upstream DOC. In contrast, in late spring and summer, DRIPs can be important sources of C for stream biota, delivering labile resources from their upstream contributing areas (UCAs) and promoting local hotspots of in-stream DOC uptake downstream confluences. These results suggest that future changes in catchment hydrology associated with global change can affect DOC exports from boreal fluvial networks by shifting the dominant mechanisms by which DRIPs drive spatial patterns of DOC concentrations and processing along headwater streams. Thus, the identification and characterization of DRIPs are essential for understanding the current and future mechanisms behind C fluxes from boreal fluvial networks.