Thermokarst amplifies fluvial inorganic carbon cycling and export across watershed scales on the Peel Plateau, Canada

The chemical weathering of minerals is a primary control on atmospheric CO2 levels and Earth’s climate over geological timescales. As climate warming and precipitation intensify at high latitudes, glaciated terrains across the circumpolar north are poised for rapid geomorphic change and associated changes in mineral weathering dynamics. Here, we determine how the effects of permafrost thaw on mineral weathering sources and inorganic 15 carbon cycling and export integrate across watershed scales (from ~2 to 1000 km2) in a permafrost terrain within a former glacial margin and dominated by relatively inorganic sediments (Peel Plateau, Canada). Our work was conducted along three nested transects with varying intensities of retrogressive thaw slump (RTS) thermokarst activity: a 550 m transect along a first-order thaw stream within a RTS; a 14 km transect along a stream which directly received RTS inputs; and a 70 km transect along a larger stream which received inputs from RTS-affected 20 tributaries. In the thaw stream, rapid sulfuric acid weathering of carbonate tills appeared to amplify CO2 efflux to the atmosphere and HCO3– export downstream, where DIC and CO2 stable isotopes revealed a shift to an abioticinorganic driven aquatic carbon cycle. Along the intermediate transect, DIC concentrations were ten times higher in the RTS-affected reach than in the undisturbed headwaters, and decreased downstream with decreasing RTS area. Along the largest transect, HCO3– concentrations increased by two orders of magnitude in association with RTS 25 activity, despite RTSs covering only ~0.5% of the landscape. Statistical modeling of hydrochemical measurements and geospatial landscape data showed that RTSs were a primary landscape driver of HCO3– export across watershed scales. Constraining sources and rates of mineral weathering across diverse permafrost terrains will help to understand future changes in Arctic aquatic carbon cycling, as our results suggest that abiotic-inorganic processes may become prevalent. 30


Introduction
The chemical weathering of minerals is a fundamental control on atmospheric CO2 levels over geological timescales and thus Earth's long-term climate (Berner, 1999). Rapid modern warming at northern latitudes (Serreze and Barry, 2011) is driving a suite of changes including the degradation of the cryosphere (Biskaborn et al., 2019), increasing vegetation productivity (Bjorkman et al., 2018), intensifying hydrologic cycles (Rawlins et al., 2010), and strengthening land-freshwater linkages (Vonk et al., 2019;Walvoord and Kurylyk, 2016). These processes are activating large amounts of mineral substrate into biogeochemical cycles, with significant implications for the rate and magnitude of mineral weathering over modern timescales (Lacelle et al., 2019;Wadham et al., 2019). In recent decades, fluxes of carbonate alkalinity (HCO3 -, CO3 2-) and solutes in rivers across the circumpolar north have increased and more than doubled in some regions (Drake et al., 2018a;Tank et al., 2016;Toohey et al., 2016; 40 Zolkos et al., 2018). This intensification of dissolved inorganic carbon (DIC = S[CO2, HCO3 -, CO3 2-]) cycling reflects enhanced mineral weathering associated with a thickening of the active layer, deepening hydrologic flowpaths into mineral soils, and greater soil acidity associated with increasing vegetation productivity (Drake et al., 2018a;Tank et al., 2016;Toohey et al., 2016). Coupled to these processes, carbonate equilibrium reactions along the land-freshwater continuum determine the balance of DIC species and therefore its susceptibility to atmospheric 45 exchange as CO2 versus export to the ocean as alkalinity (Stumm and Morgan, 1996).
From a geochemical perspective, three coupled factors primarily influence the magnitude and directionality of mineral weathering within climate feedbacks: the mineral composition of substrate; the acid responsible for chemical weathering; and the rate of weathering (Stumm and Morgan, 1996). First, half of riverine carbonate alkalinity exported to the ocean is returned to the atmosphere over geological timescales by precipitation reactions 50 within the marine carbon cycle (Calmels et al., 2007). Hence, silicate weathering by atmospheric or soil respiration CO2 dissolved in water (carbonic acid, H2CO3) represents a long-term CO2 sink, as this process generates alkalinity without liberating mineral carbon. In contrast, carbonate weathering by H2CO3 is a CO2 sink only over shorter timescales, as half of the alkalinity produced is geogenic. HCO3produced during carbonate weathering in the presence of strong acids, for instance sulfuric acid (H2SO4) from sulfide oxidation, is thus a CO2 source. Second, 55 these reactions demonstrate how the mineral composition of a substrate and the acids responsible for its weathering can influence the degree to which weathering is a CO2 source or sink. Lastly, rates of chemical weathering are orders of magnitude faster for carbonates and sulfides than for silicates (Stumm and Morgan, 1996). Further, weathering rates generally increase with mineral surface area, and therefore are generally fast in glacial environments owing to intense physical weathering (Anderson, 2007). Indeed, hydrochemical signatures of trace 60 carbonate and sulfide lithologies can dominate weathering fluxes in primarily silicate glacial environments (Anderson, 2007). The disparity is so significant that, when sediment supplies are sufficient, H2CO3 carbonate weathering in proglacial streams can consume dissolved CO2 to below atmospheric levels (Sharp et al., 1995;St. endowed former glacial margins across northwestern Canada with relatively inorganic, ice-rich permafrost (Kokelj et al., 2017b). Today, the climate-driven renewal of post-glacial landscape change is mobilizing immense stores of minerals into modern biogeochemical cycles (Kokelj et al., 2017a). On the Peel Plateau (NWT, Canada), retrogressive thaw slumps (RTSs) expose carbonate-and sulfide-bearing glaciogenic permafrost sediments that are 75 tens of meters in thickness. The chemical weathering and fluvial transport of these sediments results in increased HCO3immediately downstream and greater solute and sediment loads throughout downstream systems Malone et al., 2013;Zolkos et al., 2018). Yet, it remains unknown how thermokarst effects on mineral weathering and DIC sources and cycling integrate across watershed scales on the Peel Plateau and in relatively inorganic-rich permafrost terrains elsewhere. These terrains, which are thought to occur within former glacial limits 80 across the northern permafrost zone (Zolkos et al., 2018), represent a frontier with respect to current knowledge on carbon cycling and climate feedbacks.
In this study we evaluated trends in major ions, DIC concentration, and dual DIC and CO2 stable isotopes along transects within three nested watersheds in the Stony Creek watershed (1100 km 2 ) on the Peel Plateau. The three transects spanned gradients of thermokarst disturbance: (i) a 550 m thaw stream formed by a runoff channel within 85 an active retrogressive thaw slump (RTS FM2 runoff); (ii) a 14 km transect in a creek which originated in undisturbed headwaters, but was directly affected by RTS FM2 and additional RTSs along the studied reach (Dempster Creek); and (iii) a 70 km transect in a large stream which received inputs from multiple large RTSaffected tributaries and was itself a major tributary of the 70,000 km 2 Peel River watershed (Stony Creek) (Fig. 1).
Our nested watershed approach and comparisons of stream chemistry between undisturbed headwaters and 90 thermokarst-affected reaches enabled us to develop a novel conceptual model which details how RTS effects on inorganic carbon cycling integrate across watershed scales. These results help to bring abiotic-inorganic aquatic processes into our conceptualization of thermokarst effects on carbon cycling, which to-date have been driven by studies outside of former glacial limits and focused on organically-driven processes (Vonk et al., 2015).

