Spring coherence in dissolved organic carbon export dominates total coherence in Boreal Shield forested catchments

The wide range of forested landscapes in boreal environments store and cycle substantial amounts of carbon, although the capacity of these systems to act as either a carbon sink or source is uncertain under a changing climate. While there are clear reports of regional-scale increases in dissolved organic carbon (DOC) concentrations in streams and lakes, there remains substantial watershed-scale variability in these patterns. Coherence is a framework for examining if variables of interest within adjacent spatial units change synchronously or asynchronously through time and has been widely applied in the context of lentic hydrochemistry, and which can shed light on the relative importance of regional- vs. local-scale controls. The objective of this research was to quantify coherence in discharge, DOC concentrations, and DOC loads in forested boreal watersheds, and to what extent coherence varied by season. Coherence was assessed using data from three long-term ecological research sites spanning boreal forest environments (IISD-Experimental Lakes Area, Turkey Lakes Watershed Study, and Dorset Environmental Science Centre) that included 29 829 DOC measurements across 739 stream-years, examining correlation between stream-pairs within each site, but not between sites. Seasonal coherence in DOC export was consistent across the three sites; coherence was significantly greater in spring than all other seasons, and was strongly related to discharge coherence. Currently, the season with the greatest loads (spring) is also the most coherent season, suggesting that annual between-stream coherence may be reduced if spring becomes proportionally less important in hydrologic budgets under a changing climate. This research aids in determining which factors contribute to synchronous watershed behaviour, and which factors may contribute to the timing and extent of individual watershed-scale deviations from landscape-level patterns.


