Heterogeneous Patterns of Aged Organic Carbon Export Driven by Hydrologic Flow Paths, Soil Texture, Fire, and Thaw in Discontinuous Permafrost Headwaters

Climate change is thawing and potentially mobilizing vast quantities of organic carbon (OC) previously stored for millennia in permafrost soils of northern circumpolar landscapes. Climate‐driven increases in fire and thermokarst may play a key role in OC mobilization by thawing permafrost and promoting transport of OC. Yet, the extent of OC mobilization and mechanisms controlling terrestrial‐aquatic transfer are unclear. We demonstrate that hydrologic transport of soil dissolved OC (DOC) from the active layer and thawing permafrost to headwater streams is extremely heterogeneous and regulated by the interactions of soils, seasonal thaw, fire, and thermokarst. Repeated sampling of streams in eight headwater catchments of interior Alaska showed that the mean age of DOC for each stream ranges widely from modern to ∼2,000 years B.P. Together, an endmember mixing model and nonlinear, generalized additive models demonstrated that Δ14C‐DOC signature (and mean age) increased from spring to fall, and was proportional to hydrologic contributions from a solute‐rich water source, related to presumed deeper flow paths found predominantly in silty catchments. This relationship was correlated with and mediated by catchment properties. Mean DOC ages were older in catchments with >50% burned area, indicating that fire is also an important explanatory variable. These observations underscore the high heterogeneity in aged C export and difficulty of extrapolating estimates of permafrost‐derived DOC export from watersheds to larger scales. Our results provide the foundation for developing a conceptual model of permafrost DOC export necessary for advancing understanding and prediction of land‐water C exchange in changing boreal landscapes.

hydrologic regimes directly impact the structure and biogeochemical functioning of northern ecosystems (Rawlins et al., 2010;Wrona et al., 2016). Enhanced hydrological connection can promote greater transfer of C from terrestrial to aquatic ecosystems, leading to greater downstream organic matter decomposition that may increase emissions of carbon dioxide (CO 2 ) and methane from northern high-latitude aquatic ecosystems to the atmosphere (Schuur et al., 2015;Walter Anthony et al., 2016). Of particular concern, the thawing of ice-rich yedoma deposits in regions of Alaska, Siberia, and Canada can mobilize previously stored, ancient C into the active C cycle (Drake et al., 2015;Feng et al., 2013;Vonk et al., 2015).
Although permafrost thaw may deliver aged, dissolved organic carbon (DOC) to river networks (Feng et al., 2013;Neff et al., 2006;O'Donnell et al., 2020;Tank et al., 2016), detailed documentation of this process remains sparse for most of the arctic and boreal biomes. One understudied region of concern is the discontinuous permafrost zone, where ancient C is susceptible to release due to permafrost thaw and increasing wildfire . The Yukon River Basin (YRB) spans interior Alaska, U.S.A., and portions of the Yukon Territory and British Columbia, Canada (Figure 1), and is a representative basin within this region. In recent decades, broad scale changes in hydrologic flow paths, water source, and water residence times in the YRB have been inferred from changing river water and solute fluxes (Striegl et al., 2005;Toohey et al., 2016;Walvoord et al., 2012;Walvoord & Striegl, 2007). Hydrologic changes and permafrost thaw may impact DOC export and/or biological processing in multiple ways (Striegl et al., 2005;Tank et al., 2016;Zolkos et al., 2019); while there is currently no sign of ancient permafrost-derived DOC (hereafter "permafrost DOC") reaching larger rivers of interior Alaska (Aiken et al., 2014), elevated uranium isotope ratios, an indicator of permafrost thaw (Ewing et al., 2010(Ewing et al., , 2015, suggest that solutes released from localized permafrost thaw are entering tributaries in this region (Koch, Ewing, et al., 2013). Rapid mineralization of thawed permafrost DOC under experimental conditions and from isotopic studies (Drake et al., 2015;Ewing et al., 2015;Mann et al., 2015) indicates that ancient organic carbon (OC) cycling could be intense across headwater streams. On the other hand, large potential DOC yields (Wickland et al., 2018) and relatively low biodegradability of Holocene permafrost DOC in interior boreal Alaska (Textor et al., 2019;Wickland et al., 2018) point to a high potential for permafrost DOC transport along the aquatic continuum (Koch et al., 2021;. It is important to better quantify the patterns and drivers of ancient DOC mobilization within the headwaters of the YRB aquatic network in order to understand the fate of previously stored, aged terrestrial C, and to predict changes to aquatic ecosystems with continued thaw. Permafrost thaw and DOC mobilization in the YRB and other discontinuous permafrost regions in the boreal biome may be further accelerated by increasing fire frequency and intensity. Recent work in interior Alaska has shown that permafrost thaw is connected to recent fire occurrence using geophysical techniques (Minsley et al., 2016) and remote sensing . Linking the effects of thaw to catchment DOC export is complex, and highly dependent on soil type and fire history , and the interaction between infiltrating water and soils (Carey & Woo, 2001;Koch, Ewing, et al., 2013;Petrone et al., 2006;Prokushkin et al., 2007). Soil stratification often limits deep percolation, resulting in significant lateral flow along the organic-mineral soil boundary in the shallow subsurface Quinton & Marsh, 1999) that can quickly transport DOC to streams (Carey, 2003). The DOC exported to slower, deeper flow paths may be decomposed or sorbed to mineral particles before reaching the stream (Ågren et al., 2007;Koch, Runkel, et al., 2013;Prokushkin et al., 2007). In contrast, macropores and/or soil pipes formed as a result of thermokarst or fire may enhance rapid transport of water and DOC to streams (Carey & Woo, 2000Koch et al., 2014), and lead to thermal degradation of the frozen boundary and further enhancement in the transport of thawed materials to streams (Koch, Ewing, et al., 2013). Matrix flow along flow paths near the top of degrading permafrost may also serve as mobilization pathways of aged DOC . Because matrix flow is a slower pathway for solute transport than preferential flow, it can promote DOC sorption and microbial processing rather than export (Ågren et al., 2007;Striegl et al., 2005). Recent observational and modeling studies provide consistent support for the increasingly widespread development of shallow taliks (perennially unfrozen zones) as both depth to permafrost and the seasonal active layer depth increase (Connon et al., 2018;Parazoo et al., 2018;Streletskiy et al., 2015;Walvoord et al., 2019). Wildfire has been shown to enhance talik development in interior Alaska (Nossov et al., 2013;Rey et al., 2020). Yet, it is currently unclear how fire history interacts with other watershed features to influence aged DOC fluxes in headwaters of the YRB.
Transitional, discontinuous permafrost regions of the YRB are critical locations for permafrost DOC mobilization and processing (Aiken et al., 2014;Striegl et al., 2005). This region is highly susceptible to disturbance, given 10.1029/2021GB007242 3 of 16 climatic changes that are leading to increasing permafrost thaw, thermokarst, and talik formation Minsley et al., 2016;Nossov et al., 2013;Rey et al., 2020;Walvoord et al., 2019). Our goal is to understand how these disturbances impact the transport of DOC from terrestrial to aquatic systems. By capturing landscape heterogeneity in our sampling effort, we begin to address the potential for biases in the interpretation of aquatic DOC cycling patterns due to widespread underrepresentation of much of the circumpolar landscape in current sampling efforts (Bogard et al., 2019;Metcalfe et al., 2018). Here, we provide a broad evaluation of patterns and drivers of permafrost DOC input into aquatic networks within the transitional discontinuous permafrost zone of the YRB in interior Alaska (Figure 1). We explore how properties of individual catchments (varying combinations of permafrost extent, soil texture, and fire history) interact with seasonal soil thaw to shape the relative importance of aged DOC to the bulk surface water DOC pool. Few studies have directly linked patterns of headwater DOC radio-isotopic composition to hydrologic conditions (Campeau et al., 2019;Neff et al., 2006), especially across seasons for multiple catchments having varying soil properties, fire history, and permafrost extent. Therefore, this work provides a new level of detail to help explain environmental controls on potential global C-climate feedbacks (Schuur et al., 2015) that may facilitate greater C processing and emissions under changing climatic conditions in high-latitude, permafrost landscapes.