Study Area
The Stony Creek watershed is located southwest of the hamlet of Fort McPherson, in the northern, or lower Peel River watershed (Fig. 1). Stony Creek, a tributary of the Peel River, originates in the Richardson Mountains, where the slopes are sparsely vegetated and mainly consist of colluvium from exposed marine shale and sandstone bedrock (Norris, 1985). As Stony Creek flows eastward, the main channel and its tributaries incise ice-rich glacial deposits 100 and underlying Cretaceous bedrock, forming a stream network comprised of tundra flow tracks grading to incised gravel bed streams with increasing watershed size. The fluvially-incised valleys and increasing regional precipitation have proven conducive to thaw-driven mass wasting of ice-rich glacial deposits and formation of RTSs (Kokelj et al., 2017b). RTS growth is driven by the ablation of exposed ground ice and is perpetuated by the downslope mass wasting of thawed material via fluidized earth flows, which can accumulate large volumes of debris 105 in stream valleys (Fig. 1). Across the Stony Creek watershed, intensifying RTS activity releases large volumes of https://doi.org/10.5194/bg-2020-111 Preprint. Discussion started: 15 April 2020 c Author(s) 2020. CC BY 4.0 License. sediment and solutes into streams relative to undisturbed headwaters (Kokelj et al., 2017b;Segal et al., 2016). This substrate is transported to streams via rill runoff channels in the scar zone and debris tongue deposits in the stream valley. Impacts to Stony Creek are representative of numerous other major Peel River tributaries that have incised the ice-rich Peel Plateau (Kokelj et al., 2015). The ~60 km 2 watershed of Dempster Creek, a tributary of Stony 110 Creek, originates in willow and open spruce taiga without RTS activity, receiving large inputs of sediments and solutes from RTSs FM2 and FM3 within several kilometers of the headwaters Malone et al., 2013). Many small, non-RTS affected streams and several larger RTS-affected tributaries flow into Dempster Creek before its confluence with Stony Creek.