Introduction
Northern forested ecosystems contain vast stores of carbon (C) in soils, biomass, wetlands, streams and lakes (Pan et al 2011). Across the diverse types of northern forested ecosystems, the boreal zone is one of the largest terrestrial ecozones on the planet with a total worldwide area of 20 million km 2 (Schultz 2005). The storage and release of C between boreal forest systems and the atmosphere represents a substantial portion of regional, national, and global C budgets (Tranvik et al 2009). The capacity of these ecosystems to act as sinks or sources for C as greenhouse gasses (e.g. carbon dioxide and methane; CO 2 and CH 4 ) forms an important feedback loop with climatic change (Bonan 2008). Rates of C storage, burial, and loss are a function of climatic conditions (Davidson and Janssens 2006) and considerable uncertainty remains for assessing how a changing climate (warming temperatures, changing growing season length, and altered precipitation regimes) will impact C storage and cycling in boreal ecosystems (Friedlingstein et al 2014). Cycling of C in boreal regions also plays a role in maintaining suitable water quality for ecosystem health and human use (Regan et al 2017, Tugulea et al 2018 through movement and cycling of dissolved organic C (DOC).
Long-term research sites at forested catchments have collected monitoring data on C cycling and storage processes for decades (Schindler et al 1997), and recent advances in our understanding of hydrological and biogeochemical processes in heterogeneous ecosystems require a rethinking of the way we monitor transport and fate of C across interfaces. Specifically, edges between terrestrial landscapes and aquatic environments (McClain et al 2003), including the intersection of a stream across a catchment, are often spots of disproportionately rapid biogeochemical cycling. The integration across space and time of these interfaces has been identified as 'ecosystem control points' (Bernhardt et al 2017). At these locations and moments where terrestrial systems link to aquatic systems, a dominant mode of C flux is the hydrologic transport of DOC, given sufficient hydrologic connectivity (Pringle 2003).
The response of DOC concentrations in streams draining forested catchments to climatic changes and other anthropogenic pressures has been widely studied, and although there are widespread observations of DOC concentrations rising in lakes and streams across Europe and eastern North America, (Monteith et al 2007), studies have repeatedly demonstrated persistent differences across space, even across adjacent catchments with similar climatic forcing (Dillon and Molot 2005, Keller 2007, Eimers et al 2008. Two spatial scales of processes have been identified as contributing to these patterns. At the regional-scale of weather systems and climate patterns, the recovery of soils from the effects of elevated atmospheric sulphate deposition (i.e. 'acid rain') and increased rainfall and flushing of soil and wetland C have been identified as the dominant drivers of changes to DOC within receiving waters (Imtiazy et al 2020). At the watershed-scale, previous work across several environments have shown factors such as proportion of catchment as wetland (Schiff et al 1997), agricultural (Chow et al 2007, clearcut (Meyer and Tate 1983), or urban (Aitkenhead-Peterson et al 2009) control the magnitude of DOC export from catchments to streams. Often, these two scales of factors (regional factors and watershed factors) interact with each other, where individual watershed characteristics will modulate responses to the broader drivers. For example, the varying acid-sensitivity and buffering capacity of particular catchment soil properties (a watershed-scale factor) will determine the extent of DOC alteration due to changing deposition inputs (a regional-scale factor; Monteith et al 2007).
Watershed DOC flux is closely tied to the timing of hydrological and biogeochemical connectivity of different parts of the landscape at different times (Laudon et al 2011), including the timing of when ecosystem control points become activated or deactivated. The distribution of ecosystem control points may not be stationary in their relationship with climate; for instance, the dominance of snowmelt in controlling annual DOC export in seasonally snowcovered regions (Zhi et al 2019) may become of lesser importance in a future with less total snowpack accumulation and more midwinter rain-on-snow events (Casson et al 2010). In another example, the increasing severity, frequency, and intensity of wildfire events due to climatic change (Flannigan et al 2005, Barrett et al 2011 are projected to have significant impacts on forested C storage and cycling (Ingram et al 2019, Wang et al 2021 and generate new control points at temporal moments and spatial edges of burn boundaries and severities. In addition to the impacts of warming temperatures on C cycling in forested catchments (Davidson and Janssens 2006), changes to precipitation frequency, intensity, and duration have been shown to alter organic C quantity and quality in boreal headwater lakes (Warner and Saros 2019). Beyond the direct impacts of climatic change, boreal regions often experience other landscape-scale confounding factors to C cycling, such as the recovery from acidification impacting DOC mobility (Pschenyckyj et al 2020), or canopy losses driving alterations to total ultraviolet radiation at the surface and associated impacts on C cycling (France et al 2000). While the DOC concentrationdischarge relationship may be challenging to characterise at subdaily (Ducharme et al 2021) to daily (Jonsson et al 2007) scales in boreal headwater catchments, the addition of hydrologic metrics such as antecedent moisture storage and precipitation characteristics aid to explain DOC export at seasonal-and event-scales (Oswald and Branfireun 2014). Thus, given the wide range of catchment-scale and regionalscale factors, determining the extent and timing of local (single watershed) vs. regional (multiple adjacent watersheds) controls in hydrological and biogeochemical functioning can assist in determining how and when water bodies will either respond synchronously or synchronously to drivers of change.
Coherence is a framework of examining synchronous behaviour which has been widely applied in the context of lentic hydrochemistry , Crump and Hobbie 2005, Carneiro et al 2009. Understanding when, where and under what conditions hydrochemistry behaves coherently can help resolve questions of regional-vs watershed-scale controls and the dilemma of divergent patterns in proximal water bodies even under similar regional-scale forcings. Approaches to characterising coherence include pairwise correlation coefficients (r) between water bodies (Magnuson et al 1990, Järvinen et al 2002, spatially-and temporally-normalized hydrochemistry values (z-scores; Arnott et al, 2003, and other statistical approaches (Soranno  Seybold et al (2021). Absolute differences in the magnitude DOC concentrations between sites may be related to subwatershed scale factors, such as vegetation cover, soil type, and land use. Note that panels (b) and (d) maintain consistent DOC concentration coherence across seasons, whereas in panels (a) and (c) DOC is coherent in summer and winter but not in spring an autumn. Similarly, in panels (c) and (d) discharge is coherent in winter only, whereas discharge is coherent year-round in panels (a) and (b). et al 2019). In lakes, variations in DOC signatures largely reflect three main components: lateral DOC flows into the lake, the water residence time of the lake, and the amount of biogeochemical processing and transformations of the dissolved chemical constituents (Ocampo et al 2006, Flewelling et al 2012, Harms and Jones 2012. In contrast, the chemical signature of a lotic system at any given time reflects some mixture of 'new' and 'old' water, proportions of which tend to fluctuate relative to recent storm characteristics and antecedent moisture conditions, making the DOC signature reflective of the watershed at various water ages and transformations (Lu et al 2014). While coherence in the lentic context is indicative of both within-lake processes and flux of material into the lake, coherence in the lotic context is indicative of within-watershed processes and can provide insight into whether, and when, regional or local controls are most dominant. Periods of synchronous behaviour in DOC concentration across adjacent (non-nested) catchments indicate stronger regional controls across subwatershed boundaries, whereas asynchronous behaviour indicates local (within subwatershed) controls override these larger-scale climatic drivers (figure 1).
Understanding the extent and timing of coherence in C export from boreal landscapes is essential for predicting if the vast mass of C stored in boreal environments will converge to a common outcome across landscapes, or diverge to site-specific outcomes under a changing climate. Further, understanding patterns and processes which drive either synchronous or asynchronous behaviour will aid to inform future projections of changes to stream DOC and assessing our approaches to catchment-scale monitoring. As such, the objective of this research is to determine the extent of DOC export synchrony in these forested environments, and to what extent this synchrony varies by season and discharge. A secondary objective is to compare measures of stream coherence to published values of lake coherence for similar regions, and explore possible interpretations of lotic coherence in comparison with existing interpretations of lentic coherence. Our null hypothesis is that coherence in stream-pairs across boreal

Study sites and methods
This study leverages existing long term ecological research sites Turkey Lakes Watershed Study (TLW), International Institute for Sustainable Development-Experimental Lakes Area (ELA), and Dorset Environmental Science Centre (DESC), hereafter, 'sites' . At each site, several watersheds were selected (hereafter 'watersheds') in which the watershed outlet streams drain span a range of boreal forest-Canadian shield environments (figure 2). These patterned forestedwetland-lake landscapes represent a range of continental climates, with multiple decades of existing data on C cycling and transport processes (table 1). At the ELA, the three study streams are three of the inflows to Lake 239 which has been the subject of long-term study (Brunskill and Schindler 1971), with draining areas of 170.3 ha (EIF), 56.4 ha (NWIF), and 13.4 (NEIF). The watersheds surrounding the lake consist mostly of a 35 yr old fire regenerated jackpine (Pinus banksiana) forest, thin Brunisolic soils, and extensive areas of exposed granitic bedrock (Hall et al 2019). At the TLW, the nine study streams drain headwater catchments ranging 4-44 ha which have been selected as representative of the ecosystem (Jeffries et al 1988, Mengistu et al 2014) and have since been used as reference catchments for other experimental work, including the impacts of logging in adjacent impacted catchments within TLW (excluded from this study; Kreutzweiser et al 2004). Bedrock geology across the TLW is primarily Precambrian silicate greenstone overlayed by orthic ferro-humic podzols with dispersed pockets of ferric humisols (Wickware and Cowell 1983), dominated by sugar maple (Acer saccharum) growth. At the DESC, the seven study streams are six of the inflows to Harp Lake and one inflow to Plastic Lake, where both the watersheds and lake bodies have been long-term research locations

DOC and discharge data
In order to consider differences between analytical methods across sites, as well as changes to analytical methods over time at each site, details on analysis from each site are provided. Briefly, for the measurement of DOC concentrations at DESC, water samples were filtered through 80 µm polyester mesh, acidified, and flushed with nitrogen gas to remove inorganic C. Then, organic C was oxidized to CO 2 by UVR exposure in acid persulfate media and measured colorimetrically with phenolphthalein using an auto analyzer. Detailed methods for DOC measurement at DESC are provided by Hudson et al (2003) and Dillon and Molot (1997). At the ELA, DOC samples were filtered through precombusted 0.7 µm GF/F filters, acidified and digested with acid persulfate in an autoclave (1971-75), UV irradiance (1975-85), or heating to 102 • C (1986 and after). The resulting CO 2  (2003), the concentration of DOC was determined by filtering a 100 ml sub-sample through a 0.45 µm Pall GN-6 Metricel membrane filter, made from hydrophilic mixed cellulose esters, and then analysing the filtered samples within 48 h of collection on an auto analyzer using the potassium persulphate method (Crowther and Evans 1978). Although there are some slight differences between methods between sites, comparisons of coherence are only being made between stream-pairs within each region, so that variations in DOC concentration due to method are not being directly compared. At each of the catchments, discharge was recorded continuously and resampled to daily intervals, using 90 • or 120 • V-notch or H-flume sharp crested recording weirs and stage-discharge relationships from stilling walls and water-level recorders (Beaty and Lyng 1989, Beall et al 2001, Eimers et al 2009. Daily DOC loads (kg DOC ha −1 d −1 ) from each watershed were determined from daily measured discharge data and estimates of linearly interpolated DOC concentration data from approximately weekly DOC concentration measurements for each stream following Moatar and Meybeck (2005). The spring freshet dominated annual runoff and loads at each site (table 2), while a combination of baseflow conditions and biogeochemical dormancy (low stream DOC concentrations and low discharge) resulted in lower winter loads at each site. Summer and autumn had higher in-stream DOC concentrations and lower discharge than peak spring discharge, instead responding to individual storm events. The general patterns of these seasonal variations between sites are also generally consistent across streams (supplementary table 1 available online at stacks.iop.org/ERL/17/014048/mmedia) at each site. Automated baseflow separation was performed on the daily discharge data at each stream following Nathan and McMahon (1990). If over 50% of the daily discharge was attributable to stormflow, that day was categorized as a stormflow date, otherwise it was considered as baseflow date. The number and timing of days which were categorized as streamflow or baseflow was not highly sensitive to the threshold relative contribution of stormflow to daily flow beyond the 50% proportion (supplementary table 2).

Statistical analysis
Coherence was assessed as pairwise Pearson correlation coefficients (r values) between all possible stream-pairs within each site, but not between sites, following Magnuson et al (1990), for the complete data record (all years), split by season. Seasons were delineated via an astronomical definition (i.e. Winter as December-February, Spring as March-May, Summer as June-August, Autumn as September-November). Differences between pairwise coherence values between seasons were assessed via an analysis of variance test (ensuring that the homoscedasticity assumption was met; variance in coherence between seasons was not significantly different via Bartlett's test) and post-hoc t-tests between season-pairs. Differences were arbitrarily deemed significant at the p < 0.05 level. To determine the distribution of DOC anomalies (from long-term stream mean DOC concentrations) during different flow types (baseflow or stormflow), z-scores of DOC concentration for each stream were calculated using longterm mean values at each stream and rescaled to values between 0 and 1 (Arnott et al 2003, Imtiazy et al 2020. For each distribution of DOC anomalies during flow types at different sites, the median value and kurtosis of the distribution was determined.

Results
There was a significant effect of season on discharge coherence (F 3,236 = 21.5; p < 0.05). Discharge was significantly more coherent during spring than all other seasons and significantly less coherent in  winter than all other season (there was no significant difference between summer and autumn). Patterns in load coherence closely followed those of discharge coherence (r = 0.91; figures 2 and 3), where DOC loads were significantly more coherent in spring than all other seasons, although there was no significant difference in load coherence among any of the other three seasons (summer, autumn, and winter). There was not an significant effect of season on DOC concentration coherence between streampairs (F 3,212 = 2.139; p = 0.096). The mean coherence across all stream-pairs and seasons in discharge (r = 0.77) and loads (r = 0.74) was greater than the mean coherence across all stream-pairs and seasons of stream DOC concentrations (r = 0.49).
The most extreme deviations from synchrony (i.e. the lowest coherence outliers) in load and discharge between stream-pairs tended to be restricted to a set of stream-pairs at the DESC, specifically during the summer months (figure 2). These points are all values of coherence between a single watershed (PC01) and all other watersheds at DESC (figure 4; supplementary table 1).
Across all sites, temporal variations in normalised DOC concentration anomalies (deviations from a long-term stream DOC concentration mean) were related to discharge. Concentrations tended to be both higher and increasingly anomalous (more platykurtic, meaning less kurtosis and longer tailed distribution) during stormflow conditions, while concentrations tended to be lower, less anomalous (less platykurtic), and right-skewed during baseflow conditions ( figure 5, table 3). Further, the normalised DOC concentration anomalies during baseflow conditions did not vary between sites (table 3). However, the concentration-discharge relationship (represented by the difference between median DOC anomalies between baseflow and stormflow) was more pronounced at ELA and TLW while least pronounced at DESC.

Seasonal and hydrologic patterns in coherence
DOC loads and discharge were both significantly more coherent in spring than all other seasons. During spring melt, regional (climatic) characteristics, such as snowpack water equivalent and temperature, are relatively more important than local watershed-scale factors, such as soil, geology, and vegetation, in determining DOC export. Further, given the hydrologic dominance of the spring season in annual DOC export (table 2), the large springtime discharge and load coherence is responsible for contributing to overall annual DOC load coherence in boreal shield watersheds under the current freshetdriven regime. In contrast, the coherence of in-stream DOC concentrations was not significantly different between seasons (figure 2), suggesting that the impact of watershed-scale characteristics and regional (climatic) factors does not vary seasonally in determining DOC concentrations (figure 1).
Given the relationship of coherence with not only season but also discharge (figure 5), the apparent variation of coherence with season may be sensitive to the delineation of beginning and end dates of season  Table 3. Descriptive statistics across all study streams DOC anomalies by flow condition (baseflow/stormflow, following Nathan and McMahon 1990) at each study region. Normalised z-scores are the z-scores of standard deviations away from each stream's long-term mean, which were then scaled to the range 0-1. . While the definition of meteorological season (i.e. December-January-February winter) tends to capture important hydrological delineations across the boreal shield climatic setting (particularly so for a March-April-May spring melt season often coinciding with annual peak discharge at some point within that window), withinseason conditions are not homogenous, such as anomalously warm events during winter (Casson et al 2019). Further, DOC samples were not evenly distributed throughout the year, but, tended to be more frequent during hydrologically dynamic periods (typically during spring months, centered on melt events), and less frequent during baseflow conditions (typically during winter months; supplementary material 3). This sampling design is useful for accurate load estimation (Williams et al 2015), allowing for lesser extrapolation during periods which dominate flowweighted mean concentrations (Moatar and Meybeck 2005). Differences in concentration-discharge relationships between sites are indicative of different DOC transport mechanisms across boreal environments (Fork et al 2020). Specifically, variations in crossregional-scale hydrological factors, such as thin soil depths at the ELA (bare bedrock to ∼30 cm soil depth, Brunskill and Schindler 1971) and TLW (<1 m at higher elevations and 1-2 m at lower elevations; Jeffries et al 1988) relative to the DESC (50 cm-several meters; Eimers et al 2009) may explain the differential extent of the change in DOC anomalies associated with baseflow and stormflow conditions. Lower storage capacity within thin soils (ELA, TLW) leading to flashier hydrological responses may cause more of a flushing mechanism for DOC whereas thicker soils (DESC) with dampened hydrological responses may generate to a dilution mechanism, yielding reduced DOC concentrations during stormflow (figure 5).

Coefficient of kurtosis
These contrasting hydrologic settings resulted in a significant effect of site on coherence of each of stream DOC, stream discharge, and stream loads. Further examination of these patterns requires subwatershedscale understanding of hydrological connectivity and timing of the incorporation of different DOCrich or DOC-poor contributions to streamflow, often controlled by catchment-specific wetland coverage, soil distributions, or vegetation (Schiff et al 1997, Don andSchulze 2008). This present work shows that within-watershed-groups, these mechanisms (either dilution or flushing) tend to be trigged synchronously in adjacent watersheds. This pattern is most obvious during winter conditions (figure 2), which tend to be associated with baseflow conditions (increasingly right-skewed DOC anomalies; figure 5), where in-stream DOC concentrations may be most representative of soil conditions and relative lability of DOC from that watershed's parent material.

Coherence in the watershed context
This work examined the coherence framework in the lotic context between stream-pairs within three boreal shield long term ecological research sites. While this is not the first study to examine synchrony across lotic environments, either between-streams or within-streams (Carneiro et al 2009), the application of this method to a long-term record to understand seasonal variation across decades of natural variability allows for a stable and robust assessment of the coherence of C export from forested catchments.
Correlation between watershed-pairs was primarily positive across each of DOC concentration (237 out of 240 stream-pairs), discharge (all stream-pairs), and loads (all stream-pairs). This is consistent with primarily positive relationships between lake-pairs in DOC concentration coherence in previous studies (e.g. 83% of 190 lake-pairs in northern Michigan, Pace and Cole 2002). These results indicate that the coherence in DOC export from adjacent watersheds were driven by similar regional-scale processes (e.g. air temperature, rainfall, snowpack), which may be important than subwatershed-scale processes (for example, the weak relationship between mean coherence of a stream with the proportional wetland coverage in the watershed, supplementary table 4). However, these subwatershed-scale factors may remain important to determining differences in absolute magnitude of DOC concentrations or loads between watersheds (e.g. Meyer and Tate 1983, Schiff et al 1997, Chow et al 2007, Aitkenhead-Peterson et al 2009, rather than the extent of synchrony of their behaviour. Given the strong seasonality in discharge in boreal regions due to the prominence of snowmelt in the water budget, the relative seasonal impact of watershed coherence on coherence of downstream water bodies will be weighted by those season's relative loads (Magnuson et al 1990, Imtiazy et al 2020. Specifically, the highest coherence in watershed DOC loads occurred during the spring months (figure 2), corresponding to large melt events, moderate-low DOC concentration, and large loads (table 2), whereas the least coherent activity of DOC loads during autumn, winter, and summer seasons, corresponding with relative hydrological dormancy, lower DOC concentrations, and lesser loads. Given the periods which have the highest DOC loads are the most synchronous, suggests a mechanism (in addition to within-lake processes) for the similarity in lake DOC coherence values (Pace and Cole 2002) with watershed DOC load coherence. If the seasonal trend in coherence patterns is maintained while the relative contribution of spring melt decreases under a changing climate (Musselman et al 2017), it would imply a reduction in lake DOC coherence as a result. For example, if a greater proportion of precipitation and stormflow occurred during localized and intensive summer and autumnal rainfall, as observed and projected under a changing climate (Donat et al 2013, Huang et al 2018), unique subwatershed characteristics which determine DOC response could become more important to overall lake DOC coherence in a region as these summer and fall storms begin to make up a larger proportion of annual hydrologic budgets.
Absolute differences in DOC concentrations, discharge, and loads between watersheds (supplementary table 1) are masked by this analysis. Two watersheds that change synchronously relative to their own (distinct) mean values will have high coherence, but two watersheds which may maintain relatively similar values of DOC concentrations, discharge, and loads do not necessarily change synchronously (figure 1). For example, the grouping of very low discharge and load coherence values at DESC during the summer (figure 3) are generated by the lack of coherence between a single watershed (PC01) and all other watersheds at DESC (supplementary table 1). While PC01 drains into a different lake, about 30 km away (figure 4), from the other six streams it is climatically, geologically and ecologically very similar to all other DESC catchments (Casson et al 2010). The mean discharge at PC01 (0.39 mm d −1 ) in the summer is less than one standard deviation from the mean of the other six DESC streams (0.45 ± 0.08 mm d −1 ) while being hydrologically asynchronous at daily timescales from the other streams at DESC during the summer (low outlier values in figure 2). Further examination of spatial organisation of the maintenance of (or deviation from) these hierarchies can provide an additional dimension of a comparison between spatial and temporal variability at the landscape scale from days to decades (Abbott et al 2018, Dupas et al 2019. Further, differences between watershed areas may be responsible for lack of coherence within pairs of watersheds in each study region, given differences in both distributions of watershed transit times and potential time available for DOC processing (Harjung et al 2018). In the present study, watersheds only ranged in size from 4 to 191 ha. It is likely that on daily timescales, coherence would decrease between watersheds of increasing difference in size. However, these differences may or may not persist when the differences in when the differences in transit times become large relative to the sampling interval. For example, if two adjacent watersheds had transit times of 1 and 15 d, perhaps they would appear coherent on a monthly basis). Future studies examining the dynamics of coherence across watersheds differing in size by many orders of magnitude could consider lag correlation to account for dynamics that are simply offset between transit time distributions.
Coherence in several chemical parameters, including nutrients (C, N, and P) and major ions, across several lake-pairs was found to have a negative relationship with each of distance (markedly so with calcium), water residence time, and catchment:lake areal ratio, between lake-pairs at a set of lakes in the Toolik Lake area, Brooks Range, Alaska (Kling et al 2000). These lakes, however, were connected by considerable surface drainage, reflected in the relationship of coherence with distance and water-residence time reflecting the degree of actual shared water between lakes driving the extent of coherence, in addition to similar hydrological and biogeochemical processes acting on adjacent lakes and catchments. Generally, the structure of the relationship (structured, unstructured, uniform;following Webster et al 2000) between coherence and proximity appears to be limited to within-lake districts (Magnuson et al 2004), which are nested in the 'same hydrological flow system' either via surface or subsurface connectivity. In contrast, the relationship between proximity and coherence for watersheds, as in the present study (selected for which the catchments are not nested between each other) should be related to similar watershed-scale hydrological and biogeochemical processes working in adjacent catchments, rather than acting as a measure of degree of shared water along a flowpath between streams/lakes. The degree of hydrological connectivity need be considered when interpreting between-and withindistrict-level coherence to determine what is actually indicated by high or low coherence between pairs of watersheds (figure 1).

Conclusions
DOC loads across boreal shield environments in the study were significantly more coherent in spring than all other seasons, driven by strong coherence in spring discharge, in contrast to our null hypothesis that coherence would not vary across seasons or discharge conditions. This indicates that watersheds act most synchronously hydrologically in response to regional-scale drivers (in particular, in these regions, the accumulation and melting of snow). While discharge was highly coherent among streams (especially during spring snowmelt) in adjacent watersheds or sub-watersheds, local factors (e.g. watershed size/connectivity, soil types or thickness, contribution of wetlands, etc) were responsible for potentially divergent relationships (i.e. production/dilution/flushing) driving stream DOC concentrations. Currently, the most coherent season (spring) is also the season with the highest magnitude DOC loads. Under a changing climate with a decreased relative contribution of snowmelt to annual water and DOC budgets, with increased mid-winter melt events and larger summer storms, a shift to greater streamflow contributions to lakes outside of the currently most coherent season may result in reduced overall DOC coherence in lakes. Comparisons between the traditional applications of coherence approaches to lake-pairs in comparable forested environments shows that watershed-pairs are approximately equivalently coherent in DOC export as lake pairs are coherent in DOC concentration. This research has important implications for determining which factors contribute to synchronous watershed behaviour, and which factors may contribute to the timing of watershed-scale deviations from landscape-level patterns.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).