Description of Study Region
This study was conducted in the discontinuous permafrost zone of interior Alaska, within the boreal forest biome. For a detailed overview of the YRB, see Brabets et al. (2000). The region consists of rocky alpine zones in the Ray and White Mountains of the Yukon-Tanana Upland terrain, and silty ice-rich permafrost in valley bottoms that formed syngenetically during loess deposition in the Pleistocene (Péwé, 1975). Mean annual air temperature in the region ranges from <−5 to −2°C  and mean annual precipitation across the basin is 483 mm yr −1 (range: 254-3,302 mm yr −1 ) (Brabets et al., 2000). Soil landscapes are generally described as alpine rocky, upland rocky, upland silty, or lowland organic-rich terrain . Silty uplands and peaty lowlands tend to host open black spruce (Picea mariana) forests with groundcover dominated by moss (Sphagnum spp. and feathermosses) and lichens (Cladina spp.). Upland rocky sites are better drained and commonly host deciduous trees . Following low severity fire, spruce self-replacement is common, whereas deciduous species or mixed stand replacement tends to follow high severity fire (Johnstone et al., 2020). Thermokarst is present in the silty ice-rich landscapes in both upland and lowland settings . Differences in the hydrologic and thermal properties of rocky, silty, burned, and unburned soils were quantified by Ebel et al. (2019). The impact of thermokarst and fire on groundwater-surface water interactions and groundwater flows in the rapidly warming boreal Alaskan landscapes have been typically investigated based on the hydrology of large rivers and lowland lakes (Callegary et al., 2013). Contrary to findings for other major Arctic rivers, studies have not documented increases in total annual discharge of the Yukon River, although groundwater inflows are increasing (Toohey et al., 2016;Walvoord & Striegl, 2007). The hydrologic features of headwater streams are comparatively less studied.