115
In late July, 2017, we sampled three nested transects: (i) at five locations along a 550 m-long transect within the RTS FM2 rill runoff channel (Fig. 1a); (ii) Dempster Creek in the undisturbed headwaters, along its RTS-affected reach upstream of three major tributaries, and from the tributaries near their mouths, to characterize the downstream effects of RTS FM2 on DIC; (iii) Stony Creek in the undisturbed headwaters, along the RTS-affected reach of Stony Creek upstream of seven major tributaries, and from the tributaries near their mouths, to determine RTS effects on 120 DIC across broader scales (Fig. 1, Table A1). While the Dempster and Stony Creek headwater sites were not affected by RTSs, mainstem and tributary sampling sites were (Table A1). Minor tributaries of Dempster and Stony Creek, many of which were not affected by RTSs, were not sampled.
At all sampling sites, stream temperature and pH were measured using a pre-calibrated YSI Professional-Plus water quality meter. At most sites, water samples were collected for ions, DIC, CO2, CH4, dissolved organic carbon 125 (DOC), UV-visible absorbance, and total suspended sediments (TSS). Along the RTS FM2 runoff transect, we sampled only for DIC and CO2 concentration, and stable isotopes of dissolved CO2 (d 13 CCO2). One day prior, additional parameters were sampled at runoff site five, located near the confluence of the RTS runoff with Dempster Creek, for comparison with the full suite of chemistry parameters collected along the Dempster Creek transect. At the Dempster and Stony Creek sites, we additionally sampled water for stable isotopes of DIC (d 13 CDIC) and d 13 CCO2 130 and used dual DIC and CO2 stable isotopes to characterize abiotic and biotic processes influencing DIC sources and cycling.
Water samples were collected from the thalweg where possible, as an integrated sample from ~15 cm below the surface to ~1 m depth. An additional sample for TSS was collected in a 1 L HDPE in the same fashion. Water samples were filtered using sample-rinsed 0.45 µm polyethersulfone (PES, ThermoFisher) or cellulose-acetate (CA, 135 Sartorius) membranes. Samples for DIC were collected without headspace in airtight syringes. Samples for ions, DOC, and UV-visible absorbance were collected in acid washed (24 h, 10% v/v HCl) all-plastic syringes. Syringes were triple sample-rinsed, sealed without headspace, and stored cool and dark until processing within 10 h. Water for DIC was filtered (PES) into precombusted (5 h, 500˚C) glass vials without headspace and sealed with a butyl rubber septum for DIC concentration or two silicone-teflon septa for d 13 CDIC. Samples for cations were filtered (CA) 140 https://doi.org/10.5194/bg-2020-111 Preprint. Discussion started: 15 April 2020 c Author(s) 2020. CC BY 4.0 License. into acid-washed bottles and acidified with trace metal-grade HNO3, while anions were filtered (CA) but not acidified. Samples for DOC were filtered (PES) into precombusted glass vials and acidified to pH < 2 using trace metal-grade HCl (Vonk et al., 2015). Samples for UV-visible absorbance were filtered (PES) into non-acid washed 30 mL HDPE bottles. Samples were refrigerated (4˚C, dark) until analysis.
Dissolved gases were collected following the headspace equilibration method (Hesslein et al., 1991) and stored in 145 airtight syringes (for CO2 concentration) or over-pressurized in pre-evacuated serum bottles sealed with pre-baked (60˚C, 12 h), gas-inert butyl rubber stoppers (for d 13 CCO2, CH4). At each site, atmospheric samples for CO2 and CH4 concentration and d 13 CCO2 were stored in the same fashion. Gas samples were stored in the dark at ~20˚C prior to analysis within 10 h (CO2) or two months (d 13 CCO2). Water and air temperature, atmospheric pressure, and the volumetric ratio of sample to atmospheric headspace was recorded for correcting later calculations of CO2 partial 150 pressure (pCO2) and d 13 CCO2 (Hamilton and Ostrom, 2007).

Stream Flow
Discharge (Q) in Stony Creek tributaries was estimated from a hydraulic geometry model (Gordon, 2004) that we ) and cross-sectional area made at increments equal to 10% of stream width (Gordon, 2004;Lurry and Kolbe, 2000), and were averaged for sites with multiple measurements.

