A parsimonious model of longitudinal stream DOC patterns based on groundwater inputs and in-stream uptake

The supply of terrestrial dissolved organic carbon (DOC) to aquatic ecosystems affects local in-stream processes and downstream transport of DOC in the fluvial network. However, we have an incomplete understanding on how terrestrial DOC inputs alter longitudinal variations of DOC concentration along headwater stream reaches because groundwater discharge, 15 groundwater DOC concentration and in-stream DOC uptake vary at relatively short spatial and temporal scales. In the riparian zone, the convergence of subsurface flow paths can facilitate the inflow of terrestrial DOC from large upslope contributing areas to narrow sections of the stream. We refer to these areas of flow path convergence as discrete riparian inflow points (DRIPs). In this study, we ask how longitudinal patterns of stream DOC concentrations are affected by DRIPs, as they are major inputs of terrestrial DOC and important locations for in-stream processes. We used a mixing model to simulate stream 20 DOC concentrations along a 1.5 km headwater reach for fifteen sampling campaigns with flow conditions ranging from droughts to floods. Four sets of model scenarios were used to compare different representations of hydrology (distributed inputs of DRIPs vs diffuse groundwater inflow), and in-stream processes (passive transport vs in-stream biological uptake). Results showed that under medium (10-50 l/s) and high flow conditions (>50 l/s), accounting for lateral groundwater inputs from DRIPs improved simulations of stream DOC concentrations along the reach. Moreover, in-stream biological uptake 25 improved simulations across low to medium flow conditions (< 50 l/s). Only during an experimental drought, longitudinal patterns of stream DOC concentration were simulated best using diffuse groundwater inflow and passive transport scenarios. These results show that the role of hydrology and in-stream processes on modulating downstream DOC exports varies over time. Importantly, we demonstrate that accounting for preferential groundwater inputs to the stream is needed to capture longitudinal dynamics in mobilization and in-stream uptake of terrestrial DOC. The dominant role of DRIPs in these transport 30 and reaction mechanisms suggests that consideration of DRIPs can improve stream biogeochemistry frameworks and help inform management of near-stream areas that exert a large influence on stream conditions. https://doi.org/10.5194/hess-2021-358 Preprint. Discussion started: 28 July 2021 c © Author(s) 2021. CC BY 4.0 License.


Introduction
Running waters play a critical role in the global carbon (C) cycle because they transport dissolved organic carbon (DOC) from 35 lands to oceans (Cole et al., 2007). Accounting for stream DOC fluxes within stream networks can be valuable for understanding net in-stream C retention (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 (Kritzberg et al., 2020).
Across fluvial networks, headwater and high-order river segments have different controls of C dynamics. In headwaters, terrestrial source-transport mechanisms are considered dominant drivers of DOC dynamics (Creed et al., 2015). The 40 organization of flow paths to riparian zones (RZ) (Jencso et al., 2010;McGlynn and McDonnell, 2003;Ploum et al., 2020), riparian wetness regimes (Vidon, 2017), and local differences in soil organic carbon stocks (Grabs et al., 2012) are major controls of the spatiotemporal variability of terrestrial DOC fluxes to streams. Further, recent studies have shown that headwaters can have a reactive role as well, which can reduce downstream supply of DOC Hotchkiss et al., 2015). Mineralization by biota is a major removal mechanism of terrestrial DOC, in addition to abiotic processes such as 45 adsorption, flocculation, and photooxidation (Mineau et al., 2016). Since the upstream supply influences downstream DOC dynamics, the entire stream network relies on terrestrial DOC sources and fluxes from RZ to headwaters (Raymond et al., 2016). Therefore, it is important to understand where and when RZs hydrologically connect to headwaters, and what effect this connection has on headwaters stream biogeochemistry (Casas-Ruiz et al., 2017).