Field Sampling Approach
Sites investigated in this study are headwater catchments located in interior Alaska along the Dalton and Steese Highways over a 210 km latitudinal gradient ( Figure 1) The eight catchments were chosen based on the presence of a headwater stream and selected to encompass a wide range of landscape types and subsequently varying soil properties, vegetation types, fire history, and estimated permafrost distribution (see Section 2.5 and Table 1). Many of the sites display complex hydrology, including the presence of ephemeral streams, seeps and gullies, macropores, and soil pipes that allow rapid drainage of slopes and transport of solutes and sediments. The probability distribution of permafrost to 1 m depth (Pastick et al., 2015) is shown.
These features, hereafter lumped as "tributaries", fed headwater streams (hereafter "streams"), defined by continuous summer flow, and better-defined channels and floodplains. Streams are useful integrators of catchment hydrologic processes due to their larger size (Fisher & Welter, 2005;Williamson et al., 2008), and may contain solutes and nutrients derived from runoff throughout the catchment. We sampled the main stream in six catchments to characterize surface water discharge. Tributaries to the catchment streams, due to their small size, are closely connected to the terrestrial landscape, and are thus more likely to contain solutes, nutrients, and carbon that have not been biogeochemically altered by in-stream processes (Aiken et al., 2014;Koch, Runkel, et al., 2013). Samples were collected from more than one tributary per catchment, when present, to capture variability in runoff chemistry. We sampled each stream (n = 6) and one or more flowing tributaries (n = 0-2) in each study catchment during each site visit, resulting in a total of 28 stream and 39 tributary samples. Sites were visited three times in 2016 (May, mid-July, and late August or early September) and twice in 2017 (May and September) to capture seasonal changes in stream chemistry associated with active layer thaw. Stage was measured continuously in the six streams using Rugged Troll 100 (In-Situ Inc.) pressure transducers and discharge was measured during each site visit using a wading rod and an acoustic Doppler velocimeter (Flowtracker 1, Sontek). The relationship between stage and discharge was used to determine continuous stream discharge, which is reported in Koch et al. (2020). Uncertainty in individual discharge measurements was typically around 8%, and continuous discharge was not well constrained due to the limited number of site visits and discharge measurements. Tributary discharge was difficult to quantify due to the small size of the features. Depending on conditions, we estimated tributary discharge via tracer dilution, wading measurements, or by filling a graduated container. Water quality parameters were measured in the field with a YSI EXO 2 sensor package and included water temperature, specific conductance (SpC), pH, and dissolved oxygen. Water samples were collected and filtered in the field using a peristaltic pump connected to a 0.45 μm Gelman Versapor capsule filter that was flushed with stream or tributary water prior to sample collection. Filtrate for isotopic analysis of DOC was collected in acid washed polycarbonate bottles and frozen for later analysis. DOC was collected in 40 mL amber glass vials that had been precombusted. Stable isotopes of water were collected in 12 mL borosilicate vials ensuring that there was no headspace in the vials. All samples were kept chilled until analysis.

Laboratory Analyses
DOC was measured by the platinum catalyzed wet oxidation method on an OI Analytical 700 TM Total OC Analyzer (Aiken et al., 1992) at the USGS laboratory in Boulder, CO. Frozen samples were thawed at room temperature and prepared for Δ 14 C-DOC measurement of DOC using standard methods fully detailed by Butman et al. (2012). Briefly, samples were acidified and sparged with pure N 2 gas to remove all inorganic C. Samples were then UV-oxidized to convert DOC to CO 2 , cryogenically purified, and CO 2 was trapped and sealed in 6 mm thick, combusted Pyrex glass tubes. Samples of CO 2 were processed at the National Ocean Sciences Accelerator Mass Spectrometry facility using standard methods. Values of Δ 14 C-DOC are reported in per mil notation and the mean age of the bulk DOC pool is determined using standard calculations following Stuiver and Polach (1977), and reported in units of years B.P. Stable isotopes of water collected in 2016 and May of 2017 were analyzed on a dual-inlet isotope-ratio mass spectrometer at the U.S. Geological Survey Reston Stable Isotope Laboratory. Stable isotopes of water collected in September of 2017 were analyzed with a Picarro L2130-i Isotope and Gas Concentration Analyzer  Table 1 Catchment Properties (Picarro) at the University of Washington. All water chemistry data from this study and additional details on analytical methods are publicly available in Foks et al. (2020).

