Revealing biogeochemical signatures of Arctic landscapes with river chemistry

Riverine fluxes of carbon and inorganic nutrients are increasing in virtually all large permafrost-affected rivers, indicating major shifts in Arctic landscapes. However, it is currently difficult to identify what is causing these changes in nutrient processing and flux because most long-term records of Arctic river chemistry are from small, headwater catchments draining <200 km2 or from large rivers draining >100,000 km2. The interactions of nutrient sources and sinks across these scales are what ultimately control solute flux to the Arctic Ocean. In this context, we performed spatially-distributed sampling of 120 subcatchments nested within three Arctic watersheds spanning alpine, tundra, and glacial-lake landscapes in Alaska. We found that the dominant spatial scales controlling organic carbon and major nutrient concentrations was 3–30 km2, indicating a continuum of diffuse and discrete sourcing and processing dynamics. These patterns were consistent seasonally, suggesting that relatively fine-scale landscape patches drive solute generation in this region of the Arctic. These network-scale empirical frameworks could guide and benchmark future Earth system models seeking to represent lateral and longitudinal solute transport in rapidly changing Arctic landscapes.

We syringe-filtered water samples through 25 mm 0.7 µm GF/F filters (Whatman, 0.7 um nominal pore size) into pre-rinsed dark (DOC) and opaque (all other solutes) 60 mL Nalgene HDPE bottles. After filtration, samples were frozen (-4°C) until analysis, with the exception of aliquots for DOC (stored at 2°C until analysis). Prior to refrigeration, DOC samples were acidified with hydrochloric acid to a final normality of 0.01N and refrigerated. We measured NO2 -+NO3on a Lachat QuikChem FIA+ 8000 Series autoanalyzer (QuikChem Method 31-107-04-1-E; detection limit = 0.03 µM). We did not subtract NO2from these results, as concentrations are typically low to undetectable in these systems. SRP was measured on a Shimadzu UV-2600 spectrophotometer (method adapted from Parsons et al. 1984; detection limit = 0.05 µM). DOC was measured on a Shimadzu TOC-L CSH total organic carbon analyzer (combustion catalytic oxidation method; detection limit = 5.0 µM). When values were at or below the limit of detection (LOD), we treated values as a constant value of the LOD divided by 2 for analysis 2 .
During each synoptic campaign, we sampled stream water for chemical analysis and physiochemical variables. We collected grab samples of water from each field site in acidwashed 1L amber PCTE bottles. We then lab-filtered aliquots of the sample for specific analytes within 8 hours of collection following LTER protocols.

Calculating subcatchment metrics: variance collapse, leverage, and spatial stability
To assess the patch size of solute sources and sinks we determined the threshold of variance collapse for each solute from each synoptic sampling event. We ordered the subcatchments by watershed area and used the pruned exact linear time (PELT) method to identify significant increases or decreases in variance using the 'changepoint' package in R 3,4 .
Next, we estimated subcatchment leverage from the synoptic sampling events. or net production (negative mean leverage) at the watershed scale 4 . In addition to estimating mean subcatchment leverage, we mapped leverage estimated at each sampling point using 'ggmaps' in R 5 and spatial data available through the Toolik Field Station GIS and Remote Sensing center (https://toolik.alaska.edu/gis/data/). Lastly, we analyzed this data to quantify the spatial stability of stream nutrient concentrations and to determine the level of subgrid resolution necessary to represent controls on lateral nutrient loss. Spatial stability is calculated as the standard Spearman's rank correlation coefficient (rs) between instantaneous measurements and varies from -1 to +1. High values of spatial stability (value near 1) means that a single sampling campaign predicts the concentration taken at the next sampling campaign; conversely, a lower value indicates a mismatch between sampling dates and thus a spatial reorganization of stream hydrochemistry. We calculated spatial stability using the correlation function in R (Version 3.3.0) 4 .

Statistical analysis
We used a linear mixed effects model with site as a random effect to explore landscape predictors (e.g. slope, surface roughness) of surface water concentrations in each watershed using the package 'lme' in R 6 . Finally, we compared both mean solute concentrations and subcatchment leverage estimates between early and late season and among watersheds using 2way analysis of variance (ANOVA). When the ANOVA was significant, we performed a Tukey's Honesty Significant Differences (HSD) test to compare seasons among watersheds.  Table 2: Summary statistics for linear models for subcatchment concentrations among catchments, season, and topography (i.e., slope). Asterisks indicate parameter significance at p < 0.05 (*), p < 0.01 (**), and p < 0.001 (***), while (-) denotes no significance. Tukey's HSD for watershed differences are also reported HSD (p < 0.05).