50
To integrate the supply of terrestrial DOC fluxes and removal in aquatic ecosystems, we need to combine source-transport hydrochemical and in-stream C spiraling frameworks (Li et al., 2020). While there are hydrochemical frameworks that integrate spatial heterogeneity of terrestrial DOC fluxes to streams (Seibert et al., 2009), they are mostly suitable during high flow conditions, when stream DOC dynamics are dominated by source-transport mechanisms, and in-stream uptake processes are less significant (Raymond et al., 2016). By contrast, in-stream processes are mostly measured during low flows and within 55 reaches with no lateral inputs (Mineau et al., 2016), and therefore, little is known about the role of in-stream biota on controlling DOC fluxes during storms or snow melt periods (Talbot et al., 2018;Wollheim et al., 2018). Moreover, the few quantitative biogeochemical frameworks that consider both source-transport and in-stream processing often rely on a sparsely distributed gauging stations, which do not fully capture the spatiotemporal variability of DOC dynamics within the fluvial network (Futter et al., 2007). A major influence on local stream DOC dynamics that is not captured by gauging stations is the spatial distribution 60 of groundwater discharge (Briggs and Hare, 2018). Discrete riparian inflow points (DRIPs) convey substantial fluxes of DOC from large upslope contributing areas to narrow sections of the stream (Ploum et al., 2020). Subsequently the labile DOC that DRIPs convey to headwater streams leads to in-stream uptake directly downstream of DRIPs (Lupon et al., 2019). Depending on physical stream conditions and local microbial communities (Berggren et al., 2009;Kothawala et al., 2015;Wollheim et al., 2015), this means that DRIPs are important conveyors of terrestrial DOC to streams, but also potential hotspots of in-65 stream DOC uptake.
To represent local evasion of C from streams to the atmosphere in C budgets, it is important to consider landscape features such as DRIPs (Rocher-Ros et al., 2019). To identify and characterize the groundwater supply from DRIPs to streams, there is a growing body of research that combines field-based methods with digital elevation models, and remote sensing approaches 70 (Antonelli et al., 2020;Barclay et al., 2020;Briggs et al., 2019;Leach et al., 2017;Lidberg et al., 2019;Rosenberry et al., 2016;Selker et al., 2006). However, to quantify the supply of terrestrial DOC from DRIPs to stream networks and to consider subsequent in-stream uptake, it is important to characterize the spatial variability in groundwater chemistry as well. This remains an ongoing challenge because groundwater chemistry is highly variable in space and time (Kiewiet et al., 2019), and groundwater dynamics often need to be extrapolated from a relative small number of observations (Rinderer et al., 2019). As 75 such, the way groundwater and streams are currently monitored does not facilitate a systematic consideration of source- In this study, we integrate source-transport mechanisms and in-stream uptake in a model framework, by considering DRIPs as the primary driver of both hydrological and biogeochemical processes that influence patterns of stream DOC concentrations.
We used a mixing model to simulate stream DOC concentrations along a 1.5 km headwater reach, for fifteen sampling campaigns with flow conditions ranging from droughts to floods. With this modelling approach, our objective was to 85 disentangle the control of terrestrial DOC inputs and in-stream DOC uptake on longitudinal stream DOC concentrations during different flow conditions. In different model scenarios, we accounted for two types of transport mechanisms: 1) assuming uniform, diffuse inflow of groundwater along the reach and 2) by assuming groundwater inflow relative to upslope contributing area (UCA). These two assumptions were than combined either with the assumption that the stream is a passive pipe, or with the assumption that in-stream DOC uptake takes place directly downstream of DRIPs. 90 2 Methods

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 segment located between the gauging stations C5 and C6 ( Fig. 1)  . The gauging station C5 is the outlet of lake Stortjärn (4.2 ha), with a catchment area of 65 ha. The gauging station C6 is situated 1.5 km downstream, and has a catchment 95 area of 110 ha. The catchment 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, and in the RZ the shallow subsurface (<1.2 meter) is dominated by peat.
For the period 1981 -2010, the average temperature was 1.8 °C and the average annual precipitation in Krycklan was 614 mm, of which 35-50% fell as snow (Laudon and Ottosson Löfvenius, 2016;Laudon et al., 2013). Recently, it has been reported 100 that winter and fall temperatures are increasing, and that snow cover is decreasing (Laudon et al., 2021). On average, approximately 50% of the annual precipitation translates to streamflow. The hydrological regime at gauging stations C5 and C6 is dominated by the annual snowmelt peak, occurring around May (100-200 l/s), but at C5 peak flows are dampened by the lake . In summer and autumn, low flows dominate (5-10 l/s), but are alternated with medium to high flows (25-75 l/s) as a response to rain events. During the winter and early spring, the stream is snow and ice covered with 105 flows around 1-2 l/s. At C5, streamflow is mostly driven by lake level variations as a response to water contributed from the surrounding mire complex. As a result, peak flow events are dampened and recession limbs decrease slowly (Fig. S1). The downstream effect of the lake is a dominant aspect of the discharge and chemistry at C6 . At C6, discharge can respond much faster to hydrological events compared to C5 and is characterized by steep rising limbs (Ploum et al., 2018). This hydrological response is typically a result of rapid increases of shallow groundwater tables in the RZ, 110 enabling periods of flow through highly conducting organic soil layers in the upper decimetres of the soil profile (Ledesma et al., 2015). The lateral inputs from the RZ are dominated by discrete riparian inflow points (DRIPs), which route 60% of the upslope contributing areas to 5% of the stream length (Leach et al., 2017). At DRIPs, flow paths from the upland converge in the RZ, which leads to near-surface groundwater levels and organic-rich groundwater chemistry (Ploum et al., 2020).