Hydrologic Mixing Model
An endmember mixing model was used to determine the dominant sources of water in the streams. Previous work identified SpC and DOC as useful tracers of the source of boreal waters given that shallow soils are dominated by organic material (i.e., the dominant DOC source) and deeper soils are more mineral rich (the dominant SpC source) (Koch et al., 2014). Stable isotopes of water were also considered as a potential tracer, but ultimately rejected because of limited variability among observations and small differences in endmembers. Endmembers were chosen by plotting all of the stream chemistry in x-y space ( Figure S1 in Supporting Information S1) and comparing to endmembers from Koch et al. (2014) (precipitation, organic layer flow, and mineral layer flow) that were collected from one of the study catchments (West Twin Creek). We found that the previously determined precipitation endmember also effectively bounded the most dilute stream water measured in this data set, which is presumably derived from direct rainfall on the stream and rapid runoff through shallow soils and/or overland flow related to snowmelt and summer storms. The existing organic layer flow endmember for West Twin Creek (Koch et al., 2014) did not fully capture the wider range of variability observed in DOC concentrations and so needed to be adjusted upward. Finally, we determined that a previously unobserved solute-rich (i.e., high SpC and high DOC) endmember was needed to bound the stream data. Our highest SpC, DOC samples were measured in two seeps in the Richardson catchment, which burned in 2003 and is characterized by high silt content, has active thermokarst and shows geochemical evidence of permafrost thaw (Koch, Ewing, et al., 2013). We set the soluterich endmember equal to the mean plus one standard deviation from the Richardson seep samples, which bound all but one of the observed water chemistries in the mixing space ( Figure 2). One site, Globe 3, was excluded from the mixing analysis because the samples displayed unusually high SpC especially early in the summer, which we attributed to contamination from road salts. Upon defining the endmembers, we developed a three-component mixing model based on water and solute mass balances: where f i is the fractional contribution of each source to stream discharge and is the concentration of the tracer i (SpC or DOC) for endmember j (precipitation, organic, and solute-rich endmembers). The equations were solved in MATLAB (Mathworks) to determine the proportion of endmembers contributing to each stream sample.

Catchment Properties
Remote sensing and GIS analyses were conducted to quantify and summarize differences in catchment size, climate indices (i.e., mean annual air temperature and mean annual precipitation), soil landscapes, generalized geology, permafrost extent, topography (i.e., elevation, slope), burned area, and wetland and vegetative cover. Catchment boundaries were manually delineated using ridge crests evident on a USGS 60-m digital elevation model (DEM, National Elevation Data set (Gesch et al., 2002)). Climate indices were derived from downscaled, historical  climate data generated by the Scenarios Network for Alaska + Arctic Planning (https:// uaf-snap.org/). Mapping of soil landscapes and generalized geology was done through manual image interpretation following previously developed approaches (Jorgenson & Grunblatt, 2013;Jorgenson et al., 2009). Estimates of near-surface (within 1 m) permafrost occurrence (Pastick et al., 2015) served as a proxy for permafrost extent. Topography indices were developed using the DEM. Burned area was estimated using the fire-history-perimeter database  obtained from the Alaska Interagency Coordination Center (https://fire.ak.blm.gov/ predsvcs/maps.php). Estimates of vegetative and wetland cover were gathered from the 2016 National Land Cover Data set (NLCD [Dewitz, 2019]). Data for each property were extracted from within the catchment boundaries and summarized with ArcGIS (ESRI).

Statistical Analyses
The effects of hydrology, soil texture, and fire on values of Δ 14 C-DOC (and thus mean DOC age) were evaluated using Generalized Additive Models (GAMs). GAMs are tools to model nonlinear associations between predictor variables and responses, using the sum of unspecified smooth functions to estimate trends and have been widely used recently to model complex data sets (Webb, Hayes, et al., 2019;Wiik et al., 2018). They are a useful tool because they allow for flexibility between the response and predictor variables by not assuming a linear relationship. A Gaussian distribution was used for the response variable to account for the negative Δ 14 C-DOC values. Model fit was determined by examining the deviance and distribution of the residuals, comparing the residuals against the linear predictor, and the response versus fitted values. Basic statistical analyses (OLS regression) were conducted in R (version 4.0.3 [R Core Team, 2020]) using the lm function. GAMs were estimated using the mgcv package (version 1.8-33 [Wood, 2011]), and graphics were plotted with package ggplot2 (Wickham, 2009) for R.

Catchment and Stream Properties
The study watersheds varied across a broad range of landscapes, from alpine to lowland, rocky to silty, low to high fraction burned (Table 1). Catchment areas spanned from 48 to ∼7,000 ha, ranged over 350 m in elevation, and contained variable proportions of permafrost cover (30%-87%), burn area (0%-88%), and vegetative and wetland cover (Table 1). Mean instantaneous discharge in the six continuously monitored study streams during May to September of 2016 and 2017 ranged from 0.08 to 0.65 m 3 s −1 . Instantaneous discharge in eight tributaries ranged from 7.0 E−5 to 2.4 E−2 m 3 s −1 . Total summer discharge was 2.0 ± 0.8 times higher in 2016 relative to 2017 for the three streams with full records over both summers (from 20 May to 30 August). Based on summer precipitation measured at the Fairbanks, AK airport, 2016 was the fifth wettest summer during 1950-2021, with a total rainfall of 270 mm, which was 1.7 times higher than approximately average summer rainfall in 2017 (163 mm, relative to a 72-year mean of 152 mm).

High-Solute Loads in Endmember Mixing Model
Across sites and seasons, we observed a wide range in stream DOC, SpC, and δ 18 O-H 2 O. DOC concentrations ranged from ∼3.3 to 45 mg L −1 , SpC from ∼10 to 830 μS cm −1 , and δ 18 O-H 2 O from −23.6 to −17.9‰. We found that once the active-layer soils had thawed (i.e., excluding the May samples), streams generally fell on a mixing line between precipitation and the solute-rich endmember (except for the West Fork of Dall Creek; r 2 = 0.82, p < 5E−6; Figure S2 in Supporting Information S1), whereas tributaries fell on a mixing line between either precipitation and organic soil or organic soil and the solute-rich endmember. All chosen endmembers (Table 2) had similar δ 18 O-H 2 O concentrations ( Figure S2 in Supporting Information S1) indicating limited utility of δ 18 O-H 2 O as a tracer in this system during the summer.
The endmember mixing model indicated a broad distribution in the importance of the high solute endmember to the water chemistry of each stream (Figure 3). Eight streams or tributaries from three different catchments displayed contributions from the solute-rich endmember of over 20%. These three catchments, Erickson,   Koch et al., 2014) and Updated Values Figure 3. Endmember fractions for each stream and tributary, where streams are denoted by capital letters corresponding to watersheds in Figure 1 and Table 1, and "t" and "2" indicate a primary and secondary tributary, respectively. The red dashed line indicates the division between rocky sites (to the left), and silty sites (to the right), which tend to have higher solute-rich endmember contributions.
Richardson, and Isom, also displayed thermokarst features and channel incision. In two subcatchments (Erickson unburned tributary and Isom seeps), high solute endmember fractions existed despite no impacts of recent fire in the immediate catchment. In the Erickson watershed, we sampled two similar tributaries in burned and unburned regions, and the unburned tributary displayed a higher fraction of the high solute endmember than the burned tributary.

