Arctic concentration–discharge relationships for dissolved organic carbon and nitrate vary with landscape and season

Climate change is intensifying the Arctic hydrologic cycle, potentially accelerating the release of carbon and nutrients from permafrost landscapes to rivers. However, there are limited riverine flow and solute data of adequate frequency and duration to test how seasonality and catchment landscape characteristics influence production and transport of carbon and nutrients in Arctic river networks. We measured high frequency hydrochemical dynamics at the outlets of three headwater catchments in Arctic Alaska over 3 years. The catchments represent common Arctic landscapes: low‐gradient tundra, low‐gradient and lake‐influenced tundra, and high‐gradient alpine tundra. Using in‐situ spectrophotometers, we measured dissolved organic carbon (DOC) and nitrate (NO3−) concentrations at 15‐min intervals through the flow seasons of 2017, 2018, and 2019. These high‐frequency data allowed us to quantify concentration–discharge (C‐Q) responses during individual storm events across the flow season. Differences in C‐Q responses among catchments indicated strong landscape and seasonal controls on lateral DOC and NO3− flux. For the two low‐gradient tundra catchments, we observed consistent DOC enrichment (transport‐limitation) and NO3− dilution (source‐limitation) during flow events. Conversely, we found consistent NO3− enrichment and DOC dilution in the high‐gradient alpine catchment. Our analysis revealed how high flow events may contribute disproportionately to downstream export in these Arctic streams. Because the duration of the flow season is expected to lengthen and the intensity of Arctic storms are expected to increase, understanding how discharge and solute concentration are coupled is crucial to understanding carbon and nutrient dynamics in rapidly changing permafrost ecosystems.