Field measurements and laboratory analysis
Field measurements were collected between May 2017 and May 2019. In total, we conducted 15 sampling campaigns with different streamflow conditions, which ranged from drought conditions to peak flows (Fig. S1). Ten sampling campaigns were centred around the snowmelt periods of 2017-2019, and five around a lake damming experiment in August 2017. In this 125 experiment, the upstream lake was blocked, and after a period of artificial drought a series of controlled flows were released using a pump (Gómez-Gener et al., 2020).
For each sampling campaign, stream water was collected along the stream segment at approximately 50 meter intervals over 1200 meters, dividing the stream into 25 reaches (Lupon et al., 2019). Riparian groundwater was sampled from a paired well 130 network setup, which included pairs of DRIP and non-DRIP wells (Ploum et al., 2020). The fully screened, PVC wells (30 mm diameter) were positioned 1-5 m from the stream edge, and had a depth of approximately 1 m. Five of the stream reaches were directly associated with a DRIP groundwater well, and five directly with a non-DRIP well. Two additional DRIPs were identified after installation of GW wells, and are therefore without a well and groundwater samples (Fig. 2, dashed vertical lines). At one stream reach, two DRIPs are located (Fig. 2, marked with 513). 135 Stream water was collected from the talweg with acid-washed high-density polyethylene bottles. Groundwater wells were sampled using suction cup lysimeters and vacuumed glass bottles, or by using a peristaltic pump to fill acid-washed highdensity polyethylene bottles. The wells were pre-pumped to ensure we did not sample stagnant water. Stream water was collected using grab samples. Bottles for both stream water and groundwater were rinsed three times before filling with 140 minimal headspace. Within 24 hours, all samples were filtered (0.45 µm) and kept refrigerated at 4 °C before analysis. DOC analysis consisted of acidification of the sample, followed by combustion using a Shimadzu TOC-VCPH. The analysis was repeated at least three times per sample resulting in a DOC concentration in mg/l and a percent standard deviation. https://doi.org/10.5194/hess-2021-358 Preprint. Discussion started: 28 July 2021 c Author(s) 2021. CC BY 4.0 License. Figure 2 Relative upslope contributing area along the stream reach. The solid line is the relative gain in upslope contributing area between each of the stream sampling sites (triangles), which referred to as UCA model. The dashed line represents the DIFF model, which assumes uniform, diffuse inflow of groundwater along the entire reach. Green vertical bars indicate locations of discrete riparian inflow points equipped by GW wells (DRIPs; solid) and additional DRIP sites without GW wells (dashed). In-stream dissolved organic carbon (DOC) uptake was considered only in those reaches comprising DRIPs (i.e. after the vertical bars).

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 and the net lateral riparian groundwater flux that is gained between upstream location i-1 and i. In addition, we considered that riparian DOC inputs were subjected to in-stream uptake. 155 where , and , are the stream DOC concentration at location i and i -1, respectively; and are the streamflows at locations i and i-1, respectively; , is the groundwater DOC concentration at location i, and 160 is the in-stream DOC uptake associated with lateral groundwater labile DOC inputs (see below). We used this equation in combination with different assumptions for the terms that represent streamflow and in-stream uptake.