Linking Catchment Structure, Flow Paths, and Δ 14 C-DOC
Seasonal changes in the values of Δ 14 C-DOC and DOC mean age were not uniform across sample sites. Overall, DOC at all sites had an average mean age of 253 years B.P. (range = 0-1,900 years B.P.). Sites with the oldest mean age of DOC in May tended to have relatively aged DOC in July through September (Figure 4; r 2 = 0.68; p = 0.006; y = 1.56x + 340.1). The mean age of DOC across sites in May was 119.6 ± 230.3 (S.D.) years B.P., but in July through September it was 361.5 ± 459.0 years B.P. The shifts in Δ 14 C-DOC generally corresponded to seasonally increasing stream water SpC, consistent with active layer thaw and greater access to weathered mineral soils ( Figure  S3 in Supporting Information S1).
We observed an inverse relationship between Δ 14 C-DOC values and the relative contribution of the solute-rich endmember, burned area, and silty soil coverage ( Figure 5). These trends were generally stronger in the fall, when mean age of DOC was typically older (i.e., more depleted Δ 14 C values). Further, there was no clear relationship between the concentration of DOC and Δ 14 C-DOC ( Figure S4 in the Supporting Information S1). The difference from spring to fall appeared larger at sites with a greater relative contribution (>∼30%-50%) from the solute-rich endmember (Figure 5a). Distinct seasonal shifts in Δ 14 C-DOC were also observed across the gradient of relative catchment burn extent (Figure 5b). Among the catchments, shifts in stream Δ 14 C-DOC relative to fire extent and silt extent were nonlinear, and most pronounced in sites with >50% catchment burned area, with the most depleted Δ 14 C-DOC values (thus oldest mean age DOC) detected in the catchment with the most extensive fire disturbance (>85% catchment burned). Similarly, the most depleted Δ 14 C-DOC values were found in sites where upland silt soils cover 50% of the watershed. Seasonal shifts in Δ 14 C-DOC were also nonlinear, and most pronounced in catchments with the most extensive burned area.
The high solute endmember, area burned, and upland silt extent all had significant negative relationships with Δ 14 C-DOC ( Figure 6), with season additionally contributing significantly to each model. Interestingly, estimated   permafrost extent was unrelated to Δ 14 C-DOC (not shown). Due to collinearity between the high solute endmember, area burned, and upland silt ( Figure  S5 in Supporting Information S1), it was not appropriate to include them all in the same model, and it is not possible to tease apart the unique contributions of each to DOC isotopic composition or mean age. Instead, we present separate models for each landscape variable with season to demonstrate the shape of the relationships (Figures 6, S6-S8 in the Supporting Information S1). DOC age increased with the solute rich endmember and season (Figures 6a and 6b; GAM deviance explained = 69.5%, p < 0.001), catchment area burned and season (Figures 6c and 6d; deviance explained = 47.1%, p < 0.001), and upland silt fraction and season (Figures 6e and 6f; deviance explained = 51.1%, p < 0.001).