Geochemical Analyses
Upon returning from the field each day, CO2 was measured using an infrared gas analyzer (PP Systems EGM-4).
The EGM-4 was checked monthly for drift using a commercial standard (Scotty Gases). We calculated pCO2 using Henry's constants corrected for stream water temperature (Weiss, 1974) and accounting for the ratio of water volume to headspace during sample equilibration. DIC samples were measured by infrared CO2 detection (LiCOR 165 7000) following acidification within a DIC analyzer (Apollo SciTech model AS-C3). Calibration curves were made daily using certified reference material (CRM) from Scripps Institution of Oceanography. Samples with DIC concentrations < 400 µM were analyzed using solutions prepared from a 1000 ppm TIC stock (ACCUSPEC) that were calibrated with CRM. DIC species (CO2, HCO3 -, CO3 2-) were calculated from DIC concentration and pCO2 or pH using CO2sys (v.2.3) (Pierrot et al., 2006), using field temperature and pressure at the time of sampling, and the 170 freshwater equilibrium constants for K1 and K2 (Millero, 1979). https://doi.org/10.5194/bg-2020-111 Preprint. Discussion started: 15 April 2020 c Author(s) 2020. CC BY 4.0 License.
Cations and trace elements were measured by optical emission spectrometry (Thermo ICAP-6300) and anions by ion chromatography (Dionex DX-600) at the University of Alberta Biogeochemical Analytical Services Laboratory (BASL, ISO/EIC accreditation #17025). DOC was measured using a total organic carbon analyzer (Shimadzu 5000A). DOC standard curves were made daily with a 1000 ppm KHP solution (ACCUSPEC) and an in-house 175 caffeine standard (10 mg L -1 ) was run every 20 samples to monitor instrument drift. Absorbance spectra were analyzed using an Ocean Optics UV-VIS instrument with a Flame spectrometer module, following Stubbins et al. (2017) and corrected for Fe interference (Poulin et al., 2014). To evaluate organic carbon reactivity, we used specific ultraviolet absorbance at 254 nm (SUVA254) to infer DOC aromaticity (Weishaar et al., 2003). d 13 CDIC was determined using an isotope ratio mass spectrometer (Finnigan Mat DeltaPlusXP) interfaced to a total 180 organic carbon analyzer (OI Analytical Aurora Model 1030W) at the University of Ottawa Stable Isotope Laboratory. d 13 CCO2 and CH4 concentration were analyzed within two months using a Picarro isotope analyzer (G2201-i; < 0.2‰ precision, CH4 operational range = 1.8-1500 ppm) equipped with an injection module for discrete samples (SSIM). Commercial d 13 CCO2 and CH4 standards were used to check for drift during each run. We used mass balance to correct d 13 CCO2 values for the d 13 C and mass of atmospheric CO2 used for equilibration (Hamilton 185 and Ostrom, 2007). To assess d 13 CCO2 fractionation from calcite precipitation (Turner, 1982) and methanogenesis (Campeau et al., 2018) in RTS FM2 runoff, we calculated the saturation index (SI) and partial pressure of CH4 (pCH4). SI was calculated using the hydrochemical software Aqion version 6.7.0 (http://www.aqion.de), which uses the U.S. Geological Survey software PHREEQC (Parkhurst and Appelo, 2013) as the internal numerical solver.
Samples for atmospheric and dissolved CH4 were collected in the same fashion as d 13 CCO2. pCH4 was calculated 190 using Bunsen solubility coefficients (Wiesenburg and Guinasso, 1979) converted to the appropriate units (Sander, 2015).
TSS samples were filtered onto muffled (450˚C, 4 h) and pre-weighed glass fiber filters (Whatman GF/F; 0.7 µm nominal pore size) upon returning from the field, stored frozen, and dried (60˚C, 24 h) for gravimetric analysis following a modified version of U.S. Geological Survey Method I-3765.

Mineral Weathering and DIC Sources
We used a Piper diagram, which reflects the proportional equivalent concentrations of ions in a sample relative to mineral weathering end-members, to constrain the sources of mineral weathering and HCO3 -. The products of Eq. A1-A7 defined the mineral weathering end-members in the Piper diagram (Table A2). We further constrained mineral weathering and DIC sources using d 13 CDIC and pH. End-member δ 13 CDIC ranges for equilibrium processes 200 (mixing with atmospheric and/or biotic CO2) and kinetic reactions (mineral weathering) were derived following Lehn et al. (2017).
To evaluate potential effects on d 13 CCO2 from DIC speciation along the pH continuum (Eq. A8) (Clark and Fritz, 1997), we compared theoretical and observed d 13 CCO2 values in the Stony Creek mainstem. Theoretical d 13 CCO2 https://doi.org/10.5194/bg-2020-111 Preprint. Discussion started: 15 April 2020 c Author(s) 2020. CC BY 4.0 License. values were calculated using mass balance to obtain d 13 CHCO3 from measurements of DIC, CO2, HCO3 -, d 13 CDIC, and 205 d 13 CCO2. We then used measurements of stream temperature (T) to calculate the equilibrium fractionation between CO2 and HCO3 -(ε = -9.483 ✕ 10 3 /T + 23.89‰; Mook et al., 1974). Finally, ε was subtracted from d 13 CHCO3 to obtain theoretical d 13 CCO2. Similarity between observed and theoretical d 13 CCO2 values was interpreted as d 13 CCO2 variability driven by carbonate equilibrium reactions, whereas dissimilarity was taken to reflect effects from CO2 degassing (Zhang et al., 1995) and/or biotic CO2 production (Kendall et al., 2014). RTSs were interpreted as active where exposed sediment visibly dominated the feature surface (Cray and Pollard, 215 2015) in orthorectified SPOT multispectral imagery that we pan-sharpened to 1.6 m resolution using the ArcGIS Image Analysis tool. The satellite imagery was collected from September 9 to 25, 2016. Active RTSs that were connected to streams were manually delineated using ArcGIS. We used RivEx 10.25 software (Hornby, 2017) to quantify the number of active RTSs impacting streams in the Stony Creek watershed and to visualize RTS impact accumulation across the fluvial network. We defined RTS impact accumulation as the cumulative number of active 220 RTSs impacting upstream reaches. RTSs were interpreted to impact streams based on contact or interpreted downslope flow based on slope direction and gradient from the CDEM (Supplementary Information). Where a single RTS affected multiple streams, only the upstream segment was used for the accumulation.
We used the Geomorphic and Gradients Metrics Toolbox (Evans et al., 2014) to calculate terrain roughness, which is a measure of variance across a land surface and represents topographic complexity (Riley et al., 1999). We use 225 terrain roughness as a proxy for potential physical erosion, which is known to enhance sulfide oxidation by exposing shale regolith throughout the Peel River watershed (Calmels et al., 2007) and may therefore influence DIC. The enhanced vegetation index (EVI) was used to broadly evaluate vegetation productivity (Huete et al., 2002), which is known to influence DIC production by enhancing mineral weathering (Berner, 1999). We used the U.S. National Aeronautics and Space Administration EVI product (Didan, 2015), which is derived from gridded (250 m) moderate 230 resolution imaging spectroradiometer (MODIS). The MODIS data were collected on July 28, 2017. The ArcGIS Zonal Statistics tool was used to calculate total RTS area, mean terrain roughness, and mean EVI in Stony Creek tributary watersheds.

Statistics
We tested for downstream change in HCO3concentration and pCO2 along the Stony Creek mainstem using the non-235 parametric Mann-Kendall test from the R software (R Core Team, 2018) package zyp (Bronaugh and Werner, 2013), following the trend pre-whitening approach detailed by Yue et al. (2002) to account for serial autocorrelation. We developed a multiple linear regression model to evaluate the influence of RTS activity on HCO3export in Stony https://doi.org/10.5194/bg-2020-111 Preprint. Discussion started: 15 April 2020 c Author(s) 2020. CC BY 4.0 License.
Creek tributary watersheds relative to hydrology, terrain roughness, and vegetation productivity. To account for potential effects of varying tributary watershed areas on discharge (Q) and constituent concentration, we used 240 tributary HCO3yields in the model. Instantaneous discharge (Q, m 3 s -1 ) was estimated from the hydraulic geometry relationship between Q and stream width (Eq. 1). Discharge and HCO3flux (concentration*Q) were normalized to the respective tributary watershed area and scaled to estimate daily water yield (cm d -1 ) and HCO3yield (µmol m -2 d -1 ). Daily HCO3yields in Stony Creek tributaries were modeled as: HCO3yield = RTSn + RTSarea + Water yield + TR + EVI (2)   245 where RTSn is the number of active RTSs; RTSarea is the watershed area disturbed by RTSs (%); TR is the mean terrain roughness (m); and EVI is the mean enhanced vegetation index (-1 to 1). The multiple linear regression was trimmed using the step function in the R package lmerTest (Kuznetsova et al., 2018) to eliminate covariates which did not improve the model. Highly collinear covariates were identified using a Variance Inflation Factor > 3 (Zuur et al., 2010) and removed from the trimmed models. Model fits were inspected visually with residual plots and 250 covariates were transformed as needed to meet assumptions of independent and homoscedastic residuals (Zuur, 2009). Statistics were performed in the R programming environment (v.3.4; R Core Team, 2018) and significance was interpreted at α = 0.05. To understand potential effects from variable rainfall on water yields prior to and during the two-day sampling window of the Stony Creek tributaries, we inspected total rainfall in 24 h increments preceding the sampling of each Stony Creek tributary. Hourly rainfall data were obtained from a Government of 255 Northwest Territories Total meteorological station located ~1 km from the RTS FM2 (Fig. 1).

Rapid Inorganic Carbon Cycling within a Permafrost Thaw Stream
Hydrochemical trends in the FM2 runoff transect were distinct compared to the Dempster and Stony Creek and showed that RTS activity enhanced chemical weathering, CO2 production, and fluvial export of HCO3 -.

260
Concentrations of HCO3were up to two orders of magnitude greater in the runoff of RTS FM2 (1510 µM) than in the Dempster and Stony Creek headwaters (1-115 µM), where RTSs were absent (Table 1, Fig. 2). Building on previous findings Malone et al., 2013;Zolkos et al., 2018), the increase in conductivity along the RTS FM2 runoff transect (from 1370 to 1986 µS cm -1 ) demonstrates that, upon thaw, minerals within regional permafrost tills readily weather during fluvial transport. Mineral weathering along the RTS FM2 runoff transect 265 occurred in concert with a decrease in HCO3from 1510 to 1386 µM (Table 1). The high proportion of SO4 2relative to other ions (Fig. 3) and 13 C-enriched d 13 CDIC (-1.0‰) (Fig. 4) indicates that DIC was sourced from carbonate weathering by H2SO4 from sulfide oxidation (Clark and Fritz, 1997) and also suggests that sulfate salt (e.g. gypsum) dissolution added some SO4 2-. Together, these trends reflect a strong coupling between mineral weathering and DIC production following thaw and exposure of tills by thermokarst (Zolkos et al., 2018). Further, 270 while HCO3concentrations within the RTS FM2 runoff were high relative to undisturbed headwaters, the decrease https://doi.org/10.5194/bg-2020-111 Preprint. Discussion started: 15 April 2020 c Author(s) 2020. CC BY 4.0 License. in HCO3along the transect indicates that carbonate equilibrium reactions within RTSs on the Peel Plateau can rapidly transform HCO3to CO2 during fluvial transport.
The trends in HCO3 -, pCO2, and stable CO2 isotopes along the RTS FM2 runoff transect indicated that H2SO4 from sulfide oxidation contributed to CO2 production in the upper reach near the RTS and that degassing along the lower 275 reach resulted in CO2 efflux to the atmosphere. Within the first 220 m, pCO2 increased from 1046 to 1534 µatm and relatively 13 C-enriched d 13 CCO2 (range = -11.4 to -12.1‰) (Table 1, Fig. 5) aligned with values expected from H2SO4 weathering of regional carbonate (-0.7 to -5.6‰) (Hitchon and Krouse, 1972), when accounting for isotopic fractionation of ~8‰ between carbonate and CO2 at the temperature of FM2 runoff (18˚C) (Clark and Fritz, 1997).
These trends are consistent with some abiotic CO2 production via carbonate weathering and HCO3conversion to 280 CO2. After 220 m, pCO2 decreased from 1534 to 742 µatm and d 13 CCO2 values increased from -11.4 to -10‰, reflecting the preferential loss of 12 C in the CO2 phase via DIC fractionation and degassing (Doctor et al., 2008;Drake et al., 2018b;Kendall et al., 2014). 13 C enrichment of the CO2 pool by methanogenesis (Campeau et al., 2018), photosynthesis (Descolas-Gros and Fontungne, 1990), and/or calcite precipitation (Turner, 1982) was unlikely, as CH4 in FM2 runoff was relatively low (pCH4 = 3.6 ± 1.9 µatm, mean ± SD, n = 6), the high turbidity of 285 FM2 runoff likely inhibited photosynthesis, and calcite was below saturation (SI = -0.79). These observations of rapid CO2 production along the upper reach of the transect and efflux along the lower reach align with recent estimates of high rates of CO2 efflux within RTS runoff (Zolkos et al., 2019).
Our observations of downstream DIC transformation coupled to mineral weathering indicate that recently thawed permafrost substrate can rapidly generate substantial HCO3and some abiotic CO2 during fluvial transport within 290 RTS rill runoff. A majority of the DIC was exported to larger streams as HCO3 -, while a smaller proportion was converted to CO2 and degassed during fluvial transport within the RTS runoff. Similar to research on organic matter decomposition in Yedoma regions (Drake et al., 2018b), these results from a relatively inorganic-rich terrain demonstrate that thermokarst can promote rapid inorganic carbon cycling upon thaw and during fluvial transport in meltwater thaw streams.