Riparian groundwater DOC concentrations
Out of the 25 stream reaches, there were ten reaches where DOCgw was directly obtained from DRIP and non-DRIP wells.
From the remaining 15 reaches, 13 were classified as non-DRIP zones and were represented by the chemistry of the nearest 165 non-DRIP wells. When DOCgw in non-DRIP wells were not available, their DOCgw was weighted based on the mean of all other non-DRIP wells. The remaining 2 stream reaches were classified as DRIPs because these stream reaches had comparable UCA to DRIPs: both reaches exceeded the 75 th percentile of the gained upslope contributing area (UCA) relative to catchment area, similar to DRIP reaches. In these two cases, DOCgw was considered as the mean of all DRIP wells (Fig. 2, dashed vertical bars). Lastly, the DOCgw at the stream location indicated with 513 in Figure 2 had two DRIPs flowing into the stream at the 170 same location. Here DOCgw was a weighted average based on the DOCgw of both DRIP wells. For the averaging, the UCA of the two DRIP wells relative to the total gain in UCA at that location was used. On 2018-05-11, DOCgw of both wells was not available, and therefore we considered the mean of the other DRIP wells.

Streamflow and lateral groundwater inputs 175
Streamflow at each location ( ) was considered in two ways. Both approaches assume that all net gain in streamflow between the two hydrological stations C5 and C6 is a result of shallow lateral groundwater input from the RZ. One scenario assumed that the local gains in streamflow were driven by diffuse groundwater inflow (hereafter referred as "DIFF", Fig. 2), where the net gain in streamflow is distributed according to gained stream distance: where and is the streamflow at the gauging stations C5 and C6, respectively (both in l/s). is the total length of the C5-C6 stream segment (1200 m), is distance between C5 and the sampling location i, and the distance of the upslope stream sampling location. 185 The other scenario (hereafter referred as "UCA", Fig. 2) was based on Leach et al. (2017), in which lateral groundwater inputs were distributed proportional to the gain in upslope contributing area at each riparian reach: where is the upslope contributing area that is gained along the reach (between locations i and i-1), which was used to distribute the net gained streamflow ( − ) proportional to the total gain in catchment area between the lake and the downstream outlet ( − ). This approach emphasized the hydrological contributions of DRIPs, because of their large contributing areas relative to the rest of the RZ. 195

In-stream DOC processing
We considered two different scenarios regarding in-stream DOC processing. One scenario considered that the stream was acting as a passive pipe for DOC (i.e. no in-stream DOC uptake), while the other scenario considered that stream biota rapidly take up the DOC originating from lateral groundwater inputs. We did not consider the scenario that DOC coming from upstream locations were taken up along the stream, as previous studies in Krycklan have suggested that this rarely occurs 200 (Tiwari et al., 2014).
At each location (i), in-stream DOC uptake of lateral groundwater DOC inputs ( , in mg C/s) was estimated as follows: where , is the DOC concentration of riparian groundwater (in mg/l), the DOC uptake velocity (in mm/min) associated with riparian carbon, and ℎ and ℎ are the mean channel width and the reaction path length of each reach (both in m). Based on previous work at this particular study segment, we assumed that in-stream DOC uptake mostly occurred immediately downstream of DRIPs (Lupon et al., 2019). We accounted for this by setting the ℎ of all reaches 210 to zero, except for those where a DRIP was located (Fig. 2). At these reaches, ℎ was the distance between DRIPs and the 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 a = 0.6±0.06 mm/min. This value is the median ambient DOC obtained from a literature review and has been shown to realistically simulate in-stream DOC uptake at whole river networks (Mineau et al. 2016). Because depends on temperature, 215 streamflow, DOC composition, and microbial assemblages, we tested values for ranging between 0.25 and 1.11 mm/min yielded. These values yielded similar model results for the simulations that considered in-stream DOC uptake (Fig. S2).