warms and thaws, vast stores of terrestrial organic and inorganic nutrients are liberated from previously frozen soil (Abbott and Jones 2015). However, complex interactions among nutrient cycles, hydrologic disturbance, and soil and permafrost conditions in high-latitude ecosystems complicates our understanding of the fate of released carbon and nutrients (Bring et al. 2016;Wrona et al. 2016). While many reach-scale studies have quantified carbon and nutrient fluxes between terrestrial and aquatic compartments of permafrost ecosystems (Mack et al. 2004;Slavik et al. 2004), understanding nutrient exchange among these components remains a significant challenge. As a result, the responses of carbon and nutrient fluxes to climate change remains a central uncertainty in predicting the permafrost climate feedback (Abbott et al. 2016;McGuire et al. 2018).
Changes in Arctic systems are often expressed in surface water networks, providing a useful mechanism to constrain carbon and nutrient balance across scales (Shogren et al. 2019(Shogren et al. , 2020Wild et al. 2019). Riverine dynamics reflect the integrated response of biogeochemical and physical processes in terrestrial and aquatic portions of an ecosystem (Creed et al. 2015).
Consequently, the analysis of lateral flux in different catchments can reveal processes that control nutrient and carbon balance (Moatar et al. 2017;Minaudo et al. 2019). Indeed, circumpolar riverine concentrations and fluxes are increasing for many solutes (Frey and McClelland 2009;Tank et al. 2016;Toohey et al. 2016), though these trends are less defined for organic carbon and nutrients due to short or low-frequency time series (Shogren et al. 2020). These solute increases have primarily been attributed to expansion of the seasonally-thawed soil and deepening of the active layer, which allows greater interaction between soils and water (Lafrenière and Lamoureux 2019). Additionally, increased disturbances, such as rainfall and wildfire, can accelerate solute delivery directly to surface water (Lafreniere and Lamoureux 2013;Rodríguez-Cardona et al. 2020).
Determining catchment-specific hydrologic and biogeochemical flux responses requires relatively high frequency data, which to date have not been widely used in Arctic riverine research. Recent technological advances now allow in situ, high-frequency measurement of biogeochemical and hydrological parameters at the event scale (Vaughan et al. 2017;Burns et al. 2019). These sub-hourly measurements of concentration (C) and discharge (Q) at a catchment outlet have greatly expanded our understanding of upstream and upslope processes and seasonality that control carbon and nutrient mobilization (Moatar et al. 2017;Minaudo et al. 2019). In-stream solute concentrations may increase, decrease, or remain constant as water flows through a river network during a storm or in response to seasonal conditions. A suite of metrics has been proposed to characterize and interpret these C-Q relationships. First, the slope of the logarithmic C-Q relationship (β) has been described as the "pulse" of the ecosystem (Godsey et al. 2009). The metric β indicates whether lateral flux is limited by material supply or hydrological transport capacity, and can be described as transport-limited (β > 0), source-limited (β < 0), or chemostatic (β = 0) (Godsey et al. 2009). Transport-limitation suggests that material is available in the ecosystem but is kept from reaching the stream network by the lateral transport capacity across the terrestrial-aquatic nexus (Zarnetske et al. 2018). Conversely, source-limitation indicates that there is insufficient supply of the material to sustain flux during moments or seasons of high hydrologic connectivity (Basu et al. 2011). Chemostatic behavior reflects flow-variant availability of material, potentially associated with weathering, microbial dynamics, and human manipulation of material sources (Godsey et al. 2009Zarnetske et al. 2018).
Because C-Q responses are typically nonlinear, additional metrics like the flushing index (FI) and hysteresis index (HI) are also used to characterize the direction and magnitude of catchment responses to flow events. The FI can be interpreted similarly to β, but specifically describes the direction of solute response during only the rising limb of the hydrograph, indicating the biogeochemical response at the onset of elevated flow (Vaughan et al. 2017;Rose et al. 2018). Then, dynamic mixing of various biogeochemical or hydrologic source contributions may result in complex or "hysteretic" patterns, with looped behavior during the rising and falling limbs of a flow event (Evans and Davies 1998). The HI reflects the dominant size, magnitude, and rotation of nonlinear responses generated during a flow event (Lloyd et al. 2016a(Lloyd et al. , 2016b. Together, the β, FI, and HI metrics can suggest the solute pool size, location, and production rate in the landscape, thus revealing the dominant mechanisms governing carbon and nutrient uptake and release at catchment scales (e.g., Rose et al. 2018).
Most of the studies noted previously have focused on temperate systems, and we are unaware of any high-resolution C-Q responses in permafrost-dominated regions. Despite the dearth of high-frequency data from Arctic rivers, using C-Q approaches to reveal integrated catchment processes in high-latitude and permafrost-dominated regions could provide critical insight into fundamental ecosystem processes and improve predictions of the responses of Arctic ecosystems to climate change (Shogren et al. 2020). For example, assessing C-Q metrics from diverse Arctic catchments could identify the drivers of largescale Arctic solute increases described previously. In addition to addressing this urgent knowledge gap, the opportunity to detect Arctic C-Q behavior may also uncover emergent catchment processes that are largely unobserved in permafrost systems. Permafrost catchments generally have lower subsurface storage, higher runoff ratios, and flashier hydrographs than better-characterized temperate systems (McNamara et al. 1998).
In this context, we collected three consecutive years of highfrequency data from three Arctic catchments on the North Slope of Alaska, USA (Fig. 1, Table 1). We used field-deployable spectrophotometers to collect dissolved organic carbon (DOC) and nitrate (NO 3 − ) concentrations every 15 min across the flow season (early June through late August or mid-September). Our first objective was to determine how distinct catchment characteristics influence emergent C-Q relationships evident at the catchment outlet (Fig. 2). Our second study objective was to explore how seasonal trends may modify C-Q responses. Informed by previous work linking catchment attributes to carbon and nutrient flux from permafrost landscapes (Harms et al. 2016;Connolly et al. 2018;Tank et al. 2020), we hypothesized that topographic slope and thaw depth, terrestrial productivity, and presence of stream-lake connections would influence catchmentscale patterns in C-Q behavior. In the lower-gradient, higherproductivity catchments (Tundra and Lake), DOC concentrations are high, and NO 3 − concentrations are low, while the opposite is true for the Alpine catchment (Harms et al. 2016;Shogren et al. 2019) (Table 1). Therefore, we predicted DOC transport-limitation (+β and +FI) and NO 3 − source-limitation (−β and −FI) in the Tundra and Lake catchments (Fig. 2a,b) and the opposite in the Alpine catchment (Fig. 2c). Furthermore, we hypothesized that the presence of lakes within the river network would partially or completely obscure event-scale dynamics reflected in the HI, because lakes interrupt land-river transport and alter material inputs due to in-lake biogeochemical processing and additional residence time (Kling et al. 2000). Therefore, increased residence time from stream-lake connections in the Lake-dominated catchment would increase variability of material sourcing, resulting in an inconsistent signal of material mobilization (±HI) (Fig. 2b). Conversely, we predicted that the Tundra and Alpine catchments, which have few lakes and ponds, would show more consistent HI values (Fig. 2a,c). We expected that abundant vegetation and soil organic matter would provide abundant sources of rapidly activated DOC pools (+HI) but more distant sources or slower activation for available NO 3 − (− HI) in the Tundra catchment (Fig. 2a). We predicted that low organic matter stocks and low nutrient demand would create the opposite pattern in the Alpine catchment, with hydrologically activated DOC moving into the channel on the falling limb of the hydrograph (−HI) and NO 3 − rapidly flushing from a proximate source (+HI) (Fig. 2b).
We also hypothesized that β, FI, and HI would change for each catchment in response to seasonal differences associated with thaw depth and nutrient demand. As noted above, soil and permafrost conditions exert a strong control on Arctic surface water chemistry. Therefore, hydrochemistry must be viewed in light of the unique natural seasonal influences of a thickening active layer and evolving subsurface flowpaths (Harms and Jones 2012;O'Connor et al. 2019;Chiasson-Poirier et al. 2020), changing terrestrial and aquatic productivity and nutrient demand (Treat et al. 2016), and spatial lake dynamics (Kling et al. 2000). For example, in the higher productivity landscapes (Tundra and Lake-dominated) we hypothesized that increasing organic matter availability and nutrient demand from June through late August would result in a seasonal progression of greater DOC "flushing" ("β, "FI) and NO 3 − sequestration (#β, #FI) signals ( Fig. 2a,b). Relatedly, in these catchments through the thaw season, we predicted that HI would indicate a progressively increasingly rapid source activation and mobilization for DOC ("HI) and slower activation for NO 3 − (#HI) signals ( Fig. 2a,b). In the Alpine catchment, we hypothesized declining DOC availability (#β, #FI, #HI) and greater NO 3 − availability ("β, "FI, "HI) signals through the thaw season. Additionally, we posited that C-Q signals in catchments with lake influence would be "buffered" relative to catchments with minimal stream-lake connections, altering emergent seasonal trends in C-Q behavior observed at the river outlet. In landscapes with high lake influence, stream flow may be repeatedly interrupted by increased residence within a mixed lake "reactor" (Sadro et al. 2012). Alternatively, seasonal thermal stratification could effectively reduce the *Here, Q 50 is defined as median flow; 50% of the time flows have been measured at or above that level and 50% of the time flows have been at or below that level (Bond 2019). † Runoff ratios were calculated by dividing the volume of event discharge (L event −1 ) by area-estimated volume of precipitation (L event −1 ).
residence time of the lake by constricting inflow and outflow to the epilimnion (Cortés and MacIntyre 2020). We therefore expected that both spatial and seasonal lake dynamics would obscure terrestrial changes associated with thaw depth and plant demand as a result of lake uptake and removal of carbon and nutrients, resulting in no net seasonal change of HI (Fig. 2a,c).

Study catchments
Our three study catchments are associated with the Arctic Long-Term Ecological Research site near the Toolik Field Station on the North Slope of Alaska. The sites are all within the continuous permafrost region, where more than 90% of the area is underlain by perennially frozen ground (Hobbie and Kling 2014). Each catchment represents a common ecosystem type in the Arctic: low-gradient tundra (the Kuparuk River, hereafter Tundra), lake-dominated tundra (Oksrukuyik Creek, Lake), and highgradient alpine tundra (Trevor Creek, Alpine; Fig. 1, Table 1). The Tundra and Lake rivers drain wet acidic tundra, where the dominant vegetation consists of tussock tundra, shrubs, wet sedge, heath, and lichens. Mean residence times of small lakes within Lake-dominated catchments range from several weeks (< 1 month) to several years (> 5 years) (E. Haines pers. comm.). The Alpine catchment is a steep alpine stream network, primarily draining rocky hillslopes with narrow riparian areas vegetated by heath and willow.
We quantified a variety of catchment characteristics (Toolik Environmental Data Center Team, 2019), including topographic, stream network, and biological metrics (Table 1). Of note, we use Normalized Difference Vegetation Index (NDVI) as a proxy for terrestrial primary productivity (Boelman et al. 2003). We derived NDVI data from imagery acquired in summer 2012 by the ETM+ sensor on Landsat 7 (courtesy of the U.S. Geological Survey).

Fig. 2.
Conceptual diagram depicting predictions for our three study catchments, (a) Tundra, (b) Lake, (c) Alpine of (i) DOC and (ii) NO 3 − C-Q β, HI, and FI. If the prediction text is bolded and black, it signifies that our predictions matched our observations. If the prediction text is not in black text, it signifies that our prediction did not match our observation or was unable to be assessed.

In-situ data collection
We collected high-frequency estimates of stream flow, DOC, and NO 3 − from the three catchments for three consecutive years (2017-2019; Table 1). We monitored each catchment for the majority of the flow season from early June through early September, though these efforts did not always capture the "shoulder" seasons (Shogren et al. 2020). In each, we estimated discharge using co-located, atmospherically compensated pressure transducers (Onset HOBO, Bourne Massachusetts, USA) that recorded water depth (stage) at 10-min intervals. We converted these stage data to discharge using a discharge-stage rating curve and velocity-area calculations from weekly field measurements (Turnipseed and Sauer 2010). Concurrently, we measured water chemistry every 15-mins using submersible UV-visible spectrophotometers (s::scan Messtechnik GmbH, Vienna, Austria), which were co-located with the water level sensors. The spectrophotometers measure light absorbance at wavelengths from 220 through 750 nm, and spectra are corrected for optical path length (5 mm at the Alpine site to account for higher turbidity, and 35 mm for the Tundra and Lake sites where there was less turbidity and lower NO 3 − ) ( Table 1). In the minute prior to collecting the absorbance reading, the spectrophotometers automatically cleaned their lens with a rotating brush. We housed the sensors in protective PVC tubing anchored on the streambed parallel to flow. The deployment dates varied slightly each year, depending on weather and flow conditions (Fig. 3, Table 1).
We collected biweekly grab samples at the sensor sites to analyze DOC and NO 3 − concentrations in the laboratory to allow site-specific calibration of the absorbance spectra (Ruhala and Zarnetske 2017). Each field sample was filtered to pre-ashed 0.7 μm with a glass fiber filter (GF/F Whatman) into a clean HDPE bottle. The filtered water samples were acidified and refrigerated (DOC) or frozen (NO 3 − ) until laboratory analysis. We measured DOC on a Shimadzu TOC-L analyzer using a combustion catalytic method and NO 3 − on a QuickChem Lachat analyzer (2017) or a SEAL AA3 segmented flow analyzer (2018-2019). With the lab-measured concentrations, we used a partial least squares regression variable selection approach to generate robust calibration relationships between spectra and observed grab samples (Etheridge et al. 2014) (SI Fig. 4, SI Table 3). We used the packages pls (Mevik and Wehrens 2007) tto fit the calibration in R (R Core Team 2014).

Annual C-Q analysis
To compare flow seasons within and among our study catchments, we normalized by catchment area and plotted continuous discharge measurements as specific discharge (L km −2 s −1 ) (Fig. 3). Using the baseflow function in the R hydrostats package (Bond 2019), we calculated several annual flow statistics: mean daily discharge, median daily discharge (Q 50 ), mean baseflow, and average coefficient of flow variation (Table 1). For each year, we estimated the overall C-Q relationship across all flow conditions for DOC and NO 3 − (Table 1), fitting the relationships as C = aQ β , where the exponent β represents the slope between log-concentration and log-discharge (Godsey et al. 2009).

Event-scale C-Q analysis
For the event-scale analysis (e.g., SI Fig. 1), we defined high flow events as any hydrologic response to rainfall or snowmelt which resulted in a 10% increase in discharge. We estimated baseflow and stormflow using the hydrostats R package, which measures the central tendency of baseflow using a Lyne-Hollick baseflow filter (Bond 2019). We classified prolonged events with multiple peaks as embedded, separate events so that we could quantify biogeochemical response to each rise in discharge. In total, across the Tundra, Lake, and Alpine catchments, respectively, we completed our analysis and captured 13, 8, and 6 events (2017); 9, 8, and 7 events (2018); and 9, 6, and 11 events (2019; Fig. 3). For the Tundra and Lake catchments, where we had nearby precipitation data (from the Toolik Environmental Data Center, https://toolik.alaska.edu/ edc/), we estimated runoff ratios for each event using hourly precipitation data from the nearest weather monitoring station as the ratio between total discharge (as L km −2 ) and total precipitation (L km −2 ).
We calculated a series of C-Q metrics for each high flow event. Using the C = aQ β equation noted above, we estimated a log-linear relationship (β) between concentration and discharge. We interpreted β to infer transport limitation (β > 0), source limitation (β < 0), or chemostasis (β = 0) (Godsey et al. 2009Zarnetske et al. 2018). As C-Q β does not always capture the variability around the slope estimate inherent to these responses, the direction of β is not a complete characterization. Therefore, we also calculated FI to determine whether the relationship between concentration and discharge from the onset of flow to peakflow is generally increasing, decreasing or flat. The values of FI range from − 1 to + 1, where negative values signify declining concentration during the rising limb of an event, while positive values indicate increasing concentration or flushing (Vaughan et al. 2017).
We also calculated HI using the range of normalized C and Q values between the rising and falling limb at every 5% percent increase in discharge (Lloyd et al. 2016a(Lloyd et al. , 2016b. The HI normalizes the C-Q responses so that values range between − 1 and + 1. Negative values (counter-clockwise hysteresis) suggest that high-concentration sources are distant from the sampling location or that early contributions with low concentration are mixed with higher-concentration sources later in the event (Evans and Davies 1998). Positive values (clockwise hysteresis) suggest a limited solute supply that is likely close to the measurement location or readily hydrologically accessed (Evans and Davies 1998).

Statistical analysis
We performed several statistical tests to identify catchment-specific and seasonal trends in C-Q metrics within and among each catchment. First, we explored within-site interannual differences in C-Q β, HI, and FI across years using one-way ANOVA. Second, we compared cross-site differences in runoff ratios, β, FI, and HI across catchments and sampling years using two-way ANOVA (model structure: CQ metric Catchment + Year + Catchment * Year). For all ANOVA tests, we determined differences using Tukey's post hoc tests when significant. Then, to investigate seasonal or intra-annual trends in C-Q metrics within each catchment, we also reported simple linear regressions for each solute. Because  Seasonal trends in estimated runoff ratios from individual high flow events from for the Tundra (gray points) and Lake (white points) study catchments in 2017 (circles), 2018 (triangles), and 2019 (squares). Solid black lines represent significant regression between runoff ratio over time at p < 0.05. (c) Boxplots representing the annual interquartile range and median of estimated runoff ratio in the Tundra and Lake catchments from 2017, 2018, and 2019. We did not have reliable precipitation data to estimate runoff ratios in the Alpine catchment, so these data are not included. sequential flow events are not independent due to antecedent moisture and biogeochemical conditions, we performed Mann-Kendell tests to confirm significance of seasonal trends. Finally, because we were interested in the overall hydrologic response of each catchment, we used a Permutational ANOVA (PermaANOVA) using Distance Metrics in the adonis function in R (Oksanen et al. 2019) to detect catchment-specific similarities/dissimilarities in overall C-Q responses. To determine differences in overall biogeochemical behavior, we included all C-Q responses (β, HI, and FI) for both DOC and NO 3 − in the PermaANOVA analysis. Lastly, to visualize the difference in DOC and NO 3 − C-Q behavior, we plotted our results in a Non-Dimensional Metric Scaling (NMDS) framework using the R package vegan (Oksanen et al. 2019).

Results
Catchment flow regimes We present the hydrographs and precipitation for all monitoring seasons in Fig. 3, and a suite of flow statistics in Table 1. Except for the large spring snowmelt runoff event, most flow in the Tundra and Lake catchments was derived from rainfall (Fig. 3a,b). Unlike the Tundra and Lake catchments, flow in the Alpine catchment was strongly influenced by glacier and snowfield meltwater, which created diel patterns in discharge (Fig. 3c). Runoff ratios in both the Tundra and Lake catchments tended to decrease across the flow season (R 2 > 0.5 and p < 0.05 for both, Fig. 4a,b). Within each catchment, runoff ratios were not significantly different across monitoring years (Two-way ANOVA, F 2,47 = 0.914, p = 0.70,   . 4c). However, runoff ratios were significantly higher in the Tundra than in the Lake catchment when comparing means (Two-way ANOVA, F 1,47 = 13.56, p < 0.001, Fig. 4c). Furthermore, there was no significant interaction between catchment and year (Two-way ANOVA, F 2,47 = 0.20, p = 0.82). We were unable to estimate the runoff ratio in the Alpine catchment due to a lack of precipitation data in that area. The hydrological responses in the Lake catchment were muted compared to the Tundra catchment, with lower specific discharge for similar events (Fig. 3b).

Annual C-Q patterns
The annual, composited C-Q analysis revealed a range of hydrochemical patterns within and among catchments (Table 1). Composite C-Q relationships for DOC and NO 3 − showed a range of patterns that varied by constituent and among catchments (SI Figs. 2 and 3). For DOC, the two lowgradient tundra catchments (Tundra and Lake) showed a strong positive C-Q trend across flow conditions (β > 0), indicating a general enrichment effect of concentration as discharge increased (SI Fig. 2a,b). These patterns were consistent among the three study years. In contrast, the composite patterns for DOC in the Alpine catchment were non-significant, implying a chemostatic relationship (β = 0) (SI Fig. 2c).
However, this lack of a significant relationship could be associated with daily variation in discharge driven by diel snowmelt (Fig. 3). For NO 3 − , we generally found an overall dilution signal in the Tundra and Lake catchments (SI Fig. 3a,b), and a non-significant C-Q relationship in the Alpine catchment (SI Fig. 3c).

Event-scale C-Q metrics
We used event-scale C-Q metrics to explore landscape specific trends in hydrologic mobilization and sourcing for DOC and NO 3 − (see SI Fig. 1, e.g., C-Q plots). Our analysis revealed relationships that were not discernable in the annual C-Q analysis (SI Figs. 2 and 3), showing high variability in C-Q metrics among and within catchments for DOC (Fig. 5a) and NO 3 − (Fig. 5b) (detailed ranges are reported in SI Table 1).
First, we compared within-site β, FI, and HI and found no statistical difference across years in any of our study catchments and for either solutes (One-way ANOVA, p > 0.05 for all tests). Then, we compared among-site differences in C-Q metrics for all catchments (Tundra, Alpine, Lake) and monitoring years (2017, 2018, or 2019). Across most of these models we found no significant interaction terms between catchment and year (Two-way ANOVA, p > 0.05 for most tests). The only exception was a significant interaction term for DOC β. The (ii) FI, and (iii) HI. We note any significant seasonal trend (p < 0.05) with a solid black line, with R 2 and Mann Kendall (τ) statistics reported for each significant relationship. Across all plots, the horizontal reference line at "0" is chemostatic behavior. For C-Q slope (β) (row i), when β is > 0 the metric indicates transport limitation, β < 0 indicates source limitation or dilution signal, and β = 0 indicates chemostasis. For C-Q FI (row ii), FI > 0 implies rapid initial flush in the hydrograph, FI < 0 indicates an initial dilution, and FI = 0 indicating chemostasis. Finally, for C-Q HI (row iii), HI > 0 results from a clockwise hysteresis loop, HI < 0 from a counter-clockwise loop, and HI = 0 as non-hysteretic behavior. main effect of catchment on DOC β was significant (Two-way ANOVA, F 2,73 = 85.27, p < 0.001), while the main effect of year on DOC β was also significant (F 2,73 = 4.19, p = 0.02). These main effects were qualified by a significant interaction between catchment and year (F 2,73 = 2.58, p = 0.014), indicating that values for DOC β were dependent on sampling year. Values for DOC β tended to be similar and generally positive between the Tundra (mean β = 0.36, 0.39, and 0.33, for 2017-2019 respectively) and lower in the Lake catchments (mean β = 0.07, 0.29, 0.48) but tended to be lower in the Alpine catchment (mean β = − 0.68, − 0.62, and − 0.34). As with DOC β, we found significant differences in DOC FI by catchment (Two-way ANOVA, F 2,73 = 10.21, p < 0.001) but not year (F 2,73 = 2.50, p = 0.09). DOC FI during flow events was similar between the Tundra (mean FI = 0.40, 0.53, and 0.62) and Lake catchments (mean FI = − 0.11, 0.57, and 0.51; Tukey's HSD: p = 0.23) but lower in the Alpine catchment (mean FI = 0.07, − 0.09, and − 0.06; Tukey's HSD: p = 0.007). In contrast to β and FI, DOC HI values were highly variable, ranging from positive to negative within each catchment. HI was not significantly different among catchments (Two-way ANOVA, F 2,72 = 0.62, p = 0.54) or years (F 2,72 = 1.46, p = 0.24).
Our estimated event-scale C-Q metrics for NO 3 − were highly variable among the studied landscapes (Fig. 5b) Exploring seasonal patterns of event-scale C-Q responses After establishing differences in mean C-Q behavior, we used linear models to explore potential seasonal trends in event-scale β, HI, and FI for DOC (Fig. 6) and NO 3 − ( Fig. 7; statistics reported in SI Tables 2 and 3). The main purpose of the analysis was to test for seasonal effects on observed hydrologic − C-Q metric across the season for (a) Tundra, (b) Lake, and (c) Alpine for (i) C-Q slope (β), (ii) FI, and (iii) HI. We note any significant seasonal trend (p < 0.05) with a solid black line, with R 2 and Mann Kendall (τ) statistics reported for each significant relationship. Across all plots, the horizontal reference line at "0" is chemostatic behavior. For C-Q slope (β) (row i), when β is > 0 the metric indicates transport limitation, β < 0 indicates source limitation or dilution signal, and β = 0 indicates chemostasis. For C-Q FI (row ii), FI > 0 implies rapid initial flush in the hydrograph, FI < 0 indicates an initial dilution, and FI = 0 indicating chemostasis. Finally, for C-Q HI (row iii), HI > 0 results from a clockwise hysteresis loop, HI < 0 from a counter-clockwise loop, and HI = 0 as non-hysteretic behavior.
patterns for DOC and NO 3 − across monitoring years. We note that in some monitoring years, the high variability in C-Q responses made either increasing or decreasing trends non-significant. For DOC, we found only a few significant seasonal trends within our studied landscapes (Fig. 6a-c). In 2017, we found a significant increase in DOC β in the Tundra catchment (R 2 = 0.31, p = 0.04; Fig. 6 row i). In 2018, there was no seasonal change in DOC β in any of our study catchments. In 2019, there was a seasonal decline in β within the Tundra catchment (R 2 = 0.44, p = 0.04), though the values of β were always > 0. Generally, significant trends for DOC FI followed those of β (Fig. 5 row ii, SI Table 2). We only observed a significant seasonal trend for DOC HI in 2018 within the Tundra Fig. 8. Biplot of (a) DOC and (b) NO 3 − for (i) β × FI, (ii) β × HI, and (iii) HI × FI for the Alpine (black fill), Lake (white fill), and Tundra (gray fill) catchments for all monitoring years. Across all plots, the horizontal and vertical reference line at "0" notes chemostatic (β and FI = 0) or non-hysteretic (HI = 0) behavior. In this figure, CW and CCW indicates clockwise and counter-clockwise rotation, respectively. catchment (Fig. 5, row iii), with significantly increasing HI values as the season progressed (R 2 = 0.63, p = 0.01).

Comparing C-Q responses across different landscape types
To compare event responses across our studied landscapes, we created biplots of β, HI, and FI for DOC (Fig. 8a) and NO 3 − (Fig. 8b). These plots highlight the principal relationships between the overall direction and behavior of each storm response and allow generalizations of C-Q dynamics across the three distinct landscape types. Across the two low-gradient catchments (Tundra and Lake), C-Q responses for DOC largely clustered in the upper quadrants suggesting transport-limited hydrodynamic behavior (β > 0) resulting from a rapid initial flush of material (FI > 0; Fig. 8a, row i), though the direction of the resulting C-Q hysteresis loops were highly variable (HI ranging between − 1 and 1; Fig. 8a, row ii). Conversely, C-Q responses for DOC in the high-gradient Alpine catchment largely clustered as a dilution signal (β < 0) that varied between an initial dilution or flush response (Fig. 8a, row i). Like the variable hysteresis observed in the Tundra and Lake catchments, DOC HI values in the Alpine catchment varied from negative to positive (Fig. 8a, row ii).
C-Q dynamics for NO 3 − were distinct between the Tundra/ Lake and Alpine groupings, but opposite of DOC (Fig. 8b). Tundra and Lake NO 3 − C-Q metrics largely clustered in quadrants associated with diluting or source-limiting hydrodynamic behavior (β < 0) with a rapid initial dilution (−FI), but generally clockwise hysteresis (+HI, Fig. 8b). In contrast, NO 3 − C-Q responses for the high-gradient Alpine catchment were predominantly an enrichment signal (β > 0) with a rapid initial flush response (FI > 1) (Fig. 8b, row i). Unlike the clockwise hysteresis observed in the Tundra and Lake catchments, NO 3 − largely clustered in the counter-clockwise quadrants.
Lastly, we explored the compiled C-Q relationships using a PermANOVA approach. Comparing the distribution of eventscale responses for β, HI, and FI combined, we determined that C-Q behaviors were statistically unique among our study catchments (F 2,71 = 15.93, R 2 = 0.39, p = 0.001) and plotted the NMDS (Fig. 9, model stress = 0.12). This analysis substantiates the overall finding that the biogeochemical responses of the Tundra and Lake catchments are similar and distinct from the Alpine catchment.

Discussion
In this study, we used high-frequency Q, DOC, and NO 3 − data to assess carbon and nutrient dynamics across three flow seasons in distinct Arctic headwater landscapes. The eventbased C-Q analyses suggested that landscape characteristics, including slope, vegetation type, and stream-lake connectivity, influence the hydrologic mobilization and transport of biogeochemical solutes from headwater catchments underlain by permafrost. Our results also indicated that within-season variability influences the activation, mobilization, and transport of DOC and NO 3 − conveyed during flow events. Together, these results underscore that accurately representing Arctic riverine fluxes must account for the variation in hydrodynamic responses across landscapes and the thaw season (Shogren et al. 2020).

C-Q responses revealed high variability in solute generation
Characteristics of C-Q relationships can indicate terrestrial nutrient source availability or limitation (Minaudo et al. 2019). We therefore expected Arctic C-Q metrics to signal rapid hydrologic export of DOC and reduced NO 3 − release in catchments draining low-gradient, and high-productivity landscapes (see predictions in Fig. 2a,b). Our Tundra and Lake study catchments drain low-slope wet acidic tundra that supports productive vegetation growth by Arctic standards, where we might expect longer or more complex flowpaths through organic-rich soil layers (Neilson et al. 2018) and elevated terrestrial N demand and aquatic uptake (Townsend-Small et al. 2011;Liu et al. 2018). In these conditions, DOC is readily released from terrestrial source areas as hydrologic connectivity increases during precipitation events (Neilson et al. 2018). As water activates biogeochemical source areas along these flowpaths, DOC is rapidly transported to the stream network (Ågren et al. 2007). Conversely, N-limitation is common in tundra soils, where high rates of uptake and assimilation by vegetation and soil microorganisms limits the transport of inorganic N species (Keuper et al. 2012). Furthermore, lower topographic relief allows for the development of riparian zones and small wetland areas that promote the removal of NO 3 − by denitrification (Harms et al. 2016) and constitute a potential source of dissolved organic matter into streams (Fouché et al. 2017;Connolly et al. 2018). Accordingly, in the Tundra and Lake catchments, DOC stores were readily activated and transported to the stream network ( Fig. 2a,b, +β), while NO 3 − release was often sourcelimited as a result of strong terrestrial N-limitation (Myrstener et al. 2018) (Fig. 2a,b, −β). Our results support previous studies showing longer water residence times within a lower gradient landscape correspond with elevated organic solutes (O'Donnell et al. 2010;Harms et al. 2016) and decreased concentrations of inorganic N (O'Donnell and Jones 2006). Furthermore, the low-gradient catchments (Tundra and Lake) exhibited minimal overlap with the hydrodynamic behaviors observed in the Alpine catchment (Fig. 8), suggesting fundamentally different processes regulating water and solute transport in Arctic mountains vs. foothills in this region. The steeper Alpine catchment C-Q behavior reflected lower terrestrial productivity and shorter, faster hillslope flowpaths resulting in less DOC available for hydrologic transport (Harms et al. 2016;Connolly et al. 2018;Treat et al. 2018). Arctic alpine catchments are exposed to higher rates of wet and dry N-deposition (Choudhary et al. 2016), have shorter and faster flowpaths with fewer opportunities for biotic NO 3 − removal (Harms and Jones 2012;Harms et al. 2019), and generally experience lower terrestrial N-demand. Subsequently, in the Alpine catchment, we observed persistent NO 3 − transportlimitation and enrichment (Fig. 2c, + β), in addition to the overall pattern of DOC dilution during changes in discharge (Fig. 2c, −β). Our analyses of β indicated catchment-dependent mobilization and transport for DOC and NO 3 − . Consequently, we used the HI and FI to provide additional information about source areas and catchment mobilization processes (Fig. 8ii,iii). The HI and FI metrics have been used to gain insights into complex temperate catchments, which characteristically have multiple hydrologic storage zones and complex mixing dynamics (Harvey and Gooseff 2015). Specifically, these metrics can describe the timing and indicate the proximity of materials delivered to the stream channel (Vaughan et al. 2017;Rose et al. 2018;Burns et al. 2019). For example, a recent study found DOC was rapid flushed (+FI) but in a counter-clockwise hysteresis pattern (−HI) from forested catchments in Vermont, showing input of available DOC from a stream-proximate source during the rising limb of a storm, such as the riparian zones, followed by even higher DOC from newly activated source as hydrologic pathways connect across the catchment (Vaughan et al. 2017). In the same study, there was clear dilution (−FI) and clockwise hysteresis (+HI) for NO 3 − during the same flow events, indicating the mobilization of proximate, high concentrations at the onset of flow, followed by lower concentrations after peak flow (Vaughan et al. 2017).
We observed a surprising degree of variability in FI and HI values across our permafrost underlain landscapes. The variability in FI and HI likely indicates the presence of multiple source areas (McGlynn and McDonnell 2003), terrestrial and aquatic removal mechanisms (Kendrick et al. 2018;Harms et al. 2019), and seasonally evolving hydrologic pathways (Lafrenière and Lamoureux 2019; Chiasson-Poirier et al. 2020). Additionally, a variable response may reflect antecedent conditions altering landscape nutrient and carbon pools (Ågren et al. 2007) or seasonal changes in nutrient demand and availability (Treat et al. 2016). However, there was some similarity in DOC and NO 3 − mobilization when we considered both flushing and hysteresis behavior. C-Q patterns in the Tundra and Lake catchment tended to rapidly flush DOC (+FI), while the range of hysteresis (±HI) suggests two major carbon sources: (1) proximate, high-DOC sources that are either depleted or mixed with lower concentrations after peakflow (+HI), or (2) hydrologic activation of larger, more distant sources later during the storm (−HI). Both of these pathways have been observed in many Boreal rivers, which often have large organic matter stores near the stream channel (Ledesma et al. 2015). Conversely, FI and HI patterns in the Alpine catchment suggest an overall depletion of carbon (+HI, −FI), consistent with our predictions (Fig. 2c). NO 3 − HI was similarly variable across our study catchments, with C-Q patterns generally indicating proximate N sources that are rapidly depleted (−FI, +HI, as in Tundra), variable sources contributing to N export (+FI and ±HI, as in Lake), or high NO 3 − from proximate sources that are mixed with low-concentration water later in the event (+Fi, −HI, as in Alpine). Still, the diversity of responses in the Alpine catchment is attributable to the dominance of spatiotemporally complex glacial and snowfield melt typical of many high-alpine catchments, which can decouple river flow from overall catchment connectivity (i.e., river discharge can increase without a substantial rise in connectivity due to activation of meltwater) (Mcknight et al. 2015). While our study cannot directly pinpoint transport mechanisms or source pools of solutes (Burns et al. 2019), our results underscore that Arctic catchment processes vary by landscape type, which regulates carbon and nutrient transport from land to surface water networks (Tank et al. 2020). Identifying the source pools of landscape-derived riverine DOC and NO 3 − in high-latitude rivers will require additional analysis using combinations of methods and measurements to characterize fluxes from actively thawed portion of hillslopes (Harms and Jones 2012;Lafrenière and Lamoureux 2019).
Trends in C-Q β suggest dominance of landscape characteristics over seasonal drivers of DOC and NO 3 − mobilization Water chemistry monitoring at the catchment outlets is an integrative measure of the upstream aquatic and terrestrial processes that result in biogeochemical production, removal, or transformations as material is transported downstream (Creed et al. 2015). Therefore, within a given flow season, we expected that increasing thaw depth and changes in terrestrial and aquatic nutrient demand would alter the C-Q dynamics for both DOC and NO 3 − (Harms and Jones 2012; Chiasson-Poirier et al. 2020) (Fig. 2). Indeed, we anticipated consistent seasonal trends in hydrologic responses, at least for the Tundra and Lake catchments (Fig. 2a,c). Previous studies in lowgradient catchments have found that deeper flowpaths result in declining DOC export (Striegl et al. 2005;McClelland et al. 2007) and elevated NO 3 − (Harms and Jones 2012;Khosh et al. 2017). Though we found several instances of seasonal shifts in DOC and NO 3 − β in the Tundra and Alpine catchments, we did not observe consistent seasonal trends. In fact, our findings of late-season transport-limitation of DOC are in contrast with these previous studies, which may be due to our focus in this study on storm events, while previous studies rarely include high flow events (Shogren et al. 2020). Hence, relative to these previous studies, our results may more directly reflect the contribution of DOC from organic layers that are only activated during high flows (Ledesma et al. 2015) or the result of increased availability of carbon derived from terrestrial vegetation later in the season (Treat et al. 2018) (Fig. 2a). Our finding of NO 3 − dilution behavior during most flow events in the Tundra catchment emphasizes the relatively low N-availability and high N-demand typical of these Alaskan tundra landscapes (Fig. 2a), despite increasing NO 3 − concentrations later in the thaw season (Treat et al. 2016;Khosh et al. 2017). In contrast, in the Alpine catchment specifically, the consistent pattern of negative DOC β across all sampling years (Fig. 6a) can likely be attributed to lower ecosystem productivity generating organic carbon and, subsequently, storing organic carbon in the subsurface of this landscape (Keuper et al. 2012;Harms et al. 2016 Rather, we interpret the lack of observed seasonal dynamics to mean that seasonally-variable factors such as thaw depth and internal lake processes do not modulate hydrochemical responses as effectively as landscape characteristics (e.g., slope, vegetation, and stream-lake interactions) which are the predominant C-Q controls at the catchment outlets (Laudon and Sponseller 2018). Landscape dependence has been observed in other permafrost regions regarding the hydrochemical response to wildfire, thermokarst, and soil hydrology, further highlighting the importance of catchment context on biogeochemical exports (Laudon et al. 2017;Tank et al. 2020). Emergent C-Q properties are often used to investigate terrestrial to aquatic ecosystem linkages. However, substantial additional material processing can take place within the aquatic ecosystems when a river network is hydrologically connected to a series of lakes, as water and reactive solutes sequentially flows through successive stream-lake environments (Kling et al. 2000;Sadro et al. 2012). Stream-lake connectivity is likely a factor confounding our ability to detect emergent within-season patterns and intra-annual C-Q responses within the Lake-dominated catchment specifically (Shogren et al. 2019). While we do not have diagnostic information about the variable environmental conditions and biogeochemical reactions along consecutive stream to lake interactions in our Lake-dominated catchment, it is known that stream lake interactions can alter the terrestrial and stream signal of solute delivery measured at the catchment outlet. Lakes are a key factor due to internal processing and subsequent release of either attenuated or elevated solute concentrations to the river network (Kalinin et al. 2016). Lakes may also integrate heterogeneously distributed sources, leading to chemodynamic export regimes (Godsey et al. 2009). Additionally, the presence or absence of thermal stratification in Arctic lakes of the study region are known to alter water routing through stream-lake connections. In many small Arctic lakes, the period immediately after ice-off is often characterized by the short presence of a strong thermocline, where stream water flow may preferentially travel above the colder, nutrient-rich water layer (Cortés and MacIntyre 2020). Then, later in the season, lake turnover from solar convection and turbulence from wind rapidly allows stream water to mix into the lake depth profile (MacIntyre et al. 2006). As a result, biogeochemically-reactive material transport in the Lake catchment was likely dependent on a combination of internal lake reaction dynamics, thermal stratification or lake mixing regimes, and lake water residence times and network position (Kling et al. 2000;Sadro et al. 2012).
While internal lake processing via biotic assimilation or denitrification can reduce inorganic nitrogen concentrations, such as NO 3 − in the Lake catchment (Kling et al. 2000), warmer temperatures associated with lake outflows may simultaneously promote nitrification and therefore increase NO 3 − . Therefore, under some conditions, NO 3 − may be internally produced within the lakes of the Lake catchment such that there is ample inorganic N available for hydrologic transport (β> 0). During other conditions, perhaps as N-demand increases during the seasonal progression, lake NO 3 − may become source-limited (β < 0), such that precipitation events result in a dilution behavior. Overall, the C-Q observations from the Lake catchment underscore that the presence of stream-lake connectivity is likely a major controlling characteristic for biogeochemical fluxes from Arctic headwater streams when lakes are present. The potential influence of stream-lake interactions is corroborated by spatiallydistributed sampling from the same catchment network (Shogren et al. 2019). In that study, DOC and NO 3 − concentrations were spatially redistributed across the Lake catchment's network throughout the flow season, while the Alpine and Tundra catchments showed spatial stability. When paired with spatially-explicit understanding of biogeochemical contributions (Shogren et al. 2019), our high-frequency hydrochemical analyses revealed that the influence of streamlake interactions is propagated downstream such that nutrient and carbon C-Q responses are altered at the catchment-scale, resulting in different solute flux characteristics than other landscape types in the region.
High-frequency monitoring across changing Arctic catchments Long-term time series of Arctic hydrochemistry are rare, and are usually based on relatively low-frequency sampling at a limited number of field sites (Laudon et al. 2017). These limitations are often a result of logistical and financial constraints. Additionally, in nearly all Arctic riverine solute flux studies, individual storm or flow events are rarely incorporated or represented within biogeochemical flux estimates. The long-term data that does exist demonstrates that hydrochemical fluxes of carbon and nutrients are largely increasing from catchments across the Arctic, though this effect depends on the solute in question (e.g., McClelland et al. 2007McClelland et al. , 2016Toohey et al. 2016), in addition to landscape features that control solute export (Tank et al. 2016). The generality of increasing material flux through rivers implies that lateral nutrient losses from terrestrial environments are increasing via multiple drivers, including shifting hydrologic seasonality, increasing frequency of precipitation events, permafrost degradation, and the restructuring of both landscape topography and terrestrial vegetation (see reviews Bring et al. 2016, Wrona et al. 2016. Accurately quantifying riverine loads across variable flow conditions and landscape types now and into the future will be crucial for informing global change research and earth systems models (Kicklighter et al. 2013;Fan et al. 2019).
First, our results underscore that incorporating carbon and nutrient fluxes during high-flow events will help constrain uncertainty in Arctic biogeochemical budgets (Kicklighter et al. 2013). In our study catchments specifically, storm exports represented between 35% and 75% of total annual DOC exports, though these estimates depend primarily on capturing sporadic precipitation events across the thaw season (Shogren et al. 2020). Indeed, many study regions outside of the highlatitudes have benefitted from the implementation of highfrequency water chemistry monitoring and the analysis of C-Q responses to understand catchment hydrochemical behaviors and responses to perturbations (Burns et al. 2019;Godsey et al. 2019). Despite rapid changes to high-latitude landscapes, efforts in these regions have lagged behind those in temperate catchments (Laudon et al. 2017). In the very near future, the mobilization, processing, and export of carbon and nutrients are expected to change across high-latitude regions (Vonk et al. 2019;Tank et al. 2020). As a direct result of climate change, larger flow pulses, enhanced subsurface flow, and an extended season are likely to lead to greater lateral and longitudinal solute flux (Treat et al. 2016;Wrona et al. 2016;Vonk et al. 2019), but there is little infrastructure or means in place to assess these expectations. High-frequency frameworks offer a solution to constraining biogeochemical budgets in a changing Arctic, such as enhanced capture of biogeochemical flux variability across more of the Arctic thaw season (Shogren et al. 2020).
Second, the high-frequency observations from this study demonstrate how characteristics of the catchment influence the hydrodynamic responses of biogeochemical constituents transported within Arctic river networks (Harms et al. 2016;Connolly et al. 2018;Tank et al. 2020). Therefore, catchment and seasonal context are necessary to predict changes in export of DOC and NO 3 − at larger spatial scales, especially predictions of export under different climatic and ecosystem conditions predicted for the near future (Abbott et al. 2016;McGuire et al. 2009). The C-Q differences occurring across the landscapes of this study also imply that attempts to model or scale up Arctic riverine flux estimates must incorporate both landscape heterogeneity and seasonality into model testing and validation, rather than assuming a single C-Q relationship for DOC or NO 3 − through time and across high-latitude sites.
While additional work is needed to parse the sourcing of materials across distinct Arctic landscapes within and beyond our study catchments (Metcalfe et al. 2018), using a highfrequency C-Q framework presents an opportunity to constrain and identify the processes controlling lateral carbon and nutrient losses in diverse Arctic river networks. It will also provide new data to assess predictions of Arctic landscape change and inform future regional and earth system models. Therefore, we encourage investment in high-frequency C-Q monitoring and analysis efforts to expand understanding of biogeochemical conditions in a changing Arctic.

Data availability statement
The hydrochemistry data used for this study are freely available in the Arctic LTER Environmental Data Initiative (EDI) portal: Kuparuk (https://doi.org/10.6073/pasta/