Discussion
Boreal regions having sporadic to discontinuous permafrost are highly susceptible to permafrost thaw, the mobilization and processing of soil OC stores, and hydrologic change (Karlsson et al., 2021;Serikova et al., 2018) owing to the presence of warm permafrost and the potential for talik development Minsley et al., 2016;Walvoord et al., 2019). Such thaw-induced changes in hydrologic flow paths have major implications for the fate and transport of aged DOC released with thaw (Striegl et al., 2005). By sampling distinct catchments in the discontinuous permafrost zone of Alaska spanning gradients in season, catchment structure, fire extent, and hydrology, we provide the spatiotemporal resolution needed to better constrain the patterns and controls on aged DOC export in headwaters (Figure 7). Values of Δ 14 C-DOC and thus the mean age of DOC ranged widely across catchments and seasons (Figures 4 and 5), consistent with a synthesis of Δ 14 C data from broadly distributed northern watersheds experiencing permafrost thaw (Estop-Aragonés et al., 2020). GAMs ( Figure 6) identified that Δ 14 C-DOC values increased with hydrologic contributions from the solute-rich endmember, which was largest in catchments dominated by upland silty soils and following fire ( Figure 7). These results build on previous studies that have inferred thawing permafrost from changing solute content and composition in major rivers of the YRB (Striegl et al., 2005;Toohey et al., 2016). The combined effects of long-term warming and seasonal active layer deepening, plus increased fire disturbance (Pastick et al., 2017) will interactively mobilize greater quantities of ancient DOC into headwaters in the future, but these responses will likely be quite heterogeneous on the landscape (Figure 7).

Accounting for Seasonal Effects on Aged DOC Export to Streams
Our repeated sampling of DOC isotopic composition, coupled with measurements of SpC and DOC concentration showed that most sites underwent a seasonal shift to a more aged DOC pool in fall (Figures 4, 5b, and 5c). This is consistent with findings in other permafrost dominated landscapes and in nonpermafrost catchments having seasonal flow (Mann et al., 2015;Neff et al., 2006). Seasonality was an important predictor of Δ 14 C-DOC (Figures 6b, 6d, and 6f), and likely accounted for some of the relationship between DOC isotopic composition and increasing interaction between water and deeper, mineral soils by fall, once soils had thawed. Many of the sites exporting aged DOC had burned in the last 20 years, but the relationship between burned extent and DOC mean age and isotopic composition is not monotonic (Figure 5b). Furthermore, two of the sites with depleted Δ 14 C-DOC values have not burned within the last 100 years . Studies have observed increasing stream solute loads with seasonal thaw (Harms & Jones, 2012;Koch, Runkel, et al., 2013;Maclean et al., 1999;Petrone et al., 2007), altered river chemistry with increasing thaw (Dornblaser & Striegl, 2015;Striegl et al., 2005;Toohey et al., 2016), and export of ancient C following years of extreme thaw (Schwab et al., 2020), but this study is one of the first to quantify a seasonal trend toward depleted Δ 14 C-DOC and thus an increase in the age of DOC. Revision to the conceptual model of mechanisms regulating aged dissolved organic carbon (DOC) delivery to stream networks of interior Alaska. (a) The null model predicts that aged DOC in this boreal landscape may be widely transported into headwater streams in a uniform way. (b) Our cross-basin seasonal survey of headwater streams refines this conceptual model by showing that aged DOC export is not uniform, and the mean age of DOC increased (and Δ 14 C-DOC values become more depleted) as a function of the interactive effects of seasonal warming and active layer deepening, the fraction of the catchment burned by fire, plus the proportion of catchment underlain by upland silty deposits. The horizontal axis of the graphs in both panels represents movement through the individual underlying model stream networks from left to right.
The release of aged soil C stocks undoubtedly had a strong effect on headwater stream DOC at several of our sites, as the depth of seasonal thaw progressed. Multiple sources could contribute to this aged carbon signature. Yedoma exists in this region due to syngenetic permafrost development in an eolian depositional environment during the late Pleistocene, leading to OC ages ranging from 10,000 to 40,000 years BP and deposits that can be decameters thick (Strauss et al., 2017). A small contribution from such sources could yield the aged DOC observed at our sites. Alternatively, larger contributions from deeper active layer soils could also drive these patterns without engagement of yedoma. At one Hess Creek location near our study sites, the age of OC in the soil active layer at 30-50 cm depth ranged from 460 to 830 years BP (n = 3), while OC age in the shallow permafrost ranged from 1,205 to 1,385 years BP (Wickland et al., 2018). Thus, although age data are sparse, our observations confirm that deeper soil layers are actively contributing DOC to headwaters in the YRB.