Model scenarios, uncertainties and performance criteria
In total, our model approach resulted in four scenarios per sampling campaign (Table 1). For each sampling date, the model scenarios were informed with the same DOCgw concentrations, but assumed different representations of catchment hydrology 220 and in-stream DOC processing. "DIFF_NOBIO" assumes diffuse groundwater inputs and passive transport of DOC in the stream. "DIFF_BIO" also assumes diffuse groundwater inputs, but accounts for in-stream DOC uptake downstream from DRIPs. The "UCA_NOBIO" model assumes that groundwater inputs are distributed proportional to their UCA and no instream DOC uptake, which means that DRIPs contribute a relatively large component of the net gained stream water compared to non-DRIP stream sections. The "UCA_BIO" model assumes groundwater inputs proportional to UCA as well, and accounts 225 for in-stream DOC uptake downstream from DRIPs. To account for uncertainties in our streamflow observations, water sample analysis and model assumptions, we used a stochastic approach to evaluate our model simulations. This means that we compared the range of uncertainty in our simulations to the range of uncertainty in the observations. For each scenario we executed a 100 model runs, with each run using a randomized sample for every observation and every computation. We presented the simulated , using an 235 uncertainty band that represented 66% of the randomized model runs (mean +/-one standard deviation). The randomizations for each sample were based on normal distributions considering 10% uncertainty of streamflow observations and DOC uptake velocity, and the percentage standard deviation from the laboratory analysis (for most samples this was <2%). For each run, the model generates a single ("deterministic") value for , . This means that for the next downstream location the gained uncertainty (from the streamflow, DOC uptake velocity and laboratory analysis) is incorporated as a single value of 240 , . As such, DRIPs can have large influence on the uncertainty band and sampling days with relative large gains in streamflow from groundwater inputs can have increasing uncertainty bands further downstream. In the case that upstream from location i the uncertainty band in was large, but the uncertainty of the lateral input at location i was small, the uncertainty band can become narrower.

245
We evaluated the simulations using a number of goodness of fit metrics, computed using the "hydroGOF" R-package (Zambrano-Bigiarini, 2020). We computed the root mean squared error (RMSE), percent bias (PBIAS), Nash-Sutcliffe efficiency (NSE), and coefficient of determination (R 2 ) for each sampling. RMSE aggregates the magnitude of errors in mg/l.
We considered a 2 mg/l magnitude of error acceptable. PBIAS shows systematic errors in the simulations as a percent deviating from the observations, where positive values indicate overestimation and negative values underestimation. We considered a 250 5% bias acceptable. NSE ranges from -∞ to 1, and shows how the variance of the simulation corresponds to the variance in observations. R 2 compares variances as well, but the simulated and observed variance is considered a cloud of points, which does not account for the ability of to simulate patterns in the observations. For both NSE and R 2 , we considered that values > 0.5 were indicating that the model was capturing the general direction and variance of stream DOC concentrations. For R 2 , these assumption was only true if the relation between observed and simulated stream DOC concentrations was positive. 255

Results
During the study period, the net gain in streamflow along the segment ranged between 8% and 90%, with an average of 37% ( Fig. S1). Stream DOC concentration ranged between 15 to 32 mg/l and generally decreased along the C5-C6 stream segment (11 out of 15 sampling campaigns) (Fig. 3). The variability in stream DOC concentration along the segment and across different flow conditions was reasonably well captured by all models because, although the models rarely captured the entire 260 longitudinal pattern, the order of magnitude of the simulations was within the range of the observations, and the general direction of the simulated values corresponded with the observations. For the majority of the simulations, RMSE were in the order of 2 mg/l and PBIAS was within a range of -5% to 5% (Table 2). Further, for 10 out of 15 sampling campaigns, at least one model had an R 2 greater than 0.5 (Table 2). However, only in four sampling campaigns, the NSE of at least one model was greater than 0.5 (Table 2, Fig. 3F, K-M). 265

Snowmelt 2017
The sampling campaigns of spring 2017 captured an early snowmelt peak, and the receding limb of a rain-on-snow event ( Fig   3A-E). At the onset of the snowmelt peak, stream DOC concentration ranged between 25 and 30 mg/l and showed a marked longitudinal pattern (Fig. 3A-B). In the first sampling campaign, the two scenarios assuming diffusive groundwater inputs captured the decreasing trend of stream DOC concentration in the first 700 m, but underestimated the concentrations at the 270 downstream section of the stream segment (Fig. 3A). In the second sampling campaign, none of the models correctly captured the longitudinal patterns of stream DOC concentrations (NSE<0), as all them simulated uniform concentrations around 27 mg/l (Fig. 3B). Around the time of the rain-on-snow event, stream DOC concentrations were more uniform along the section than during the two previous sampling days, ranging between 15 and 20 mg/l ( Fig. 3C-D). The headwater lake contributed the majority of the stream water during these sampling campaigns (Fig S1). During sampling campaign C, between the snowmelt 275 and rain-on-snow event, the UCA_NOBIO scenario yielded the lowest RMSE and PBIAS, and highest R 2 (Fig. 3C). For all scenarios, NSE < 0. The second rain-on-snow sampling campaign had similar flow conditions as C, but none of the scenarios resulted in an accurate simulation (Fig. 3D, Table 2). During the post-snowmelt low flow on 2017-06-02 (Fig. 3E), step changes in stream DOC concentrations that coincided with DRIPs were well captured by the two scenarios that consider in-stream uptake (NSE=0.81 and NSE=0.5 for the "DIFF_BIO" and "UCA_BIO" models, respectively). For those scenarios, the RMSE 280 and PBIAS were also acceptable. By contrast, the scenarios without in-stream uptake ("DIFF_NOBIO" and "UCA_NOBIO") did not meet our criteria, except for R 2 .

