Multi-year, spatially extensive, watershed scale synoptic stream chemistry and water quality conditions for six permafrost-underlain Arctic watersheds

. Repeated sampling of spatially distributed river chemistry can be used to assess the location, scale, and persistence of carbon and nutrient contributions to watershed exports. Here, we provide a comprehensive set of water chemistry measurements and ecohydrological metrics describing the 20 biogeochemical conditions of permafrost-affected Arctic watersheds. These data were collected in watershed-wide synoptic campaigns in six stream networks across northern Alaska. Three watersheds are associated with the Arctic Long-Term Ecological Research site at Toolik Field Station (TFS), which were sampled seasonally each June and August from 2016 to 2018. Three watersheds were associated with the National We

may not reconcile (Kareiva and Andersen, 1988) because connectivity among patches can create 60 emergent patterns and processes (Sivapalan, 2003;McDonnell et al., 2007;Covino, 2017). To understand and predict ecosystem behavior in the Anthropocene, we need to observe how biogeochemical patterns are produced and propagate across scales. Here, we describe a medium scale watershed chemistry dataset that includes spatially distributed hydrological, ecological, and geochemical properties. Using a synoptic experimental design, we measured these parameters across 65 medium-scale watersheds (<1 to >1000 km 2 ) multiple times over several years. We hope this dataset will help bridge the gap between plot-level and regional investigation of ecosystem change in the permafrost zone.
Most water chemistry and flow assessments conducted in the Arctic and elsewhere are based on observations at river outlets (McClelland et al., 2006(McClelland et al., , 2007Tank et al., 2016;Toohey et al., 2016;70 Shogren et al., 2020;Zarnetske et al., 2018). The flow of water integrates biogeochemical signals, such that river chemistry at the watershed outlet contains information about both terrestrial and aquatic biogeochemical processes that occurred upstream in the network (Temnerud et al., 2010;Vonk et al., 2019;Tank et al., 2020). Indeed, using sampling and monitoring approaches that capture the watershed outlet response over time has logistic and safety advantages for site access. Further, the recent 75 application of novel sensor technology has enabled high-frequency watershed-scale studies (Shogren et al., 2021;Ruhala and Zarnetske, 2017;Khamis et al., 2021). For example, the paired high-frequency flow and a limited set of chemical properties for the watersheds in this data paper are available at the Arctic Data Center (Zarnetske et al., 2020b, c, a). While these watershed outlet measurements can provide insight into possible upstream and upslope processes (Laudon et al., 2017;Shogren et al., 2021;80 Moatar et al., 2017), they often do not diagnose primary drivers of lateral transport of materials (Burns et al., 2019;Appling et al., 2018;Temnerud et al., 2010;Hoffman et al., 2013;Collier et al., 2018).
These large-scale measurements are the result of variable inputs which are buffered and blurred as multiple spatiotemporal signals are mixed and propagated through the surface and subsurface network (Creed et al., 2015;Abbott et al., 2018;Kolbe et al., 2019). To identify the processes behind those 85 signals, we need to venture into the headwaters, extending our observations into smaller subcatchments that match the spatial scale of mechanisms controlling carbon and nutrient uptake and release (Shogren et al., 2019;Gu et al., 2021;Abbott et al., 2017).
Spatially extensive or "synoptic" sampling frameworks, such as contained in this data paper, provide multiscale information about the source of signals across the entire watershed network, creating 90 a direct complement watershed outlet monitoring. With a synoptic sampling design, researchers can capture the spatial extent of nested subcatchments and therefore assess terrestrial-aquatic transfer of material and stream network processing Shogren et al., 2019;Gu et al., 2021).
Because permafrost degradation is triggering both large-scale deepening of the active layer and discrete permafrost collapse (thermokarst) features (Gao et al., 2021;Turetsky et al., 2020;Farquharson et al., 105 2019), synoptic snapshots could be invaluable in detecting the degree, location, and type of climate response. Therefore, measuring the spatial distribution of water chemistry in high latitude river networks could advance understanding of permafrost ecosystems and improve estimates of ecosystem feedbacks to climate change (Bring et al., 2016;Wrona et al., 2016;Schuur et al., 2015;Mu et al., 2020). 110 The datasets presented here were derived from repeated synoptic samplings in six Arctic watersheds in northern Alaska occurring on three distinct high latitude ecosystem types: Arctic tundra, Boreal forest, and Alpine tundra ( Figure 1). In this paper, we illustrate the utility of such data via a set of initial watershed chemistry analyses for ecologically significant reactive solutes including dissolved organic carbon (DOC), nitrogen (e.g., nitrate, N-NO3 -; ammonium, N-NH4 + ; dissolved organic nitrogen, 115 DON; total dissolved nitrogen, TDN), phosphorous (soluble reactive phosphorus, SRP; total dissolved phosphorus, TDP), as well as a suite of geochemically significant anions and cations (e.g., calcium, Ca 2+ ; total iron, Fe; dissolved silica, DSi; see Table 1 for full list of analytes). In addition, we use these datasets to introduce simple metrics for biogeochemical solutes: variance collapse, subcatchment leverage, and spatial persistence Shogren et al., 2019;Gu et al., 2021;Dupas et al., 120 2019;Frei et al., 2020). These new metrics seek to extract information more fully from rich spatiotemporal water chemistry datasets. Specifically, these metrics characterize what spatial scale is the most relevant in explaining terrestrial-aquatic material flux, how much influence or leverage each sampling site has on the watershed budget, and whether individual samplings are adequate to capture temporal variation. In this light, synoptic sampling frameworks provide robust information about how to 125 scale plot-and reach-level observations while also providing a multi-scale target for remotely sensed data and numerical models. Ultimately, the information gleaned from these metrics is desired by a range of disciplines from ecologists to natural resource managers.
First, we use subcatchment leverage to identify nested areas within the network that exert a disproportionate influence on flux at the watershed outflow Shogren et al., 2019). 130 Subcatchment leverage can be interpreted as the contribution of the subcatchment to watershed mass flux where the value can be negative (indicating the subcatchment has lower areal flux than the outlet, decreasing watershed flux), positive (indicating the subcatchment has higher areal flux than the outlet, a net increase in flux), or zero (no influence because it is the same as the outlet). Estimating leverage allows identification of specific subcatchments with disproportionate influence on material export, 135 defined here as high leverage. Subcatchments with high leverage behave as a strong source or sink within the watershed network, strongly influencing the resulting concentrations at the outflow, and can be selected as sites for further mechanistic study or monitoring. Likewise, the direction and magnitude of leverage averaged across the entire watershed contains information about net solute removal and production in the stream network (Shogren et al., 2019). For example, if the mean leverage for the 140 watershed is above zero, this indicates there are more solute sources than can be accounted for at the watershed outlet, implying there has been solute removal during transport through the network. Second, we examine what spatial extent or patch size controls solute production and removal by identifying thresholds of concentration variance collapse . We generally expect the amplitude of solute variability to decrease moving downstream from headwaters to larger systems (Creed et al., 145 2015) with greater variability among headwaters, whereas downstream reaches are less likely to have extremely high or low concentrations because they integrate multiple upstream source or sink processes (Wolock et al., 1997;Temnerud and Bishop, 2005;Burt and Pinay, 2005;Abbott et al., 2017).
Therefore, the size of nutrient sources and sinks in the landscape can be assessed by the spatial scale of the variance collapse of concentration among watershed reaches Shogren et al., 150 2019). The threshold of variance collapse is similar to the elementary representative area concept (Zimmer et al., 2013, p.20), where the threshold represents the spatial scale at which landscape "patches" or processes throughout the watershed network that produce and remove solutes are effectively integrated. Lastly, the spatial persistence metric can be used to assess whether a given site is representative (i.e., the same pattern continues through time), or if patches restructure in space between 155 sampling campaigns (i.e., reorganization of patches requires greater frequency in sampling) Dupas et al., 2019;Gu et al., 2021). Spatial persistence effectively quantifies the temporal representativeness of an instantaneous measurement at a given site, potentially indicating the type of process creating the patterns and informing future watershed study design and data analysis of extant data (Kling et al., 2000;Shogren et al., 2019). Creek. The three study watersheds were chosen because they spanned dominant circumarctic vegetation types, permafrost characteristics, and hydrologic conditions (Table 1). Further, the climate, morphology, and ecology of the sites and region have been previously described (Hobbie and Kling, 2014).
• The Kuparuk River (68.64816, -149.41152, Figure 2A) is a meandering stream flowing 170 through primarily tundra vegetation, located about 10 km northeast of TFS. The Kuparuk River includes a long-term monitoring site for the Arctic LTER, used as a site for ecological study and monitoring since 1979. From 1983-2016, the 4 th order reach of the Kuparuk River was used for a whole-stream fertilization study (Peterson et al., 1993;Slavik et al., 2004;Iannucci et al., 2021), where phosphorous (H3PO4) was continuously added to assess 175 response to nutrient fertilization. As the Kuparuk River continues north, it meets a large aufeis (ice) field (Yoshikawa et al., 2007;Terry et al., 2020).
• Trevor Creek (68.28482, -149.350063, Figure 2C) is a mountainous alpine stream, draining into the Atigun River watershed, located 30-km south of TFS. Trevor Creek drains primarily steep, rocky slopes with limited heath and willow vegetation. The majority of stream runoff is generated by precipitation and snowmelt. 185 As a result of long-term study and a sustained commitment to data stewardship, the Arctic LTER and TFS hosts an extensive catalogue of terrestrial, aquatic, and atmospheric data that are complementary to the data presented in this publication. For more information, please see the LTER data catalogue (https://arc-lter.ecosystems.mbl.edu/data-catalog), in addition to the abiotic and biotic monitoring data from the TFS Spatial and Environmental Data Center (https://toolik.alaska.edu/edc/index.php). 190

National Parks Service and U.S. Geological Survey Sites
We also sampled three watersheds associated with the National Park Service (NPS) Arctic Inventory and Monitoring Network and a project funded by the U.S. Geological Survey's (USGS) Changing Arctic Ecosystem program. The Agashashok and Cutler River watersheds are within Noatak National Preserve and the Akillik River watershed is within Kobuk Valley National Park. All three watersheds 195 are situated near the northern extent of Alaska's boreal forest, where tree line is expanding (Suarez et al., 1999), and subcatchments vary in areal extent of forested versus tundra land cover. The study sites vary with respect to permafrost characteristics, including soil texture, ground ice content, and subsurface hydrology (O'Donnell et al., 2016). Evidence suggests stream chemistry varies across these watersheds, including the form, amount, and age of dissolved carbon (O'Donnell et al., 2020). 200 • The Cutler River (67.845, -158.316, Figure 3A) flows north out of the Baird Mountains through gently rolling tundra into the upper Noatak River. The watershed is underlain by icerich glaciolacustrine deposits (O'Donnell et al., 2016), and soils tend to be organic-rich and poorly drained. Vegetation is dominated by moist acidic tundra and wet sedge meadows.
• The Akillik River (67.201, -158.572, Figure 3B) flows south out of the Baird Mountains and 205 into the Kobuk River downstream of the village of Ambler, Alaska. The river passes through alpine terrain in the headwaters before draining terrain comprised of ice-rich loess in the lower reaches. Vegetation is a mixture of boreal spruce forests and tundra.
• The Agashashok River (67.268, -162.636, Figure 3C) is a braided, clearwater river that flows from the northeast to southwest into the lower Noatak River north of Kotzebue, Alaska. The 210 headwaters drain rocky, alpine tundra terrain of the western Brooks Range. Downstream, the river drains broader valleys with a mixture of boreal spruce forest and tundra vegetation.

Arctic LTER Sites
Our sampling of the TFS watershed networks was designed to capture 30-50 "nested" subcatchments within the Kuparuk River, and Oksrukuyik and Trevor Creeks. Site selection was based primarily on (1) presence of flowing surface waters, (2) representation across varying subcatchment drainage areas, and (3) site accessibility. Often, we a priori chose sites located at subcatchment confluences, sampling both 220 upstream locations and then downstream of river mixing. In each of the TFS watersheds, we performed 5 repeated synoptic campaigns, sampling each stream network in August 2016, June 2017, August 2017, June 2018, and August 2018 (exact dates in Table 2). We accessed sampling sites either on-foot or by helicopter within a 6-hour period.

NPS/USGS Sites 225
Sampling of the NPS/USGS watershed networks was designed to capture ~5-10 subcatchments within the Agashashok, Cutler, and Akillik Rivers. Sites were selected to span a gradient of size (subcatchment area, stream order), vegetation (forest vs. tundra), and permafrost characteristics (parent material, ground ice content). Due to variation in watershed aspect, streams also spanned a spatial gradient in permafrost ground temperatures, areal extent, and active layer thickness Sjöberg et 230 al., 2021). In addition to stream chemistry parameters, stream discharge was measured, and samples were collected to characterize stream biota (benthic biofilm, macroinvertebrates, and resident juvenile fish).
In each of the NPS/USGS watersheds, we performed 4-10 repeated synoptic campaigns,  Table 2). We accessed sampling sites by helicopter within a 24-to 96-hour period.