Abiotic-Inorganic Carbon Cycling in a Thermokarst-Affected Intermediate Watershed
Along the Dempster Creek transect, trends in HCO3concentration, pCO2, and stable isotopes reflected a clear shift in carbon cycling driven by biotic-organic processes in the undisturbed headwaters to abiotic-inorganic processes in the RTS-affected reach (Fig. 1). In the undisturbed Dempster Creek headwaters (site one), relatively low pH (5.82) (Table 1), high pCO2 (2467 µatm) (Fig. 2b), and low d 13 CDIC and d 13 CCO2 values (Figs. 4,5) reflected DIC that was 300 primarily sourced from inputs of soil CO2 and mixing with atmospheric CO2 (Campeau et al., 2018;Kendall et al., 2014). The concentration of HCO3was ten times lower in the undisturbed headwaters (115 µM) than at the downstream site nearest the FM2 runoff (site two, 1321 µM) (Table 1, Fig. 2a), where d 13 CDIC reflected H2SO4 carbonate weathering and inputs from the RTS runoff described above (Fig. 4). In conjunction with the increase in HCO3between Dempster Creek headwaters and the first RTS-affected sampling site, pCO2 decreased (from 2467 to 305 https://doi.org/10.5194/bg-2020-111 Preprint. Discussion started: 15 April 2020 c Author(s) 2020. CC BY 4.0 License. 686 µatm) and d 13 CCO2 increased (from -21.6 to -16.0‰) ( Table 1). Although DOC from RTS runoff is known to be relatively biolabile within streams on the Peel Plateau (Littlefair et al., 2017), these trends between the undisturbed headwaters and the first RTS-affected sampling site suggest stronger effects on DIC from CO2 degassing to the atmosphere than from biotic CO2 production associated with inputs of permafrost organic carbon from RTS activity (Doctor et al., 2008;Drake et al., 2018b;Kendall et al., 2014). This may be partly due to the protection of DOC 310 from microbial oxidation via adsorption to RTS sediments (Gentsch et al., 2015), which aligns with our observations of high TSS (11800 mg L -1 ) and a decrease in DOC concentration downstream of RTS FM2 (Table 1, Fig. 6) (Littlefair et al., 2017). These trends are consistent with observations that substantial CO2 is lost in headwaters via efflux to the atmosphere (Drake et al., 2018b) and also show that abiotic-inorganic processes can dominate carbon cycling where thermokarst releases inorganic substrate into fluvial networks.

315
The hydrochemical trends along the RTS-affected reach of the Dempster Creek transect (sites 2-4) showed that the balance between inputs from RTS-affected and undisturbed streams moderated the magnitude of HCO3export within downstream environments. The decrease in HCO3 -(1321 to 946 µM; Fig. 2) and conductivity (958 to 416 µS cm -1 ; Table 1) along the entire RTS-affected reach of the transect occurred in conjunction with inputs from dozens of small, non RTS-affected tributary streams, suggesting that these streams were partly responsible for diluting the 320 inputs from RTS activity. Further, the area of RTS disturbance within successive downstream watersheds decreased from 3.2 to 1.2% (
Together, our observations along the RTS FM2 runoff (Sect. 3.1) and Dempster Creek (Sect. 3.2) transects suggest that rates of carbonate weathering and CO2 degassing are most pronounced following thaw and within headwaters.
These effects are less pronounced in mid-order streams like Dempster Creek, which mix with inputs from https://doi.org/10.5194/bg-2020-111 Preprint. Discussion started: 15 April 2020 c Author(s) 2020. CC BY 4.0 License. undisturbed tributaries, may receive and/or generate some biotic CO2, and export HCO3downstream at a magnitude 340 coupled to the area of RTS disturbance.

Thermokarst Effects on HCO3 -Export and CO2 Production across a Major Watershed
Hydrochemical trends along the larger Stony Creek transect resembled those in Dempster Creek and showed that RTS activity strengthens signals of H2SO4 carbonate weathering and amplifies HCO3export across watershed scales. In the mountainous, non-RTS affected headwaters of Stony Creek, where bedrock lithologies contain limited 345 carbonate, the Ca 2+ -Mg 2+ -SO4 2--type waters at site one (Fig. 3a) reflected the denudation of sulfide-and gypsumbearing bedrock exposed in the Richardson Mountains (Norris, 1985). In these headwaters, pCO2 near atmospheric equilibrium (Fig. 2), low HCO3 -(33 µM) and pH (5.66), and intermediate d 13 CDIC (-11.6‰) values indicated that mixing with atmospheric and biotic CO2 were the primary sources of DIC (Fig. 4). Along the entire Stony Creek transect, the striking increase in HCO3 -(from 33 to 461 µM) (Tables 1, A1, Fig. 2a) together with ionic and d 13 CDIC 350 signals of intensified carbonate weathering (Figs. 3a, 4) reflect amplified inorganic carbon cycling associated with inputs from RTS-affected tributaries. Along Stony Creek, the RTS disturbance area increased from 0% in the undisturbed headwaters to 0.36 ± 0.03% (mean ± SD, n = 5) of the watersheds of sites four through eight (Table 1).
This builds on the findings that solutes and sediments from RTSs cascade through fluvial networks Malone et al., 2013) and shows that RTS activity, despite covering less than 0.5% of the greater Stony Creek 355 watershed, greatly increases fluvial HCO3export to the Peel River.
While the concentration of HCO3increased downstream along the entire Stony Creek transect, pCO2 was slightly above atmospheric and relatively constant (479 ± 35 µatm, mean ± SD), suggesting a greater balance between CO2 production and degassing to the atmosphere than in the RTS FM2 runoff or undisturbed headwaters. Yet, d 13 CCO2 values decreased downstream from -13.8 to -18.1‰ along Stony Creek (Fig. 5), indicating a shift in the processes 360 influencing CO2. To evaluate the relative influence of biotic CO2 production, CO2 conversion to HCO3 -, and CO2 degassing to the atmosphere on d 13 CCO2 values, we compared our measurements of d 13 CCO2 with theoretical values reflecting DIC controlled by speciation along the pH continuum (Eq. A8) (Sect. 2.5). From ~5 to 35 km downstream, the good agreement between measured and theoretical d 13 CCO2 values reflected equilibrium fractionation (ε = 9.7‰ at 9˚C; Mook et al., 1974) between CO2 and HCO3 - (Fig. 7). This suggests that carbonate 365 equilibrium reactions coupled with increasing HCO3from RTS-affected tributaries was a primary control on CO2 along the upper reach of the Stony Creek transect. Along the lower reach of the transect, as DOC availability increased (Fig. 6), d 13 CCO2 values were lower than would be expected from degassing and/or CO2 conversion to HCO3 - (Fig. 7). This suggests that biotic CO2 sustained the elevated pCO2 downstream (Fig. 3b). Yet, this signal of biotic CO2 production was not evident in the 13 C-enriched d 13 CDIC along Stony Creek (Table 1, Fig. 5), which 370 reflected a geogenic DIC pool dominated by inputs of HCO3from RTS-affected tributaries. These trends show that, in environments with abundant carbonate weathering, measuring both d 13 CCO2 and d 13 CDIC can provide a more complete understanding of DIC sources and cycling (Horgby et al., 2019). From a carbon cycling perspective, these https://doi.org/10.5194/bg-2020-111 Preprint. Discussion started: 15 April 2020 c Author(s) 2020. CC BY 4.0 License. trends in Stony Creek show that HCO3export across watershed scales is among the most striking legacies of regional thermokarst activity.

Landscape Controls on Inorganic Carbon Cycling in a Thermokarst-Affected Fluvial Network
In RTS-affected tributaries of Stony Creek (Fig. 1), HCO3yields (1070 ± 985 µmol m -2 d -1 , mean ± SD, n = 6) were several orders of magnitude higher than in the non-RTS affected headwaters (1.6 µmol m -2 d -1 ) ( Table 2) and approached summertime HCO3yields in watersheds with carbonate rock weathering by glacial activity (2700 to 3300 µmol m -2 d -1 ) (Lafrenière and Sharp, 2005;Striegl et al., 2007). To determine the importance of RTS activity 380 on HCO3export in Stony Creek tributary watersheds, we paired estimates of HCO3and water yields with geospatial data for active RTSs, terrain roughness, and vegetation productivity in a multiple linear regression model (Sect. 2.7). Water yield, the area of RTS disturbance, and terrain roughness were retained during automated covariate selection for the final model (F3 = 45, p < 0.01, R 2 = 0.96). In addition to water yield, HCO3yield was positively correlated with RTS disturbance area (Table 2). Terrain roughness was anticorrelated with HCO3yield, 385 likely owing to relatively more limited carbonate weathering in the mountainous headwaters of the Stony Creek watershed, where exposed shale bedrock predominates (Norris, 1985) (Table 1). Low HCO3 -: water yield in the Stony Creek headwaters ( Table 2) further suggested that HCO3export was limited by carbonate availability rather than by water. Higher HCO3 -: water yields in RTS-affected tributaries aligned with the model results indicating that RTS activity increases inorganic carbon availability and export across watershed scales on the Peel Plateau.

390
These findings align with the observation that climate warming and intensifying precipitation has rejuvenated deglaciation-phase geomorphic and associated mineral weathering dynamics across the western Canadian Arctic (Kokelj et al., 2017a). Mineral weathering will likely come to dominate biogeochemical cycles in these regions, which may bear increasing hydrochemical resemblance to landscapes with glacial coverage (Anderson, 2007;St. Pierre et al., 2019;Striegl et al., 2007;Wadham et al., 2019). Further, the considerable HCO3export within the 395 Stony Creek watershed, where RTS spanned < 0.5% of the watershed, demonstrates that thermokarst development across a fractional proportion of ice-and sediment-rich landscapes can substantially intensify inorganic carbon cycling. Our observations, which capture the effects of mid-summertime thaw and relatively moderate rainfall within the four days prior to sampling (Fig. A1), reflect hydrochemical trends that are likely to vary during the summertime. Constraining these trends across a broader range of hydrologic and permafrost thaw conditions will 400 help to understand variability in mineral weathering effects on inorganic carbon cycling.

Inorganic Carbon Cycling in Thermokarst-Affected Fluvial Networks
Our study is among the first to evaluate inorganic carbon cycling across gradients of thermokarst disturbance and nested watersheds, enabling an assessment of how catchment chemical characteristics and thermokarst effects integrate across watershed scales. In this framework, our findings support a novel conceptual model of land-405 freshwater linkages and carbon cycling from headwaters to intermediate scales in thermokarst terrains (Fig. 8).
Broadly, these findings reveal fast DIC cycling in headwaters and substantial DIC export downstream. Thermokarst https://doi.org/10.5194/bg-2020-111 Preprint. Discussion started: 15 April 2020 c Author(s) 2020. CC BY 4.0 License. activity increased CO2 production in headwaters and fluvial HCO3export across scales by unearthing large amounts of reactive inorganic substrate previously sequestered in permafrost (see also Kokelj et al., 2013;Malone et al., 2013). The striking carbonate weathering and DIC production that we documented in the RTS FM2 runoff indicates 410 that rapid carbon cycling can be expected where thermokarst releases inorganic substrate with limited prior modification (Lacelle et al., 2019). In the larger Stony Creek watershed, DIC cycling was characterized by relatively slower processes associated with HCO3export, which reached magnitudes comparable to watersheds with carbonate denudation by glacial activity (Lafrenière and Sharp, 2005;Striegl et al., 2007). Our findings directly link accelerating thermokarst activity on the Peel Plateau (Segal et al., 2016) with signals of intensifying carbon 415 weathering across the broader Peel and Mackenzie River watersheds (Tank et al., 2016;Zolkos et al., 2018) and suggest that accelerating thermokarst activity will intensify inorganic carbon cycling across broad swaths of the circumpolar north (Zolkos et al., 2018).
Permafrost terrains susceptible to hillslope thermokarst like RTSs are abundant within and outside of former glacial limits across the circumpolar north (Zolkos et al., 2018). Owing to regional variability in geology, glacial activity, 420 climate, and ecosystem history, the mineral composition of permafrost in these regions is likely heterogeneous (Lacelle et al., 2019;Zolkos et al., 2018). The degree to which carbonate weathering is coupled with sulfide oxidation in thermokarst terrains will determine if mineral weathering is a CO2 sink or source, and thus the strength of abiotic components of the permafrost carbon-climate feedback (Schuur et al., 2015;Zolkos et al., 2018). Current models suggest that mineral weathering across the Peel River watershed is a modest source of CO2 to the 425 atmosphere, owing to widespread sulfide oxidation coupled with carbonate weathering (Beaulieu et al., 2011). Yet these models do not account for enhanced weathering associated with RTS activity (Tank et al., 2016) and are based on relatively coarse lithologic data which under-represent the abundance of carbonate in exposed bedrock in the southern Richardson and Mackenzie Mountains (Norris, 1985). In the Eurasian north, increasing alkalinity fluxes in recent decades are thought to reflect enhanced H2CO3 carbonate weathering coincident with a suite of environmental 430 changes including intensifying permafrost thaw, precipitation, and vegetation activity (Drake et al., 2018b).
Integrating thermokarst activity and the mineral composition of permafrost into the next generation of dynamicnumerical biogeochemical models (Beaulieu et al., 2011) will facilitate a better understanding of the effects of permafrost thaw on carbon cycling across rapidly changing northern environments.

435
Few studies have traced the effects of thermokarst on carbon cycling across watershed scales or within a relatively inorganic-rich permafrost terrain. Here, we find that the climate-driven renewal of deglaciation-phase geomorphic and associated mineral weathering dynamics in the western Canadian Arctic is amplifying CO2 release in headwaters and exporting substantial HCO3primarily from H2SO4 carbonate weathering across watershed scales, despite RTS disturbance across only a fractional proportion of the landscape. Mineral weathering in our study region 440 is both a contemporary and long-term positive feedback to climate change, albeit minor at the scale of inquiry https://doi.org/10.5194/bg-2020-111 Preprint. Discussion started: 15 April 2020 c Author(s) 2020. CC BY 4.0 License. (Zolkos et al., 2019). Constraining DIC sources and fluvial export across diverse permafrost terrains is critical for understanding a rapidly changing arctic carbon cycle, as our results demonstrate that accelerating thermokarst activity may cause abiotic-inorganic processes to dominate aquatic carbon cycling across broad swaths of the circumpolar north.

445
Data availability. All data used in this study are available in the supplement.

Supplement. The supplement for this article is available online.
Author contributions. SZ and SET designed the study with contribution from RGS and SVK. SZ led the field research, laboratory analyses, and manuscript writing. JK contributed to geospatial analyses. CEA contributed to 450 laboratory analyses. All authors contributed to manuscript writing.
Competing interests. The authors declare no conflicts of interest.