Experimental drought 2017
The sampling campaigns in the summer of 2017 (Fig. 3F-J) consisted of a rain event, an experimental drought, a controlled lake water flood pulse and post-flood base flow conditions. The absolute gain in streamflow was small during these sampling 285 campaigns (0-5 l/s), but the relative gain varied between 9% and 90% (Fig S1). During the rain event and the drought, stream DOC concentration increased along the segment (between 17 and 28 mg/l) and distinct step changes in DOC concentrations occurred at the locations where DRIPs flow into the stream (Fig. 3F-H). For the rain event, none of the models accurately https://doi.org/10.5194/hess-2021-358 Preprint. Discussion started: 28 July 2021 c Author(s) 2021. CC BY 4.0 License. simulated the observed patterns (Fig. 3F). For the first drought sampling campaign, the model "DIFF_NOBIO" most accurately simulated the longitudinal variability in stream DOC concentration (NSE=0.48, Table 2), while the two models that accounted 290 for in-stream DOC uptake ("DIFF_BIO" and "UCA_BIO") underestimated stream DOC concentrations and were partially out of bounds (Fig. 3G). "UCA_NOBIO" overestimated stream DOC concentrations along the entire stream section ( Table 2). The second drought sampling campaign showed similar model results, but the pattern produced by scenario "DIFF_NOBIO" was less successful (Fig. 3H). The following flood pulse and post-flood base flow conditions resulted in uniform stream DOC concentrations along the stream segment. Subsequently, all models had low RMSE and PBIAS, but NSE indicated that patterns 295 were not well reproduced (Fig. 3I-J, Table 2).