Synoptic Site Characterization 3.1.1 Subcatchment Delineation for Drainage Area 240
The location of each stream sampling site was recorded in a spreadsheet and imported into GIS software (ESRI ArcGIS v. 10.4). These sites served as starting points ('pour points') from which watersheds and subcatchments were delineated following the general procedure described here: (https://support.esri.com/en/technical-article/000012346). The following two digital elevation models (DEMs) were needed to cover the spatial distribution of the stream sampling sites and were used to 245 create the necessary flow direction and flow accumulation layers: ArcticDEM from the Polar Geospatial Center (Porter et al., 2018) and ASTER GDEM v.2 (NASA/METI/AIST/Japan Spacesystems and US/Japan ASTER Science Team, 2009). A Python script was written to iterate over the list of sample sites and execute the watershed delineation procedure.

Estimation of terrestrial catchment characteristics for TFS sites 250
We characterized the terrestrial environment of the TFS sites using remotely sensed data pertaining to the vegetation and topography of each subcatchment. For each subcatchment polygon, we extracted the mean, standard deviation, and range of the elevation, slope, and topographic position index (i.e., the elevation of a given pixel relative to surrounding pixels, sometimes known as slope position). These metrics were calculated from 25-meter-resolution elevation data retrieved from the USGS National Map 255 website (https://viewer.nationalmap.gov/basic/). The normalized difference vegetation index (NDVI), which indicates the presence of green vegetation, was derived from imagery acquired in summer 2012 by the ETM+ sensor on Landsat 7 (courtesy of the USGS). We also extracted percent cover of vegetation classes in each subcatchment from the 30-meter-resolution Jorgenson northern Alaska ecosystems map (Muller et al., 2018). All data extraction was performed using zonal statistics via 260 ArcPy (ESRI, 2016) in Python.

Field sample collection & preparation 3.2.1.1 Arctic LTER
During each synoptic campaign, at each site we measured in-situ physiochemical variables (this section) 265 and sampled stream surface water for chemical analysis (section 3.2.2). All physical water samples were "grab" sampled directly from the stream thalweg, or as close to mid-channel as could be safely accessed. We collected samples in acid-washed and triple-rinsed 1-L amber PCTE bottles. We used handheld YSI ProPlus multiparameter probes (YSI Instruments Part No: 626281) and YSI ProODO Dissolved Oxygen Meter (YSI Instruments Part No: 6050020) to measure specific conductance 270 (µS/cm), pH, temperature (ºC), and dissolved oxygen (DO, in % saturation and mg O2/L) at each sampling site. We placed the probe into the water column where the water sample was taken and waited for the temperature and DO readings to stabilize before recording the final value.
Upon returning to the lab at TFS, we processed each water sample into aliquots for specific analytes within 8 hours of collection. We lab-filtered samples for dissolved water chemistry and 275 nutrients using handheld 60 mL syringes. We triple-rinsed syringes with unfiltered sample water. Then, we sparged each filter cartridge with ~10 mL of sample water prior to sample filtration; we used the sparge volume as the initial bottle rinse. We filtered samples for DOC/TDN into triple-rinsed amber 60-mL HDPE bottles using a 25 mm 0.2 µm cellulose acetate filter (Sartorius CA membrane, 11107-25-N).
We filtered samples for dissolved nutrients, anions, and cations into triple-rinsed clear HDPE 60-mL 280 bottles using a 47 mm 0.7 µm glass fiber filter (Whatman GF/F, 1825-047). Additionally, we placed ~60-mL of unfiltered sample water into a clear HDPE bottle for analysis of turbidity (NTU) and alkalinity (mg CaCO3/L). After processing, we froze samples at -4 ºC until analysis, except for aliquots for DOC and total dissolved nitrogen (TDN). We stored DOC/TDN samples at 2 ºC until analysis.
Samples were shipped express to the University of Vermont (UVM) and Brigham Young University 285 (BYU) for further analysis.

NPS/USGS
While sample collection and processing were similar between the TFS and NPS/USGS field sites, the filtration step varied slightly. For NPS/USGS samples, we followed standard USGS protocols. We filtered all samples for nutrient, anion, and cation analysis using 0.45-µm capsule filters (Geotech 290 Veraspor dispos-a-filter) into 250-or 500-mL HDPE bottles. We filtered samples for DOC and TDN into 125-mL amber glass bottles. Samples for alkalinity and total Fe were left unfiltered. DIC samples were collected without filtering or any headspace in 60-cc luer-lock syringes fit with two-way stopcocks. After processing, we froze samples at -4 ºC until analysis, with the exception of aliquots for DOC, TDN, and DIC that were stored at 2 ºC until analysis. Samples were shipped express to Oregon 295 State University's Cooperative Chemical Analytical Laboratory (CCAL; http://ccal.oregonstate.edu/) or the USGS in Boulder, Colorado, for further analysis.

Estimation of ecohydrological metrics
In addition to reporting solute concentrations for each synoptic campaign (e.g., Figures 4-5), we estimated ecohydrological metrics for each nested site and watershed. Across these analyses, we assigned any value below detection as the values of half the limit of quantification and kept these data points in the analysis. When the sample was not run for a specific solute, the cell was left blank. 335

Subcatchment Leverage
First, we estimated subcatchment leverage from each of the synoptic sampling events for each solute.
Subcatchment leverage is calculated as the difference in terms of concentration at each site (Cs) from the concentration at the watershed outlet (Co), subcatchment area (As) relative to the entire watershed area (Ao), and specific discharge at the sampling location (q = Qs/As, where Qs and As are the discharge 340 and subcatchment area at the sampling point): In the case of Eqn. 1, leverage is expressed in units of flux (mass/volume/time). However, if specific discharge is unavailable for each sampling location leverage can be estimated using only variability in concentration and subcatchment area, so long as specific discharge (q) is similar between 345 subcatchments (Asano et al., 2009;Karlsen et al., 2016). With the exception of the Agashashok River, which has flow generated from deeper flowpaths, our study watersheds have very little regional groundwater influence (Lecher, 2017), and the synoptic campaigns were performed near base-flow conditions. Therefore, for the purposes of this study, we assumed that q was similar for subcatchments within a study watershed, but not necessarily across the six study watersheds. This assumption was 350 tested at all Arctic LTER sites using dilution gauging at a subset of sites in summer 2018 and 2019, where we found that values of specific discharge were similar across subcatchment sizes (Shogren, unpublished data). We used Eqn. 2 to estimate subcatchment leverage for all sampling locations across sampling events: Here, subcatchment leverage has units of concentration, or percentage when normalized to outlet concentration. In other words, the distributed mass balance is relative to a pre-determined outflow point; with the presented analysis, we used the furthest downstream point with the largest drainage area as our "outlet" location. The interpretation of leverage values is opposite at the site and watershed scales Shogren et al., 2019). For example, a site with a positive value for subcatchment 360 leverage is contributing more than the typical subcatchment in the watershed. Conversely, a watershed with a mean leverage value that is positive is indicative of a net removal in the stream network because there is more solute in the tributaries than can be accounted for at the watershed outlet, while a negative value suggests solute production in the network Shogren et al., 2019). We report both mean leverages for each catchment (presented in Figures 6 and 7) and site-specific subcatchment 365 leverages for each solute (Figure 10 for DOC and NO3 -, but all other solutes can be found within the ecohydrological metrics datasets).