Aged DOC Input Greatest in Catchments With High Stream Solute Concentrations
Streamflow containing aged DOC was sourced from waters containing high concentrations of DOC and SpC (Figures 5a and 6a). High DOC and SpC in Alaska streams may serve as an indicator of thermokarsting of ice-rich soils. Elevated SpC indicates that this water has had contact with mineral soils, and near-surface permafrost tends to have higher OC content and more leachable OC than overlying active layer soils (Wickland et al., 2018). Large proportions of the high-solute endmember were limited to three locations-Richardson, Erickson, and Isom catchments ( Table 1), all of which are dominated by silt loess . Although these sites have incised stream channels and seeps flowing through gullies on the hillsides, no relationship was found between catchment slope and DOC isotopic composition (not shown). Catchments with a greater extent of silt-rich soils may be a dominant source of aged DOC and high solute loads across northern landscapes, because the easily eroded soils transfer leachable solutes to the runoff and allow for the creation of soil pipes, channel incision, and subsequently elevated interactions between runoff and ice-rich permafrost (Koch, Ewing, et al., 2013). This interpretation is supported by the collinearity between the high solute endmember and the proportion of upland silts in each catchment ( Figure 5), and the importance of upland silt as a predictor of Δ 14 C-DOC (Figure 6c). The highest proportion of the solute-rich endmember was found in two thermokarst seeps in the Richardson catchment. One of these seeps is silt laden, flowing through an actively thawing area containing thermokarst slumps and pits. The second one emanates from a steep hillslope and is slightly incised. Both seeps flow through a burned area, are deeply incised at their confluence with Richardson Tributary, and contain elevated uranium isotope concentrations that substantiates export of permafrost-derived solutes (Koch, Ewing, et al., 2013). The Isom catchment is similar in that it is dominated by silty, low-hydraulic conductivity soils , silt-laden soil pipes and contains several pingos, which provide additional evidence of preferential subsurface flow and the mobilization of thaw-derived sediments. Two of the four locations with the high solute endmember (Isom and Erickson Trib) have not burned recently, highlighting that aged DOC may be mobilized in the absence of recent fire. While catchments with the greatest upland silt fraction tended to have the most depleted Δ 14 C-DOC, there was considerable variability in isotopic values across a gradient of catchment silt content (Figure 5b). This indicates that while the variables predicting DOC isotopic composition, such as seasonality, high-solute endmember fraction, and fire extent are important in predicting headwater DOC mean age and Δ 14 C-DOC, they also reflect some degree of collinearity with soil silt content.
Differences in the mixing tendencies between tributaries and streams suggest that aged DOC can be derived from shallow subsurface flow paths rather than deep groundwater. Later in the summer (i.e., July through September), streams roughly fall along a mixing line from precipitation to the solute-rich endmember (except for the West Fork of Dall Creek, which emanates from a large, lowland wetland complex; r 2 = 0.82, p < 5E−6; Figures 2b and S3 in the Supporting Information S1), with a larger fraction of the solute-rich endmember in streams draining silty catchments. Tributaries on the other hand are either a mixture of precipitation and the organic endmember, or the organic endmember and the high-solute endmember. The lack of a high SpC, low DOC samples, or seep samples along the stream mixing line suggests that only water that has had substantial contact with organic soils can attain a high SpC signature in these systems. Although deeper groundwater tends to have elevated SpC, it does not necessarily contain high DOC concentrations, and thus deep groundwater flow alone cannot explain the presence of depleted Δ 14 C-DOC in the tributaries. Insignificant differences in δ 18 O concentrations of the endmembers ( Figure S1 in Supporting Information S1) further indicates limited influence of a unique, deep, or regional groundwater source on tributaries or streams. The combination of high DOC and high SpC may be indicative of the elevated leachability of solutes from eroding and thawing silty yedoma soils. While we did not differentiate yedoma as a specific subset of conditions of upland silty soils, the prevalence of deep thermokarst lakes within or adjacent to the Erickson, Richardson, and Isom catchments indicates that these basins are dominated by yedoma. Silty yedoma is broadly distributed across the Arctic and Alaska (Kanevskiy et al., 2011;Péwé, 1975) and has low permeability, which may lead to a substantial fraction of runoff flowing along the organic-mineral boundary , thereby picking up a high DOC concentration. Yedoma soils also tend to have ice-rich permafrost and store substantial quantities of aged C (Strauss et al., 2017) that can be exported to streams if fire or thermokarst increase erosion and subsequent percolation to the frozen boundary.
Although disturbance may lead to the release of older OC, we do not find evidence that permafrost thaw will lead to an increase in stream DOC. There was no relationship between DOC concentration and age ( Figure S4 in Supporting Information S1), and the high solute endmember had lower DOC (38.9 mg/L) than the shallow runoff through organic soils (i.e., the organic endmember, DOC = 44.2 mg/L), suggesting that OC contribution from shallower soils outweighs those from deeper, thawing permafrost. DOC yields have been shown to be greater from permafrost soils relative to active layer soils (Wickland et al., 2018), so we infer that lower DOC concentrations from permafrost soils are related to transport limitation. Hydrology plays an important role in C transport (Koch, Runkel, et al., 2013), and the largest annual flush of water occurs with the spring snowmelt when soils are still mostly frozen and deep soils cannot be accessed by runoff (Carey, 2003;Koch, Ewing, et al., 2013). Permafrost-derived DOC can generally only be mobilized later in the year after the active layer has thawed and new permafrost thaw has initiated, which also tends to coincide with lower hydrologic fluxes.