Snowmelt 2018 and 2019
In the spring of 2018 and 2019, mostly snowmelt peak flow conditions were captured (Fig S1), which were characterized by flows exceeding 100 l/s and diurnal snowmelt patterns. In all sampling campaigns, stream DOC concentrations generally decreased along the stream segment, yet there were step changes indicating either dilution or enrichment of DOC at the DRIP 300 locations (Fig. 3K-O). Our modelling results show that, for the 2018 spring flood, pattern simulations met all our criteria for at least one scenario (Fig 3K-M, Table 2). The scenario that better fit the data were the "UCA_NOBIO", the "DIFF_BIO" and the "UCA_BIO", for the first, second and third sampling campaigns, respectively (Fig 3K-M, Table 2). During the 2019 spring peak flow, both "UCA" scenarios overestimated stream DOC concentrations at two DRIPs, while the "DIFF" scenarios more closely represented the observations (Fig. S1 and Fig. 3N). Further, none of the models performed well (NSE and R 2 <0.5; 305 Table 2). In the receding limb of the snowmelt peak flow, none of the models were able to capture the increments in stream DOC concentration between DRIPs (Fig. 3O).    Our observations showed that longitudinal stream DOC concentrations between the two gauging stations varied spatially during different flow conditions. While during some sampling campaigns stream DOC concentration was uniform along the 320 stream segment (around 20 mg/l), at other occasions it varied between 15 and 30 mg/l. Especially at DRIPs, we observed step changes in stream DOC concentrations, indicating that DRIPs can both decrease and increase stream DOC concentrations. We considered that these changes can be attributed to dilution and enrichment of DOC, as well as to in-stream removal of DOC.
To understand what controlled patterns of stream DOC concentrations at the various flow conditions, we simulated stream DOC concentrations based on terrestrial DOC fluxes and local in-stream DOC uptake. We compared four different model 325 scenarios with the observed stream DOC profiles. Our results showed that including DRIPs in a stream-based mixing model for DOC helps to explain longitudinal patterns through both terrestrial source-transport mechanisms as well as in-stream uptake of DOC.
Our modelling exercise demonstrated that accounting for spatial variability in lateral groundwater inputs can improve 330 simulations of stream DOC dynamics that occur along the stream segment. Most step changes in DOC coincided with DRIP locations, which suggests that supply of DOC was an important influence on stream DOC concentrations. The different patterns simulated by model scenarions DIFF_NOBIO and UCA_NOBIO shows that acknowledging spatial variability in hydrology can improve stream DOC simulations, or at least indicate where step changes occur. For example during the snowmelt event of 2018 the model scenario UCA_NOBIO was able to represent step changes in stream DOC concentrations that coincide with 335 the location of DRIPs ( Fig. 3K and Table 2). This corroborates that spatial variability in groundwater hydrology related to landscape organization has a major influence stream DOC patterns (Covino et al., 2021;Dupas et al., 2021;Rocher-Ros et al., 2019). As such, model frameworks that integrate the spatial variability in hydrology from hillslopes to stream reaches and catchments can help represent variability in hydrological connectivity between gauging stations (Jencso et al., 2009;Seibert et al., 2009). However, the other sampling campaigns show that for most of the flow conditions, the hydrology is only one of 340 the components that explains stream DOC patterns. As such our findings suggests it is equally important to represent of spatial variability in C sources in relation to landscape organization, as well as the associated biogeochemical dynamics (Grabs et al., 2012;Hale and Godsey, 2019;Wollheim et al., 2018).
We showed that accounting for in-stream DOC uptake improved predictions of longitudinal stream DOC concentrations during 345 medium flow conditions (from 10 to 50 l/s). Especially, the example in June 2017 (Fig 3E) shows how step changes are accurately reproduced by the "UCA_BIO" scenario. However, during extremely high and low flows, in-stream DOC uptake did not contribute to longitudinal stream DOC concentrations. These results concur with previous studies recently performed in headwaters (Lupon et al., 2019;Seybold and McGlynn, 2018), and suggest that aquatic biological activity becomes enhanced at the transition between low and high flows due to increases in labile DOC supply from terrestrial systems. However, our 350 results also support the idea that, beyond a certain flow threshold (> 100 l/s), the benefits of increased substrate supply become overwhelmed by low water residence times and physical disturbance (Raymond et al. 2016). Interestingly, we observed that the influence of in-stream processes on water chemistry was minimal during low flows. A possible explanation for the observed behaviour is that, during the experimental drought, stream reaches were longitudinally disconnected from each other and https://doi.org/10.5194/hess-2021-358 Preprint. Discussion started: 28 July 2021 c Author(s) 2021. CC BY 4.0 License. therefore, the role of in-stream processes to affect downstream reaches was limited (Gómez-Gener et al., 2020). Overall, our 355 results indicate that the pipe vs reactor behaviour of boreal headwater streams is closely tied to changes in lateral inputs of DOC from RZ.
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 to distinguish between 360 the spatial variability in riparian groundwater chemistry associated with different soil wetness regimes (Vidon, 2017). For example during snowmelt conditions, when occasionally DRIP water can be routed over ice (Ploum et al., 2018), local decreases in DOC concentrations in stream DOC concentrations were captured by the model. Also during the experimental drought, the representation of spatial variability in groundwater was more important that hydrology or in-stream uptake: during the experimental droughts, the scenario that approximated the stream DOC patterns the most was the "DIFF_NOBIO". This 365 result suggests that, under conditions as these, stream DOC patterns might not be directly related to the quantity of groundwater inputs, but rather to the thermal and chemical conditions that groundwater discharge creates at the local level (Briggs and Hare, 2018). However, we also see limitations of our groundwater sampling approach. For example, our groundwater sampling campaigns were 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 370 from different soil layers (Ledesma et al., 2018).
Apart from the limitations of our groundwater sampling campaigns, we identified other limitations in the hydrological components and the in-stream uptake components of the model as well. The sampling campaign of a rain event on 2017-08-02 (Fig. 3F) is a clear example of mismatch between our simulations and the observations: while stream DOC concentrations 375 increased along the reach, our model scenarios simulated a decreasing pattern. The observed step changes in stream DOC concentrations at the DRIPs is in agreement with the consistently high DOC loads that are characteristic for DRIP groundwater (Ploum et al., 2020). Also, during the rising limb of the snowmelt event (Fig. 3B), discrepancies were found between our simulations and the observed patterns: our model scenarios produced a uniform DOC pattern, while the observations showed increasing DOC concentrations along the reach which abruptly decreased at DRIPs. This observed pattern can be explained 380 by activation of organic top soils at non-DRIP reaches by rising groundwater levels (Ledesma et al., 2018) and by overland or over-ice flow at DRIPs (Ploum et al., 2018). These examples and the observed patterns that can be explained by our previous work, suggests that the model framework and the data input have some limitations representing the dynamics of streamflow generation along the stream segment. It is likely that there is discrepancy between our spatial distribution of gained streamflow based on UCA and the effective runoff to streams during different flow conditions (Ambroise, 2004;Klaus and Jackson, 2018). 385 Especially, during the spring snowmelt, our model scenarios do not capture the complexity of flow paths and the timing of lateral inputs that are associated with snowmelt hydrology (Laudon et al., 2004;Lyon et al., 2010). Besides the source-transport mechanisms on hillslope scale, other catchment characteristics affected the performance of the model as well. For example, during some high flow conditions (Fig. 3K-M, Fig S1) the lateral groundwater inputs were relatively small compared to the streamflow contributions of the upstream lake . In these cases, the observed DOC pattern was fairly 390 uniform, which made it less challenging for the different scenarios to meet our performance criteria. Also, for event-scale dynamics during rain or snowmelt events, the lake was potentially an important influence on stream conditions along this segment: while the lateral groundwater inputs respond quickly to rain or snowmelt, the lake response can be delayed .