Concentration Variance Collapse
Next, to assess the representative "patch" size where concentration variance is reduced, we determined the threshold of concentration variance collapse for each solute from each synoptic sampling event 370 (shown in Figure 8). Using concentrations plotted over watershed area, we used the 'changepoint' package in R (Killick and Eckley, 2014) to determine the collapse in variance of concentration across the whole watershed area. To determine the reduction in variance statistically, we used the pruned exact linear time (PELT) method, which compares differences in data points to determine statistical breakpoints Shogren et al., 2019). We performed this analysis using scaled 375 concentrations, which were scaled by subtracting the whole watershed mean and dividing by the standard deviation to facilitate comparison of changes in variance and evaluate convergence towards the watershed mean. The variance collapse threshold is therefore expressed in units of area (here as km 2 ). A non-significant variance collapse threshold can be interpreted to mean either the processes controlling lateral fluxes are operating at too small or too large of a scale to be captured using a subcatchment 380 sampling approach.

Spatial Persistence
Lastly, we analyzed this spatially rich synoptic data to quantify the spatial persistence of stream nutrient concentrations and to determine the level of sub-grid resolution necessary to represent controls on lateral nutrient loss. The spatial persistence metric indicates whether spatial sampling is 385 representative or whether spatial patterns "reshuffle" over time. Spatial persistence (rs) is calculated as: Where rgx is the rank of subcatchments at the time of synoptic sampling, rgy is the rank of the long-term flow weighted concentrations, while srgx and srgy are the standard deviation of the rank variables. We calculated spatial persistence using the correlation function in R (Version 3.3.0), using the Spearman 390 method Shogren et al., 2019). Significance was tested using a student's t distribution test. Additional methods for calculating spatial persistence have now been proposed that do not require discharge data for the flow weighting (Gu et al., 2021). For the purposes of the Arctic LTER analysis, we estimated spatial persistence as the Spearman's correlation between Early (June) and Late (August) site concentrations, resulting in a single spatial persistence metric (rs) for 2017 and 2018. For 395 the NPS/USGS sites, spatial persistence was calculated as the correlation between site locations sampled in the Early (June) and Mid (July) and the Mid to Late (August or September) seasons.