Complex Link Between Fire Disturbance Extent and Aged DOC Export
Our results demonstrate that modeling the long-term effects of fire on aged DOC export to streams is complex, and interactive with soil texture and thermokarst, seasonal active layer thaw, and resulting hydrologic flow paths (Figures 5-7). Although fires immediately mobilize massive quantities of terrestrial C into the atmosphere and adjacent headwaters, the longer-term, sustained link between catchment fire disturbance and permafrost C export varies and is dependent on fire severity and time since the fire event (Rodríguez-Cardona et al., 2020;Yoshikawa et al., 2003), which is closely correlated to vegetation type, soil moisture, and topography . Work in peatlands indicates that we may currently be underestimating the role of fires in mobilizing previously stored OC (Dean et al., 2019;Gibson et al., 2018). On the other hand, new evidence shows no connection between fluvial C exports via atmospheric emissions, and increasing fire in northern Canadian watersheds (Hutchins et al., 2020), suggesting that the export and processing of soil OC may be less impacted by fire than predicted.
Here, our detailed sampling shows a nonlinear relationship between Δ 14 C-DOC (mean DOC age) and the estimated extent of past fire disturbance of each study catchment (Figures 5c and 6b). This observation supports the notion that fires are a significant mobilization mechanism for permafrost OC in the YRB, although the GAM that included fire had less explanatory power than GAMs including the solute-rich endmember and upland silt fraction. However, the relative importance of fire may change because fire disturbances are expected to increase in the YRB linked to climate-driven warming and drying (Pastick et al., 2019) and increased lightning (Poujol et al., 2021). A greater burn extent in the future will likely enhance permafrost DOC export in YRB catchments over mid-to long-term time scales (years to decades), but these influences will be heterogeneous across catchments.

Future Projections for Disturbance and Implications for Stream C Cycling
Our comprehensive sampling provides a new level of understanding of aged DOC mobilization in northern headwaters ( Figure 7). Specifically, it indicates that aged DOC is released primarily from catchments dominated by runoff that interacts with deeper, mineral soils especially in silty landscapes, following disturbance by thermokarst and/or fire. Furthermore, this DOC export is regulated by annual thaw cycles. These factors render the mobilization of aged DOC extremely heterogeneous among headwater catchments in the discontinuous permafrost zone of Alaska, and likely elsewhere given the ubiquity of silty soils, thermokarst, and fire in permafrost landscapes.
Our interpretations are based on a mixing model driven by empirical observations of endmembers, and thus may be improved by a broader sampling effort. The high-solute endmember is the least certain, given that it is based on a small sample size and not bound by theoretical considerations. For comparison, the precipitation end member cannot have a much lower concentration. The high-solute endmember could be theoretically bound by permafrost OC concentration and yield, but this may not be relevant without understanding the interactions between water and thawing carbon that ultimately determine leachate chemistry and amount. The fact that the high solute endmember falls close to the extrapolated relationship between stream DOC and SpC ( Figure S3 in Supporting Information S1) supports the use of empirical, runoff data to define the high-solute endmember. New samples from these environments may result in small changes in the high solute endmember concentration and thus in the fraction of the high-solute endmember relative to the other two (Figure 3), but this would not impact the relative trend between catchments (Figure 3) or the trend between the endmember and Δ 14 C-DOC (Figures 5  and 6). Increasing Arctic precipitation and runoff (Rawlins et al., 2010;Wrona et al., 2016) could also impact the high-solute endmember concentration by increasing interactions between runoff and permafrost soils.
Although long-term warming and increased fire disturbance (Pastick et al., 2019) could interactively mobilize great quantities of both modern (Dornblaser & Striegl, 2015) and aged DOC with depleted Δ 14 C values into headwaters, our findings show that these responses will be heterogeneous on the landscape. In line with this conclusion, ancient C mobilization into nearby Alaskan lakes has also been shown to be heterogeneous (Bogard et al., 2019;Elder et al., 2018). The mobilization of aged C into Alaskan headwaters due to warming and fire disturbance is somewhat analogous to changes underway in temperate regions experiencing increased loading of aged, OC-rich soils into aquatic networks in response to hydrologic and land use change Drake et al., 2019Drake et al., , 2020Moore et al., 2013). When aggregated at the global scale, the mobilization of aged DOC observed in the YRB appears to be part of a broader trend toward increased cycling of ancient C in aquatic networks. This phenomenon may increase the rate of C exchange between land and the atmosphere, but we caution against broad extrapolation of observations from individual sites throughout boreal and arctic biomes, given the extreme heterogeneity of aged DOC inputs among regional watersheds.

Data Availability Statement
All original data generated in this study are available in Foks et al. (2020) (water chemistry data) and Koch et al. (2020) (continuous discharge data). Generalized Additive Model model and statistical analysis code are available at https://github.com/MattBogard/stream14C-DOC.