395
For the in-stream uptake component of our model, we did not take into account processes that produce (i.e. resuspension) or remove (i.e., photodegradation, adsorption) DOC from the water column. Further, in-stream DOC uptake was assumed to https://doi.org/10.5194/hess-2021-358 Preprint. Discussion started: 28 July 2021 c Author(s) 2021. CC BY 4.0 License. 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 Vf resulted in similar model output (Fig. S2), we cannot discard the idea that Vf 400 varied among DRIPs and/or over time due to changes in groundwater DOC composition and temperature. Regardless, our model framework demonstrated that UCA can be useful to identify "reactive" reaches by considering DRIPs as major locations of biological activity. Similarly, others have demonstrated that spatial patterns in DOC concentrations along headwaters are associated with biological activity and stream water permanence driven by terrestrial flow path organization (Lupon et al., 2019, Hale andGodsey, 2019). 405 For future studies, we have identified two directions that can be useful to improve our 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 410 useful to include, especially at non-DRIPs because here groundwater levels are more dynamic compared to DRIPs (Seibert et al., 2009). Furthermore, it can be of interest to downscale the number riparian groundwater chemistry samples to understand what minimum set of groundwater samples is required to represent the spatial heterogeneity in sources of lateral DOC inputs from RZs to streams. A preliminary analysis indicated that the most optimal strategy to reduce model uncertainty was to monitor DRIPs individually, while averaging DOC concentrations at non-DRIPs (Ploum, 2021). However, given that non-415 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.

Conclusions
The objective of this study was to understand the influence of discrete riparian inflow points (DRIPs) on longitudinal patterns 420 of stream DOC concentrations along a headwater reach. More specifically, our aim was to integrate source-transport mechanisms and in-stream uptake in a model framework, by considering DRIPs as the primary driver of both hydrological and biogeochemical processes that influence patterns of stream DOC concentrations. The purpose of this modelling exercise was to disentangle the role of terrestrial DOC inputs and in-stream uptake during different flow conditions. We found that sourcetransport mechanisms as well as the in-stream uptake associated with DRIPs had occasionally a large influence on longitudinal 425 DOC patterns. The different model scenarios showed that depending on flow conditions, the dominant influence of lateral groundwater inputs and in-stream uptake on stream DOC patterns shifted. As such, the identification and characterization of groundwater chemistry of DRIPs can be used to represent major lateral groundwater inputs to stream channels, and to highlight reactive reaches within stream networks. This study contributes to the greater goal of integrating hydrological and biogeochemical models (Li et al., 2020), explicit consideration of groundwater discharges in quantitative frameworks (Briggs 430 and Hare, 2018), and the spatial assessment of removal mechanisms of DOC in stream networks (Mineau et al., 2016).