Use and interpretation of ecohydrological ecosystem metrics
The original intent of this manuscript was to present our unique Arctic datasets and showcase the utility of a synoptic framework in combination with metrics that describe the spatial distribution of river 400 chemistry. To further highlight how these metrics can inform future sampling design and address fundamental ecological questions, below we describe patterns for DOC and NO3in the TFS watersheds.
For solutes, the spatial variability in concentration depends on the strength and connectivity of both source and sink patches superimposed on the structure of the stream network .
When we plot solute concentration against subcatchment area, we find more variability water chemistry 405 in smaller subcatchments (<30 km 2 ). This can be interpreted as a spatial "fingerprint" and is shown most clearly in Figure 10, which displays the spatial distribution of DOC and NO3concentrations across watersheds and sampling campaigns. Generally high concentration variability in smaller headwaters, which converges to mean watershed behaviour towards the catchment outlet holds with the conceptualizations of large rivers as "chemostats" (Creed et al., 2015). In the context of Arctic 410 watersheds, these concentration/area relationships reveal consistently high DOC and low NO3concentrations in the low-gradient tundra watersheds (Kuparuk River and Oksrukuyik Creek), despite high variability in smaller contributing subcatchments. In contrast, the alpine watershed Trevor Creek, has relatively low DOC and high NO3concentrations, likely due to shorter and faster hydrologic flowpaths and lower terrestrial biomass (Shogren et al., 2019). Overall, these findings are consistent 415 with studies that indicate that slower, longer flowpaths and productive terrestrial vegetation control carbon and nutrient transfer and mobilization in lower-gradient tundra watersheds (Shogren et al., 2019(Shogren et al., , 2021). If we assume that spatial variability in stream network water chemistry depends primarily on the extent and connectivity of upstream sources/sinks, then the patches sizes that control solute fluxes can be assessed by the spatial scale of the variance collapse Shogren et al., 2019). 420 Across all three TFS watersheds, the generality of variance collapse at intermediate scales is indicative that subcatchment scale "patches" (~10-50 km 2 ) control whether carbon and inorganic nitrogen is produced or removed at the watershed scale ( Figure 10). In addition, the consistency of the thresholds across sampling campaigns (Figure 8 and 10) highlights the importance of capturing intermediate scale biogeochemistry to bridge understandings from plot-level experimentation to larger more regional-scale 425 observations (Shogren et al. 2019).
When we convert concentrations into estimates of subcatchment leverage (e.g., Figure 6, 7, 10 and Figure 11), patterns emerge that further contextualize the spatial distribution of DOC and NO3concentrations. First, we can investigate whole watershed ("net") behaviour by calculating the mean leverage and examining the distribution of values with boxplots (as in Figures 6 and 7). As a more 430 specific example, mean NO3leverage within the Kuparuk watershed ( Figure 6D, second row) were consistently above zero (note the reversed axis), revealing strong removal or retention before it reached the watershed outlet, which is consistent with high biotic N demand. Within this same watershed, DOC leverage values were often at or just above the zero line ( Figure 6D, first row), representing primarily conservative transport of DOC (i.e., no net production or uptake). Within the lake-influenced 435 Oksrukuyik watershed, NO3leverage values were more variable (i.e., leverage above/below zero-line; Figure 6E, second row), implying a combination of removal and production mechanisms acting across the watershed network. When visualized as "net" behaviour, the watershed and season-dependent directionality of net leverage patterns are congruent with emerging evidence that landscape template exerts strong control on biogeochemical signals in Arctic rivers (Vonk et al., 2019;Tank et al., 2020;440 Shogren et al., 2021). As a compliment to the first approach, we can additionally examine individual subcatchment leverage values to reveal the effect of each contribution on what we observe at the watershed outlet. This can be interpreted similarly to statistical leverage, where one or more points may exert high influence on a linear regression.
Across all TFS watersheds, there are a few select subcatchments that contribute 445 disproportionately to DOC fluxes, while the more variable patterns for NO3suggest additional spatial and seasonal controls (Figure 11). For example, patterns in the Kuparuk River and Oksrukuyik Creek (Figure 11a-b) could be interpreted to mean that DOC is transported conservatively in lower gradient landscapes, while lateral fluxes of NO3are more tightly controlled by biotic demand (Harms et al., 2016;Khosh et al., 2017;Connolly et al., 2018;Kendrick et al., 2018;Iannucci et al., 2021). Across 450 solutes and watersheds, the information gleaned from the leverage metric is useful in several ways.
First, subcatchment leverages allow for the direct identification of watershed areas that are disproportionately driving carbon and nutrient exports. For any chosen solute or suite of materials, sites identified as "high leverage" indicate strong source/sink behaviour, which could be (1) validated with regular field observations that relate riparian or terrestrial conditions with empirical measurements of 455 water chemistry, (2) selected for further study designed to identify the abiotic and biotic mechanisms that drive patterns of riverine chemistry, and/or (3) identified as non-representative sites relative to proximal subcatchments of similar size and terrestrial characteristics. Relatedly, estimating subcatchment leverage enables researchers to identify sites that are representative of watershed-scale behaviour, which could be used to more effectively scale biogeochemical dynamics in Arctic rivers 460 relative to outlying subcatchments (Kicklighter et al., 2013;Pinay et al., 2015;Aguilera et al., 2013).
Finally, the application of the simple spatial persistence metric can help researchers determine whether a sampling location is behaving consistently, or if solute contributions are moving in space across sampling events Dupas et al., 2019). In the context of work in remote watersheds, the ability for researchers to identify both stable and unstable processes presents an exciting 465 opportunity to ask questions about the consistency of subcatchment contributions and optimize sampling or experimental design. For example, DOC concentrations are generally spatially stable between early and late sampling events (rs > 0.50), particularly in the Kuparuk River and Trevor Creek watersheds (Figure 9). In these landscapes, a high rank correlation indicates that repeated sampling of the same location will result in a similar spatial distribution of concentrations. While sampling 470 repeatedly in the early and late seasons may reveal increases or decreases in solute concentrations (Shogren et al. 2019), the high degree of relatedness indicates that these patterns will be maintained across the watershed network. However, the low persistence (rs < 0.50) for DOC in the Oksrukuyik Creek watersheds signifies substantial spatial shifts across the early and late thaw season (Shogren et al. 2019). While there was variability in the persistence across watersheds and solutes, the stability metric 475 can be used by future researchers to identify whether sampling the same location repeatedly does or does not represent the spatial dynamics across sampling events.

Conclusions
With this work, we provide a detailed characterization of physical, chemical, and biological parameters that are essential to using river network chemistry to infer ecosystem-level carbon and nutrient balance.
We apply novel metrics to these data that describe the spatiotemporal patterns of watershed 485 biogeochemistry in six permafrost-underlain Arctic watersheds. These data represent a high-resolution and temporally replicated river chemistry dataset from understudied permafrost-dominated regions.
Combining these measures with remotely sensed data, plot-level experiments, and numerical models could advance our understanding of permafrost ecosystems in the face of climate change and other disturbance. 490 Author contributions: All co-authors participated in the field collection, laboratory analysis, and/or curation of the dataset, and contributed to writing and/or editing this paper. AJS was primarily responsible for writing this paper and assembly of the archival database.

495
Competing Interests: The authors declare that they have no conflict of interest.
Financial Support: The authors would also like to acknowledge their funding support: AJS (NSF-DBI-1906381, NSF-OPP-1916567); JPZ and TW (NSF-EAR-1846855, NSF-OPP-1916567); BWA (NSF-505 OPP-1916565), WBB (NSF-OPP-1916576, NSF DEB-1637459);and FI and AM (NSF DEB-1637459    and Cutler Rivers, as these thresholds were often non-significant.    Tables  Table 1. Summary of site characteristics for the watersheds where synoptic samplings were conducted. The descriptions are considered representative of the major landform types within the 865 TFS and NPS/USGS watersheds.