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

As climate warming and precipitation increase at high latitudes, permafrost terrains across the circumpolar north are poised for intensified geomorphic activity and sediment mobilization that are expected to persist for millennia. In previously glaciated permafrost terrain, ice-rich deposits are associated with large stores of reactive mineral substrate. Over geological timescales, chemical weathering moderates atmospheric CO₂ levels, raising the prospect that mass wasting driven by terrain consolidation following thaw (thermokarst) may enhance weathering of permafrost sediments and thus climate feedbacks. The nature of these feedbacks depends upon the mineral composition of sediments (weathering sources) and the balance between atmospheric exchange of CO₂ vs. fluvial export of carbonate alkalinity (Σ[${\mathrm{HCO}}_{\mathrm{3}}^{-}$, ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$]). Working in the fluvially incised, ice-rich glacial deposits of the Peel Plateau in northwestern Canada, we determine the effects of slope thermokarst in the form of retrogressive thaw slump (RTS) activity on mineral weathering sources, CO₂ dynamics, and carbonate alkalinity export and how these effects integrate across watershed scales (∼ 2 to 1000 km²). We worked along three transects in nested watersheds with varying connectivity to RTS activity: a 550 m transect along a first-order thaw stream within a large RTS, a 14 km transect along a stream which directly received inputs from several RTSs, and a 70 km transect along a larger stream with headwaters that lay outside of RTS influence. In undisturbed headwaters, stream chemistry reflected CO₂ from soil respiration processes and atmospheric exchange. Within the RTS, rapid sulfuric acid carbonate weathering, prompted by the exposure of sulfide- and carbonate-bearing tills, appeared to increase fluvial CO₂ efflux to the atmosphere and propagate carbonate alkalinity across watershed scales. Despite covering less than 1 % of the landscape, RTS activity drove carbonate alkalinity to increase by 2 orders of magnitude along the largest transect. Amplified export of carbonate alkalinity together with isotopic signals of shifting DIC and CO₂ sources along the downstream transects highlights the dynamic nature of carbon cycling that may typify glaciated permafrost watersheds subject to intensification of hillslope thermokarst. The balance between CO₂ drawdown in regions where carbonic acid weathering predominates and CO₂ release in regions where sulfides are more prevalent will determine the biogeochemical legacy of thermokarst and enhanced weathering in northern permafrost terrains. Effects of RTSs on carbon cycling can be expected to persist for millennia, indicating a need for their integration into predictions of weathering–carbon–climate feedbacks among thermokarst terrains.

1 Introduction

Riverine export of carbonate alkalinity (Σ[${\mathrm{HCO}}_{\mathrm{3}}^{-}$, ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$]), generated by the chemical weathering of silicate and carbonate minerals, is a key component of the global carbon cycle and Earth's long-term climate (Berner, 1999; Gaillardet et al., 1999; Hilton and West, 2020; Torres et al., 2017). The degree to which carbonate alkalinity production involves CO₂ (as carbonic acid, H₂CO₃ = H₂O + CO_2(g,aq)) from atmospheric or soil respiration sources and liberates mineral carbon influences whether dissolved inorganic carbon (DIC = Σ[CO₂, carbonate alkalinity]) in fluvial networks represents a carbon sink or source. Rapid warming at northern latitudes (Serreze and Barry, 2011) is thawing permafrost (Biskaborn et al., 2019), increasing vegetation productivity (Bjorkman et al., 2018), intensifying hydrologic cycles (Rawlins et al., 2010), and strengthening land–fresh water 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 DIC cycling (Lacelle et al., 2019; Wadham et al., 2019). In recent decades, increasing riverine fluxes of carbonate alkalinity and solutes across the circumpolar north reflect enhanced mineral weathering associated with active layer thickening, deepening hydrologic flow paths into mineral soils, and greater soil acidity from increasing vegetation productivity (Drake et al., 2018a; Tank et al., 2016; Toohey et al., 2016). Glaciated permafrost terrain hosting ice-rich deposits of reactive sediments are thought to be distributed across the northern permafrost zone, raising the prospect that terrain consolidation following thaw (thermokarst) and associated carbonate alkalinity production and export may have stronger influence on climate feedbacks in such regions (Zolkos et al., 2018).

Three coupled factors primarily influence the degree to which carbonate alkalinity represents a carbon sink or source: first, the weathering source, which accounts for both the mineral composition of substrate subjected to chemical weathering and the acid(s) responsible for weathering. Silicate weathering by H₂CO₃ generates alkalinity without liberating mineral carbon and thus represents a long-term CO₂ sink. In contrast, H₂CO₃ carbonate weathering is a CO₂ sink only over ∼ 10²–10³ yr as half of the alkalinity produced is geogenic. ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ produced during carbonate weathering in the presence of strong acids, for instance sulfuric acid (H₂SO₄) from sulfide oxidation, is a CO₂ source over longer timescales (∼ 10⁶ yr; Calmels et al., 2007) and can also produce CO₂ over shorter timescales when H₂SO₄ is present in excess (Stumm and Morgan, 1996). The second factor is the rate of mineral weathering and processes that further alter this rate. 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 often fast in glacial environments owing to intense physical weathering (Anderson, 2007). Indeed, hydrochemical signatures of trace 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, H₂CO₃ carbonate weathering in proglacial streams can consume dissolved CO₂ to below atmospheric levels (Sharp et al., 1995; St. Pierre et al., 2019). The third factor is the magnitude of carbonate alkalinity export, which is influenced by its production via weathering of minerals during fluvial transport (e.g., Striegl et al., 2007) and its loss via carbonate equilibrium reactions and CO₂ degassing along the land–fresh water–ocean continuum. From a climate perspective, the magnitude of carbonate alkalinity export is particularly relevant over geological timescales because half of riverine carbonate alkalinity exported to the ocean is returned to the atmosphere as CO₂ via precipitation reactions within the marine carbon cycle (Calmels et al., 2007). Together, these three controls on carbonate alkalinity highlight the nonconservative nature of DIC and its susceptibility to transformation within fluvial networks. Hence, to constrain carbonate alkalinity export in rapidly changing permafrost terrains, nested-watershed sampling designs are critical for capturing DIC transformation along the land–fresh water–ocean continuum and resolving drivers and sources of carbon cycling across scales (Drake et al., 2018b).

Figure 1 Sampling sites on the Peel Plateau (NWT, Canada). Water samples were collected along the main stem of Dempster and Stony creeks (n=12) and major tributaries (n=10) and from the rill runoff at retrogressive thaw slump (RTS) FM2. Numbers within symbols are sampling sites (Tables 1 and A1). RTS impact accumulation represents the number of active RTSs affecting upstream reaches (n=109; see Methods Sect. 2.6). (a) Aerial photograph of Stony Creek where it was first impacted by RTS activity. (b) RTS FM2 runoff transect sampling scheme. RTS FM2 spans ∼ 40 ha, its headwall (c) reaches ∼ 25 m in height, and the debris tongue contains 2 × 10⁶ m³ of sediment (van der Sluijs et al., 2018). (d) Aerial photograph of Stony Creek (lower left) flowing into the Peel River. Satellite image of RTS FM2 in September 2017 (b) obtained from Copernicus Sentinel data (European Space Agency, https://sentinel.esa.int/, last access: 25 March 2020). Base map: Esri ArcGIS Online© OpenStreetMap contributors, GIS User Community.

Glaciated permafrost terrains are poised for rapid geomorphic and associated biogeochemical change as the climate warms and precipitation intensifies (Kokelj et al., 2017b). Despite glacial retreat across much of the circumpolar north, permafrost within these landscapes preserves biogeochemical legacies of glaciation across northern Canada, Alaska, and western Siberia (Kokelj et al., 2017b). In North America, the comminution of carbonate and shale bedrock during expansion of the Laurentide Ice Sheet (LIS) and the climate and vegetative protection of ice- and sediment-rich tills in the wake of its retreat endowed former glacial margins across northwestern Canada with thick inorganic tills held in ice-rich permafrost (Kokelj et al., 2017b). Today, the climate-driven renewal of postglacial landscape change is mobilizing immense stores of minerals into modern biogeochemical cycles via hillslope thermokarst features, the largest of which include retrogressive thaw slumps (RTSs; Kokelj et al., 2017a). On the Peel Plateau (NWT, Canada), RTSs expose carbonate- and sulfide-bearing glaciogenic permafrost sediments that are tens of meters thick. The chemical weathering and fluvial transport of these sediments results in increased ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ immediately downstream of RTSs and greater solute and sediment loads throughout downstream systems (Kokelj et al., 2013; Malone et al., 2013; Zolkos et al., 2018). RTS activity has been suggested – but not previously proved to be – partly responsible for increasing carbonate alkalinity fluxes in the larger Peel River during recent decades (Zolkos et al., 2018). Yet, it remains unknown how hillslope 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 across the circumpolar north. In this study we evaluated trends in major ions, DIC concentration, and dual δ¹³C-DIC–δ¹³C-CO₂ isotopes along transects within three nested watersheds in the Stony Creek watershed on the Peel Plateau. Our nested watershed approach was intended to allow us to determine how RTS effects on carbon cycling integrate across scales from ∼ 1 to 1000 km².

2 Methods

2.1 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 slopes are sparsely vegetated and mainly consist of bedrock colluvium (Duk-Rodkin and Hughes, 1992). Exposed marine shale and sandstone bedrock contain sulfide- and gypsum-bearing lithologies but limited carbonate (Norris, 1985). As Stony Creek flows eastward, the main channel and its tributaries incise ice-rich glacial deposits 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). Growth of RTSs 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 in stream valleys (Fig. 1). Across the Stony Creek watershed, intensifying RTS activity releases large volumes of 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² watershed of Dempster Creek, a tributary of Stony 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 (Kokelj et al., 2013; 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.

2.2 Stream sampling

In late July 2017, we sampled along transects within three nested watersheds (Fig. 1, Table A1) to understand how the effects of RTSs integrate across watershed scales. (i) The RTS FM2 runoff transect included five sampling locations along a 550 m long thaw stream formed by a runoff channel within an active RTS. The RTS FM2 runoff received no observable hydrologic inputs during the sampling period. (ii) A 14 km transect along the main stem of Dempster Creek, which received inputs directly from RTS FM2, was sampled at one location in undisturbed headwaters and at three sites downstream of RTS FM2. Sites downstream were located on the main stem, immediately upstream of three major tributaries. We also sampled from the tributaries near their confluence with Dempster Creek to characterize tributary chemistry. (iii) A 70 km transect along the main stem of Stony Creek, a sixth-order stream, was sampled at eight locations: one in undisturbed headwaters and seven on the RTS-affected reach upstream of major tributaries. We additionally sampled from one tributary of the undisturbed headwaters and from six RTS-affected tributaries near their confluence with the main stem. Stony Creek, a major tributary of the 70 000 km² Peel River watershed (Fig. 1), was sampled to determine how the effects of RTS activity on DIC integrate across broader scales.

At all sampling sites, stream temperature, specific conductance (henceforth, “conductivity”), and pH were measured using a precalibrated YSI Professional Plus water quality meter. At most sites, water samples were collected for ions, DIC, CO₂, CH₄, dissolved organic carbon (DOC), UV–visible absorbance, and total suspended solids (TSS). Along the RTS FM2 runoff transect, we sampled only for DIC and CO₂ concentration and stable isotopes of dissolved CO₂ (δ¹³C-CO₂). Additional parameters were sampled 1 d prior at RTS FM2 runoff site 5, 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 (δ¹³C-DIC) and used dual δ¹³C-DIC–δ¹³C-CO₂ isotopes to characterize abiotic and biotic processes influencing DIC sources and cycling across watershed scales.

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 high-density polyethylene (HDPE) bottle in the same fashion. Water samples were filtered using sample-rinsed 0.45 µm polyethersulfone (PES; ThermoFisher) or cellulose acetate (CA; 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 δ¹³C-DIC. Samples for cations were filtered (CA) into acid-washed bottles and acidified with trace metal-grade HNO₃, 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 nonacid-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 airtight syringes (for CO₂ concentration) or overpressurized in pre-evacuated serum bottles sealed with prebaked (60 ^∘C, 12 h), gas-inert butyl rubber stoppers (for δ¹³C-CO₂, CH₄). At each site, atmospheric samples for CO₂ and CH₄ concentration and δ¹³C-CO₂ were stored in the same fashion. Gas samples were stored in the dark at ∼ 20 ^∘C prior to analysis within 10 h (CO₂) or 2 months (δ¹³C-CO₂, CH₄). Water and air temperature, atmospheric pressure, and the volumetric ratio of sample to atmospheric headspace were recorded for correcting later calculations of CO₂ partial pressure (pCO₂) and δ¹³C-CO₂ (Hamilton and Ostrom, 2007).

Table 1 Mineral weathering equations used to create Piper diagram end-members. H₂CO₃: carbonic acid; H₂SO₄: sulfuric acid. H₂CO₃ includes dissolved CO_2(g).

2.3 Hydrochemical analyses

Upon returning from the field each day, CO₂ was measured using an infrared gas analyzer (PP Systems EGM-4), which was checked monthly for drift using a commercial standard (Scotty Gases). We calculated pCO₂ 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 CO₂ detection (LiCOR 7000) following acidification within a DIC analyzer (Apollo SciTech model AS-C3). Calibration curves were made daily using certified reference material (CRM) from the 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 (CO₂, ${\mathrm{HCO}}_{\mathrm{3}}^{-}$, ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$) were calculated from DIC concentration and pCO₂ or pH using CO2SYS (v.2.3; Pierrot et al., 2006) as well as field temperature and pressure at the time of sampling and the freshwater equilibrium constants for K1 and K2 (Millero, 1979).

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 no. 17025). DOC was measured using a total-organic-carbon analyzer (Shimadzu TOC-V). DOC standard curves were made daily with a 1000 ppm KHP solution (ACCUSPEC), and an in-house caffeine standard (10 mg L⁻¹) 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 (SUVA₂₅₄) to infer DOC aromaticity (Weishaar et al., 2003).

δ¹³C-DIC was determined using an isotope ratio mass spectrometer (Finnigan Mat DeltaPlusXP) interfaced to a total-organic-carbon analyzer (OI Analytical Aurora Model 1030W) at the University of Ottawa Stable Isotope Laboratory. δ¹³C-CO₂ and the concentration of CH₄ were analyzed within 2 months using a Picarro isotope analyzer (G2201-i; < 0.2 ‰ precision, CH₄ operational range = 1.8–1500 ppm) equipped with a Small Sample Introduction Module. Commercial δ¹³C-CO₂ and CH₄ standards were used to check for drift during each run. We used mass balance to correct δ¹³C-CO₂ values for the δ¹³C and mass of atmospheric CO₂ used for equilibration (Hamilton and Ostrom, 2007). To assess δ¹³C-CO₂ 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 CH₄ (pCH₄). SI was calculated using the hydrochemical software Aqion version 6.7.0 (https://www.aqion.de, last access: 15 October 2018), which uses the US Geological Survey software PHREEQC (Parkhurst and Appelo, 2013) as the internal numerical solver. Samples for atmospheric and dissolved CH₄ were collected in the same fashion as δ¹³C-CO₂. pCH₄ was calculated 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 preweighed 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 US Geological Survey Method I-3765.

2.4 Mineral weathering and DIC sources

We used a Piper diagram (Piper, 1944), which reflects the proportional equivalent concentrations of ions in a sample relative to mineral weathering end-members, as one method to constrain the sources of mineral weathering and ${\mathrm{HCO}}_{\mathrm{3}}^{-}$. The products of Eqs. (1)–(7) defined the mineral weathering end-members in the Piper diagram (Table 1). We further constrained mineral weathering and DIC sources using δ¹³C-DIC and pH. End-member δ¹³C-DIC ranges for equilibrium processes (mixing with atmospheric and/or biotic CO₂) and kinetic reactions (mineral weathering) were derived following Lehn et al. (2017) and using published isotopic fractionation factors (Zhang et al., 1995).

To evaluate potential effects on δ¹³C-CO₂ from DIC speciation along the pH continuum (Eq. 8, Table 1; Clark and Fritz, 1997), we compared theoretical and observed δ¹³C-CO₂ values in the Stony Creek main stem. Theoretical δ¹³C-CO₂ values were calculated using mass balance to obtain δ¹³C–${\mathrm{HCO}}_{\mathrm{3}}^{-}$ from measurements of DIC, CO₂, ${\mathrm{HCO}}_{\mathrm{3}}^{-}$, δ¹³C-DIC, and δ¹³C-CO₂. We then used measurements of stream temperature (T) to calculate the equilibrium fractionation between CO₂ and ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ (ε = −9.483 × 10${}^{\mathrm{3}}/T+\mathrm{23.89}\mathrm{‰}$; Mook et al., 1974). Finally, ε was subtracted from δ¹³C–${\mathrm{HCO}}_{\mathrm{3}}^{-}$ to obtain theoretical δ¹³C-CO₂. Similarity between observed and theoretical δ¹³C-CO₂ values was interpreted as δ¹³C-CO₂ variability driven by carbonate equilibrium reactions, whereas dissimilarity was taken to reflect effects from CO₂ degassing (Zhang et al., 1995) and/or biotic CO₂ production (Kendall et al., 2014).

2.5 Geospatial analyses

Stream networks and watershed areas were delineated using the ArcHydro tools in ArcGIS 10.5 from the gridded (30 m) Canadian Digital Elevation Model (CDEM). CDEM data were reconditioned using National Hydro Network stream vectors, which were first modified as needed to align with stream flow paths visible in Copernicus Sentinel-2 multispectral satellite imagery from 2017 (European Space Agency, https://sentinel.esa.int/, last access: 25 March 2020). To statistically assess landscape controls on DIC cycling (Sect. 2.7), we delineated active RTSs and derived terrain roughness and vegetation productivity in the major tributary watersheds of Stony Creek. RTSs were interpreted as active where exposed sediment visibly dominated the feature surface (Cray and Pollard, 2015) in orthorectified Satellite Pour l'Observation de la Terre (SPOT) multispectral imagery that we pan-sharpened to 1.6 m resolution using the ArcGIS Image Analysis tool. The satellite imagery was collected from 9 to 25 September 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 the accumulation of RTS impacts across the fluvial network. We defined RTS impact accumulation as the cumulative number of active RTSs impacting upstream reaches. RTSs were interpreted to impact streams based on contact with the channel or interpreted downslope flow based on slope direction and gradient from the CDEM (Supplement). 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 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 US National Aeronautics and Space Administration EVI product (Didan, 2015), which is derived from the gridded (250 m) moderate-resolution imaging spectroradiometer (MODIS). The MODIS data were collected on 28 July 2017. The ArcGIS Zonal Statistics tool was used to calculate total RTS area, mean terrain roughness, and mean EVI in Stony Creek tributary watersheds.

2.6 Stream flow

Water discharge (Q) in Stony Creek tributaries was estimated from a hydraulic-geometry model (Gordon et al., 2004) that we developed using flow measurements made in Peel Plateau streams during 2015–2017, and width (W) was estimated from on-site measurements or photos from 2017 with a known scale. The model reflected measurements spanning diverse stream morphologies (W = 0.4–6.6 m) and flow conditions (Q = 0.005–0.91 m³ s⁻¹; Fig. A1):

$$\begin{array}{}\text{(1)}& Q=e^{\ln (W/\mathrm{6.258})/\mathrm{0.661}}(p<\mathrm{0.001},R^{\mathrm{2}}=\mathrm{0.89},F_{\mathrm{1},\mathrm{18}}=\mathrm{150}).\end{array}$$

Discharge values from 2015 to 2017 were calculated from measurements of stream flow (RedBack Model RB1, PVD100) and cross-sectional area made at increments equal to 10 % of stream width (Gordon et al., 2004; Lurry and Kolbe, 2000) and were averaged for sites with multiple measurements.

2.7 Statistics

We tested for downstream change in ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ concentration and pCO₂ along the Stony Creek main stem using the nonparametric Mann–Kendall test from the R software (R Core Team, 2018) package zyp (Bronaugh and Werner, 2013) following the trend prewhitening 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 ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ export in Stony Creek tributary watersheds relative to other landscape variables known to influence DIC production, including hydrology, terrain roughness, and vegetation productivity (Berner, 1992; Drake et al., 2018a). To account for potential effects of varying tributary watershed areas on discharge (Q) and constituent concentration, we used tributary ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ yields in the model. Instantaneous discharge (Q; m³ s⁻¹) was estimated from the hydraulic-geometry relationship between Q and stream width (Eq. 1). Discharge and ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ flux (concentration ⋅ Q) were normalized to the respective tributary watershed area and scaled to estimate daily water yield (cm d⁻¹) and ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ yield ($\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$). Daily ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ yields in Stony Creek tributaries were modeled as

$$\begin{array}{}\text{(2)}& {\mathrm{HCO}}_{\mathrm{3}}^{-}\text{yield}={\text{RTS}}_n+{\text{RTS}}_{\text{area}}+\text{Water yield}+\text{TR}+\text{EVI},\end{array}$$

where RTS_n is the number of active RTSs, RTS_area 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 covariates were transformed as needed to meet assumptions of independent and homoscedastic residuals (Zuur, 2009). To understand potential effects from variable rainfall on water yields prior to and during the 2 d 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 Northwest Territories total meteorological station located ∼ 1 km from the RTS FM2 (Fig. A2). Statistics were performed in the R programming environment (v.3.4; R Core Team, 2018), and significance was interpreted at α = 0.05. Summary statistics are reported as mean ± standard error unless noted.

Table 2 Geochemistry of main stem and tributary sites along Dempster and Stony creeks. Retrogressive thaw slump (RTS) FM2 runoff samples collected on 31 July 2017, except where noted. RTS FM2 runoff site 5 was nearest to the confluence with Dempster Creek (Fig. 1). Area: watershed area; SE: standard error.

^* Not RTS-affected. ^†30 July 2017.

3 Results

3.1 pH, ions, and weathering sources across watershed scales

Geochemistry of the main stem and tributary sites are summarized in Table 2. Among sites, pH was generally circumneutral, and conductivity was higher in proximity to RTS activity. pH was highest in the RTS FM2 runoff (7.69 ± 0.05, mean ± SE), intermediate in Dempster Creek (7.07 ± 0.42), and lowest in Stony Creek (6.86 ± 0.21). Along the RTS FM2 runoff transect, pH decreased from 7.72 to 7.51 between sites 1 and 2, and increased thereafter to 7.80 at site 5. pH in the Dempster Creek headwaters (5.82) was lower than in the RTS-affected reach (7.48 ± 0.1). In Stony Creek, pH increased from 5.66 in headwaters to ∼ 7.30 at sites 6–8.

Similar to pH, conductivity was higher in the RTS FM2 runoff (1799 ± 111 µS cm⁻¹) than in Dempster Creek (520 ± 191 µS cm⁻¹) and Stony Creek (320 ± 19 µS cm⁻¹). Conductivity in RTS FM2 runoff increased from 1370 to 1990 µS cm⁻¹. Along Dempster Creek, conductivity increased from 52 µS cm⁻¹ in the undisturbed headwaters to 958 µS cm⁻¹ at the first site downstream of RTS FM2 and decreased downstream thereafter. In Stony Creek, conductivity decreased between the headwaters and the fourth downstream site and was relatively constant at ∼ 285 µS cm⁻¹ along the lower reach of Stony Creek (sites 5–8).

Figure 2 Piper diagrams (modified to show the upper half of the diamond plot) showing stream chemistry of the (a) main stem and retrogressive thaw slump (RTS) FM2 runoff sites and (b) tributary sites. Axes and corresponding text in gray and black reflect the proportions of cations and anions, respectively. Mineral weathering end-members were derived from the proportional concentration (meq L⁻¹) of solutes generated by H₂CO₃ carbonate weathering (CACW; Eq. 1), H₂SO₄ carbonate weathering (SACW; Eq. 3), H₂SO₄ silicate weathering (SASW; Eq. 4), sulfate salt (e.g., gypsum) dissolution (SSD; Eq. 5), and carbonate weathering by H₂SO₄ in excess (SA_exCW; Eq. 7). Site numbers given within symbols (Table A1). ^* Not RTS-affected.

Streams were characterized by Ca²⁺–Mg²⁺–${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$-type waters (Fig. 2) with low concentrations of Cl⁻ relative to ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$, reflecting a predominance of H₂SO₄ carbonate weathering and sulfate salt (e.g., gypsum) dissolution over other mineral weathering sources. A relatively greater proportion of ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ than ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ in the RTS FM2 runoff and along the upper reach of Stony Creek (sites 1–4; Fig. 2a) suggests greater sulfate salt dissolution and/or that carbonate weathering at these sites buffered less H₂SO₄ (Eq. 7) than in Dempster Creek headwaters and its tributaries (Eq. 3). Along the Stony Creek main stem (sites 1–8), increasing ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ (Fig. 2a) reflected inputs from RTS-affected tributaries (sites 2–7) having relatively more ${\mathrm{HCO}}_{\mathrm{3}}^{-}$-type waters (Fig. 2b) from H₂SO₄ and potentially H₂CO₃ carbonate weathering.

Figure 3 (a) ${\mathrm{HCO}}_{\mathrm{3}}^{-}$, (b) pCO₂, (c) dissolved organic carbon (DOC), and (d) SUVA₂₅₄ along the RTS FM2 runoff transect and the main stem of Stony and Dempster creeks (see locations in Fig. 1). Note different x-axis scales for the FM2 runoff transect (0–0.55 km, upper x axis), Dempster Creek (0–14 km, below RTS FM2 x axis), and Stony Creek (0–70 km, lower x axis). For the RTS FM2 runoff, DOC and SUVA₂₅₄ were sampled only at 0.55 km. Regression lines in (a) and (c) are from a Mann–Kendall test (details in Sect. 2.7). Bars on x axes indicate where RTS FM2 runoff enters the Dempster Creek transect (3.3 km) and where Dempster Creek enters Stony Creek (56 km). Site numbers are given within symbols (Table A1). ^* Not RTS-affected.

3.2 ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ concentration and pCO₂

Carbonate alkalinity (${\mathrm{HCO}}_{\mathrm{3}}^{-}$ + ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$) was primarily ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ (> 99 %) at all sites. ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ was highest in the RTS FM2 runoff (1429 ± 23 µM), intermediate in Dempster Creek (864 ± 261 µM), and lowest in Stony Creek (255 ± 59 µM). Along the RTS FM2 runoff transect, ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ decreased from 1510 to 1386 µM. In Dempster Creek and Stony Creek, ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ concentrations were relatively low in undisturbed headwaters (115 and 33 µM, respectively) and 2 to 10 times higher at the first RTS-affected site (1321 and 69 µM, respectively). ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ decreased along the entire RTS-affected reach of Dempster Creek (from 1321 to 946 µM) in conjunction with inputs from dozens of tributary watersheds without active RTSs. In contrast, ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ increased significantly along Stony Creek (p < 0.01, Mann–Kendall test; Fig. 3a) in conjunction with inputs from RTS-affected tributaries.

CO₂ was oversaturated at all sites (Fig. 3b) and a minor component of DIC (< 10 %) at most sites except the undisturbed headwaters of Dempster Creek (site 1) and upper Stony Creek (sites 1–3). pCO₂ was highest in the Dempster Creek headwaters (2467 µatm), relatively high in the RTS FM2 runoff (1023 ± 137 µatm), and consistently near atmospheric levels along Stony Creek (479 ± 12 µatm). Along the RTS FM2 runoff transect, pCO₂ increased from 1046 to 1534 µatm within the first 220 m and then decreased from 1534 to 742 µatm over the final 330 m. Along Dempster Creek, pCO₂ decreased from 2467 µatm in the headwaters to 686 µatm at the first RTS-affected site and further decreased to 600 µatm by the end of Dempster Creek. pCO₂ in Dempster and Stony Creek tributaries were generally similar to the main stem sites.

3.3 DOC concentration and SUVA₂₅₄

DOC concentrations were highest in Dempster Creek (933 ± 83 µM), intermediate in the RTS FM2 runoff (758 ± 152 µM), and lowest in Stony Creek (303 ± 54 µM). Along the Dempster Creek transect, DOC decreased between the undisturbed headwaters (960 µM) and the first RTS-affected site (790 µM) and increased thereafter along the transect (to 1156 µM; Fig. 3c). Along Stony Creek, DOC increased significantly (from 102 to 551 µM; p < 0.001, Mann–Kendall test).

SUVA₂₅₄ values were lowest in the RTS FM2 runoff (1.85 ± 0.4 $\mathrm{L}{\mathrm{mgC}}^{-\mathrm{1}}{\mathrm{m}}^{-\mathrm{1}}$), highest in Dempster Creek (3.10 ± 0.2 $\mathrm{L}{\mathrm{mgC}}^{-\mathrm{1}}{\mathrm{m}}^{-\mathrm{1}}$), and intermediate in Stony Creek (2.51 ± 0.3 $\mathrm{L}{\mathrm{mgC}}^{-\mathrm{1}}{\mathrm{m}}^{-\mathrm{1}}$). SUVA₂₅₄ values along the Dempster Creek transect followed a similar pattern to DOC, and along Stony Creek SUVA₂₅₄ values doubled (Fig. 3d). DOC and SUVA₂₅₄ increased in consecutive downstream tributaries of Stony Creek but not Dempster Creek.

Figure 4 The pH and composition of dissolved inorganic carbon stable isotopes (δ¹³C-DIC) in streams. The upper and lower reference lines depict theoretical end-members for equilibrium reactions (mixing with atmospheric and biotic CO₂, respectively). Gray boxes span theoretical end-member values for kinetically controlled mineral weathering reactions (SACW: H₂SO₄ carbonate weathering; CACW: H₂CO₃ carbonate weathering; CASW: H₂CO₃ silicate weathering; see Sect. 2.4 for derivation of end-members). The vertical line corresponds to the pH at which ≥ 90 % of DIC is ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ for the mean observed stream water temperature (11.7 ^∘C). At pH < 7.4, δ¹³C-DIC values primarily reflect equilibrium (rather than kinetic) controls on DIC cycling. Arrows reflect increasing downstream distance from the headwaters in Stony Creek and from the first retrogressive thaw slump (RTS)-affected site in Dempster Creek. Site numbers given within symbols (Table A1). ^* Not RTS-affected.

3.4 Stable isotopic composition of carbon in DIC and CO₂

δ¹³C-DIC values were highest in the RTS FM2 runoff (−1.0 ‰) and lower, on average, along the main stem of Dempster Creek (−7.5 ‰ ± 2.5 ‰) and Stony Creek (−8.4 ‰ ± 0.5 ‰). In the undisturbed headwaters of Dempster Creek and Stony Creek, relatively negative δ¹³C-DIC values (−11.6 ‰ to −15.6 ‰) reflected DIC sourced from a combination of atmospheric and biogenic (soil) CO₂ (Fig. 4). In the RTS FM2 runoff, relatively ¹³C-enriched δ¹³C-DIC (−1.0 ‰) aligned with H₂SO₄ carbonate weathering. δ¹³C-DIC decreased from −4.2 ‰ at the first site downstream of the RTS FM2 runoff to −5.7 ‰ at the end of Dempster Creek. Along Stony Creek, δ¹³C-DIC increased from the undisturbed headwaters (−11.6 ‰) to the most downstream site (−7.8 ‰). δ¹³C-DIC signals of H₂SO₄ carbonate weathering diminished slightly downstream along the Dempster Creek transect and intensified along Stony Creek (Fig. 4).

Figure 5 The composition of DIC and CO₂ stable isotopes at varying DIC concentrations along the Dempster Creek and Stony Creek main stems and in the rill runoff of retrogressive thaw slump (RTS) FM2. Arrows reflect increasing downstream distance from headwaters in Stony Creek, from the first RTS-affected site in Dempster Creek, and from the start of the FM2 runoff transect. Site numbers given within symbols (Table A1). ^* Site in headwaters and not affected by RTSs.

Similar to δ¹³C-DIC, δ¹³C-CO₂ values were higher in the RTS FM2 runoff (−11.0 ‰ ± 0.4 ‰) than along the main stem of Dempster Creek (−17.9 ‰ ± 1.4 ‰) and Stony Creek (−16.7 ‰ ± 0.6 ‰; Fig. 5). δ¹³C-CO₂ values were relatively low in the undisturbed headwaters of Dempster Creek (−21.6 ‰) and intermediate in the headwaters of Stony Creek (−13.8 ‰; Fig. 5). Along the RTS FM2 runoff transect, δ¹³C-CO₂ values increased from sites 1 to 4 (−12.1 ‰ to −10.0 ‰) and decreased at site 5 (−11.2 ‰). Along the RTS-affected reach of Dempster Creek, δ¹³C-CO₂ values decreased from −16.0 ‰ to −18.5 ‰ in conjunction with inputs from non-RTS-affected tributaries having relatively low δ¹³C-CO₂ (−18.7 ‰ ± 1.4 ‰) that was more similar to values from soil-respired CO₂. Along Stony Creek, δ¹³C-CO₂ values decreased from −13.8 ‰ to −18.1 ‰, showing a trend opposite that of δ¹³C-DIC (Fig. 5). Among sites, atmospheric δ¹³C-CO₂ values were relatively consistent (−9.5 ‰ ± 0.4 ‰, mean ± SD).

Variance in δ¹³C of CO₂ and DIC could be influenced by biotic production, CO₂ conversion to ${\mathrm{HCO}}_{\mathrm{3}}^{-}$, and/or mixing with atmospheric CO₂. To evaluate the relative influence of these processes, we compared measured δ¹³C-CO₂ for Stony Creek with theoretical values reflecting DIC controlled by speciation along the pH continuum (Sect. 2.4). In the undisturbed headwaters, δ¹³C-CO₂ indicated stronger influence from atmospheric CO₂ (Fig. 6). Along the upper, RTS-affected reach of Stony Creek (sites 2–5, from ∼ 5 to 35 km), the good agreement between measured and theoretical δ¹³C-CO₂ values reflected equilibrium fractionation (ε = 9.7 ‰ at 9 ^∘C; Mook et al., 1974) between CO₂ and ${\mathrm{HCO}}_{\mathrm{3}}^{-}$, indicating greater influence from DIC speciation (Fig. 6). Along the lower RTS-affected reach of the transect (sites 6–8), δ¹³C-CO₂ values more strongly reflected biotic CO₂ production, with potential effects from degassing and/or CO₂ conversion to ${\mathrm{HCO}}_{\mathrm{3}}^{-}$. These trends in δ¹³C-CO₂ values along Stony Creek show a downstream change in the processes influencing DIC source, which may be related to inputs of weathering solutes and organic matter from RTS-affected tributaries.

Table 3 Characteristics of Stony Creek tributary watersheds (upper panel) and results from the multiple linear regression model (lower panel). The Stony Creek watershed contained 109 retrogressive thaw slumps (RTSs), 92 of which were in the major tributaries of the Stony Creek main stem. TR: terrain roughness. EVI: enhanced vegetation index. Values were used in the multiple linear regression model to determine the drivers of ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ yields in Stony Creek tributary watersheds. Model results are shown in the lower panel. Covariates eliminated during model selection (RTS_n, TR, EVI) are not reflected in the lower panel or final model: $\ln {\mathrm{HCO}}_{\mathrm{3}}^{-}\text{yield}$ = 1.04ln water  yield + 0.35ln RTS_area + 8.76.

3.5 Stony Creek tributary carbonate alkalinity yields and watershed characteristics

Carbonate alkalinity yields in RTS-affected tributaries of Stony Creek (1558 ± 1135 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, mean ± SD) were 3 orders of magnitude higher than in the non-RTS-affected headwaters (1.8 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$; Table 3). Consecutive downstream tributary watersheds exhibited no clear trends in the number of RTSs, the area disturbed by RTSs, terrain roughness, or EVI. In the Stony Creek headwater tributary, which had no active RTSs, terrain roughness (16.2 m) and vegetation productivity (EVI = 0.28) were higher than in the other six tributary watersheds (4.3 ± 1.3 m and 0.46 ± 0.01, mean ± SD, respectively). In the other tributary watersheds, the number of active RTSs reached 50 (15 ± 17, mean ± SD), and RTS disturbance area reached 3.5 % (0.91 % ± 1.29 %, mean ± SD; Table 3).

To elucidate landscape controls on carbonate alkalinity export in Stony Creek tributary watersheds, we paired geospatial data for active RTSs, terrain roughness, and vegetation productivity with estimates of carbonate alkalinity and water yields in a multiple linear regression model (Sect. 2.7). Water yield and the area of RTS disturbance were retained during automated covariate selection for the final model ($F_{\mathrm{2},\mathrm{4}}=\mathrm{63}$, p < 0.001, R²=0.95). In addition to the expected relationship between water yield and carbonate alkalinity yield, RTS disturbance area was a clear, significant predictor of carbonate alkalinity yield and formed a stronger relationship with alkalinity than did water yield (Table 3).

Figure 6 Observed and expected δ¹³C-CO₂ values along the Stony Creek main stem. Theoretical δ¹³C-CO₂ values were calculated as detailed in Sect. 2.6 and reflect changes in CO₂ due to DIC speciation (i.e., H₂CO₃ ⇌ H⁺ + ${\mathrm{HCO}}_{\mathrm{3}}^{-}$; Eq. 8). Deviation from theoretical δ¹³C-CO₂ values by observed values thus indicates isotopic effects from degassing and/or microbial oxidation of organic matter (OM), as indicated by the arrows. Site numbers given within symbols (Table A1). ^* Site was not affected by retrogressive thaw slumps.

4 Discussion

4.1 Rapid carbon cycling in fluvial network headwaters

Within undisturbed headwaters and RTS runoff on the Peel Plateau, rapid carbon cycling enhanced fluvial CO₂ efflux to the atmosphere. In undisturbed headwaters, δ¹³C-CO₂ values indicate inputs of primarily biogenic CO₂ from soil respiration into Dempster Creek. In the Stony Creek headwaters, intermediate δ¹³C-CO₂ and pCO₂ saturation suggested influence from exchange with atmospheric CO₂ and from some biogenic CO₂. In the undisturbed Dempster Creek headwaters, a 70 % decrease in pCO₂ within several kilometers downstream likely reflected degassing and diminishing inputs of respired CO₂ from soils to streams relative to headwaters (Hutchins et al., 2019). These trends resemble headwater streams elsewhere in that hydrologic inputs of respired CO₂ from riparian soils can drive CO₂ supersaturation in fluvial network headwaters (Campeau et al., 2018; Crawford et al., 2013), which is rapidly effluxed to the atmosphere over short distances downstream (Hotchkiss et al., 2015). In contrast, trends in hydrochemistry and stable isotopes within RTS FM2 runoff demonstrate that drivers of carbon cycling within RTSs are starkly different from those in undisturbed headwaters on the Peel Plateau.

Along the RTS FM2 runoff transect, the increase in conductivity corroborates experimental evidence (Zolkos and Tank, 2020) that permafrost sediments on the Peel Plateau can rapidly weather during fluvial transport within runoff. In the upper reach of the runoff transect, near RTS FM2, the decrease in ${\mathrm{HCO}}_{\mathrm{3}}^{-}$, increase in CO₂, and relatively enriched δ¹³C-CO₂ (Fig. 5) indicate rapid production of geogenic CO₂ via H₂SO₄ carbonate weathering (Eq. 7) and carbonate equilibrium reactions (Eq. 8). In Yedoma terrains in Siberia and Alaska, where mineral soils are relatively more organic-rich, thermokarst is associated with production of biogenic CO₂ (Drake et al., 2018b). While respiration likely produced some CO₂ in RTS FM2 runoff (Littlefair et al., 2017), observed δ¹³C-CO₂ (−11 ‰) more strongly reflected H₂SO₄ weathering of regional carbonate bedrock (−0.7 ‰ to −5.6 ‰; Hitchon and Krouse, 1972) when accounting for isotopic fractionation of ∼ 8 ‰ between carbonate and CO₂ at the temperature of FM2 runoff (18 ^∘C; Clark and Fritz, 1997). Along the lower reach of the FM2 runoff transect, the increase in δ¹³C-CO₂ aligned with the preferential loss of ¹²C in the CO₂ phase via DIC fractionation and degassing (Doctor et al., 2008; Drake et al., 2018b; Kendall et al., 2014). ¹³C enrichment of the CO₂ pool by methanogenesis (Campeau et al., 2018), photosynthesis (Descolas-Gros and Fontungne, 1990), and/or calcite precipitation (Turner, 1982) was unlikely as CH₄ in FM2 runoff was relatively low (pCH₄ = 3.6 ± 1.9 µatm, mean ± SD, n = 6), the high turbidity of FM2 runoff likely inhibited photosynthesis (Levenstein et al., 2018), and calcite was below saturation (SI = −0.79). These trends demonstrate that weathering of sediments during fluvial transport within RTS runoff can result in rapid CO₂ production and efflux to the atmosphere, in agreement with recent estimates of high rates of CO₂ efflux within RTS runoff (Zolkos et al., 2019).

High rates of weathering within RTS FM2 runoff aligns with observations of substantial solute production via the exposure and weathering of carbonate flour in glacial foreground environments (Anderson, 2007; Sharp et al., 1995; St. Pierre et al., 2019). Because minerals exposed by deeper RTSs are generally reactive, and sediment concentrations increased by 3 orders of magnitude between the undisturbed Dempster Creek headwaters and the first RTS-affected site, we reasoned that H₂CO₃ weathering of these sediments during fluvial transport would measurably influence pCO₂ along Dempster Creek (Eq. 1; St. Pierre et al., 2019; Striegl et al., 2007). Although pCO₂ decreased along the RTS-affected reach of the Dempster Creek transect (sites 2–4; Fig. 2b), coincident decreases in conductivity, ${\mathrm{HCO}}_{\mathrm{3}}^{-}$, and pH (Table 2, Figs. 2a and 4) suggest that degassing and dilution associated with inputs from non-RTS-affected tributaries had stronger effects on pCO₂ than did H₂CO₃ carbonate weathering, even at the relatively short scale of this 14 km transect. From a carbon cycling perspective, biogeochemically reactive mineral substrate appears to be rapidly transformed in headwaters on the Peel Plateau; geogenic CO₂ production is relegated to within RTSs; and more stable weathering products, including alkalinity, are exported downstream.

4.2 RTS activity in headwaters amplifies carbonate alkalinity production and accumulation across scales

Similar to CO₂, alkalinity production on the Peel Plateau was strongly coupled to primarily H₂SO₄ carbonate weathering mediated by RTS activity. This was reflected by a modest decrease in ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ along Dempster Creek in tandem with decreasing RTS disturbance area (from 3.2 % to 1.2 %) and some dilution by inputs from non-RTS-affected tributaries. Multiple linear regression results further indicated that RTS activity was a primary terrain control on carbonate alkalinity yields. In the Stony Creek headwaters, low carbonate alkalinity yield relative to water yield suggested that ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ export was limited by carbonate availability rather than by water. In RTS-affected tributaries, higher carbonate alkalinity yields relative to water yields aligned with the model results indicating that RTS activity increases carbonate weathering and alkalinity export beyond what would otherwise be expected on the Peel Plateau. ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ yields in RTS-affected tributaries were comparable to summertime ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ yields in watersheds with carbonate rock weathering by glacial activity (∼ 3000 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$; Lafrenière and Sharp, 2004; Striegl et al., 2007), emphasizing that unmodified sulfide- and carbonate-bearing sediments in regional permafrost are highly reactive (Zolkos and Tank, 2020) and primary sources of DIC within intermediate-sized (1000 km²) fluvial networks. This aligns with stable sulfur isotopes in RTS runoff and near the Stony Creek outflow that strongly reflected sulfide oxidation (Zolkos et al., 2018). Unlike CO₂, the increase in ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ by orders of magnitude along Stony Creek in association with inputs from RTS-affected tributaries shows that more chemically stable (i.e., nongaseous) weathering products accumulated across scales. This aligns with previous findings that solutes and sediments from RTSs propagate through fluvial networks (Kokelj et al., 2013; Malone et al., 2013) and suggests that future intensification of RTS activity (Segal et al., 2016) will increase ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ export to downstream environments.

Figure 7 (a) Conceptual model of retrogressive thaw slump (RTS) activity and mineral weathering effects on carbon cycling in glaciated thermokarst terrains like the Peel Plateau. Bio: biogenic; atm: atmospheric; geo: geogenic. (b) RTS effects on CO₂ and DIC (Σ[CO₂, carbonate alkalinity]) observed in this study (dark solid line), projected across broader scales in the modern-day (dark dotted line), and under hypothetical future scenarios of increasing RTS activity (medium- and light-gray dashed lines). Shaded regions along x axis depict relative RTS area approximated for modern-day (black) and for hypothetical future increases in RTS area (medium- and light-gray).

4.3 Integration of RTS effects on carbon cycling across watershed scales

These findings enable us to develop a conceptual model of catchment chemical characteristics and how the effects of RTS activity on carbon cycling integrate across watershed scales on the Peel Plateau (Fig. 7). This model may be generalized to permafrost terrains elsewhere for testing hypotheses related to thermokarst effects on carbon cycling across the land–fresh water–ocean continuum (Tank et al., 2020).

In undisturbed headwaters on the Peel Plateau, DIC was primarily CO₂, and sources of CO₂ varied from relatively more atmospheric in the sparsely vegetated and mountainous Stony Creek headwaters (Fig. 7ai) to more biogenic in the tundra–taiga headwaters of Dempster Creek (Fig. 7aii). Downstream, CO₂ loss and mixing of streams resulted in undisturbed headwaters having relatively modest DIC comprised of a relatively large proportion of CO₂ sourced from mixing with the atmosphere and likely some inputs from soil respiration (Fig. 7aiii). Underlying the trends in CO₂ concentration, measurements of δ¹³C-CO₂ revealed shifting sources of CO₂ across scales (discussed below).

Thaw and exposure of reactive tills (Lacelle et al., 2019; Zolkos and Tank, 2020) by RTS activity in Peel Plateau headwaters (see also Kokelj et al., 2013; Malone et al., 2013) promotes mineral weathering, rapidly generating CO₂, and substantial alkalinity. Alkalinity, along with large amounts of sediment (van der Sluijs et al., 2018) and organic matter (Shakil et al., 2020), are exported from RTSs into fluvial networks (Fig. 7aiv). Similar to other locations, DOC in RTS runoff on the Peel Plateau is known to be relatively biolabile (Littlefair et al., 2017), suggesting inputs from RTS FM2 to larger streams (e.g., Dempster Creek) could stimulate biotic CO₂ production.

CO₂ degassing is most pronounced within RTSs and in undisturbed headwaters that are strongly coupled with soil respiration, and active mineral weathering is less pronounced in midorder streams (e.g., Dempster Creek). Hence, midorder streams, which also mix with inputs from undisturbed tributaries, export ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ downstream at a magnitude coupled to the area of RTS disturbance (Fig. 7av). Further, immediately downstream of the RTS FM2 inflow to Dempster Creek, the decrease in CO₂ and shift in δ¹³C-CO₂ away from a biotic source suggest that CO₂ degassing to the atmosphere was more prominent than respiration of permafrost DOC (Doctor et al., 2008; Drake et al., 2018b; Kendall et al., 2014). Thus, immediately downstream of RTSs, microbial respiration of permafrost DOC does not appear to generate substantial CO₂. This may be due to lower rates of DOC mineralization than degassing and/or the protection of DOC from microbial oxidation via adsorption to RTS sediments (Gentsch et al., 2015). The latter aligns with the observed decrease in DOC concentration and increase in TSS (to 11 800 mg L⁻¹) between Dempster Creek sites 1 and 2 (Table 2, Fig. 3c; see also Littlefair et al., 2017). However, these effects may diminish farther downstream in midorder streams. Along the lower reach of Dempster Creek (sites 3–4), the decrease in δ¹³C-CO₂, increase in DOC, and SUVA₂₅₄ resembling terrestrial-origin DOC from tributary streams suggest that undisturbed tributary streams may deliver biogenic CO₂ and/or stimulate organic matter respiration in RTS-affected streams. Thus, effects of RTS sediments on CO₂ are attenuated downstream as DOC inputs increase.

Up to and likely beyond scales of ∼ 10³ km² (e.g., Stony Creek), the largest scale of this study, ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ concentrations are likely to increase significantly downstream, reflecting the export of relatively stable weathering products (see also Kokelj et al., 2013; Malone et al., 2013; Zolkos et al., 2018) and accumulation of carbonate alkalinity (Figs. 7avi, b). These effects were primarily driven by inputs of ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ from RTS-affected tributaries, which also increased DOC significantly along Stony Creek. Potentially owing to organic matter limitation, CO₂ in the undisturbed headwaters of Stony Creek appeared to be driven by relatively faster carbonate equilibrium reactions (Eq. 8; Stumm and Morgan, 1996). In contrast, along the lower RTS-affected reach of Stony Creek as ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ and DOC increased and pH stabilized, δ¹³C-CO₂ measurements suggest that respiration of organic matter from RTS-affected tributaries contributed to CO₂ oversaturation (Fig. 6). In higher-order streams within RTS-affected fluvial networks, biotic CO₂ production may increase together with ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ concentrations. This trend was not evident in δ¹³C-DIC, which primarily reflected inputs of geogenic DIC from RTS-affected tributaries. Thus, sources of CO₂ may shift across scales in RTS-affected fluvial networks, and measurements of δ¹³C-CO₂ highlight a decoupling between the drivers of CO₂ and ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ at larger scales (Horgby et al., 2019; Hutchins et al., 2020). A stronger signal of biogenic CO₂ production in larger streams than in permafrost thaw streams within thermokarst, as we observed, is opposite to common trends in Yedoma terrains (Drake et al., 2018b) and may partly reflect limitation of organic substrate in Stony Creek headwaters that is relieved by RTS inputs farther downstream (Shakil et al., 2020). Underlying these trends, RTS disturbance area increased along the Stony Creek transect, from 0 % in the undisturbed headwaters to 0.36 % in the tributary watersheds (sites 4–8). Despite RTS activity occupying a small proportion of the landscape, carbonate alkalinity propagated through fluvial networks. These findings directly link intensifying RTS activity on the Peel Plateau (Segal et al., 2016) with signals of increasing weathering and carbonate alkalinity export in the broader Peel and Mackenzie River watersheds (Tank et al., 2016; Zolkos et al., 2018).

4.4 Implications for carbon cycling in northern permafrost regions

Permafrost terrains susceptible to hillslope thermokarst like RTSs occur within and outside of former glacial limits across the circumpolar north (Olefeldt et al., 2016; Zolkos et al., 2018), and variability in geology, glacial activity, climate, and ecosystem history cause permafrost mineral composition to vary between regions. The degree to which carbonate weathering is coupled with sulfide oxidation will determine if mineral weathering is a CO₂ sink (Eq. 1) or source (Eqs. 3 and 7) over the coming millennia (Zolkos et al., 2018). Where thermokarst releases inorganic substrate with limited prior modification – as in larger RTSs on the Peel Plateau (Lacelle et al., 2019; Zolkos and Tank, 2020) – carbon cycling can be expected to be rapid and driven by inorganic processes and strengthen abiotic components of the permafrost carbon–climate feedback (Schuur et al., 2015). Current dynamic-numerical biogeochemical models for the Mackenzie River basin suggest that the ubiquity of sulfide minerals reduces weathering consumption of atmospheric CO₂ by half (Beaulieu et al., 2012). These models do not account for enhanced H₂SO₄ carbonate weathering associated with RTS activity, which our results show is significantly and positively correlated with alkalinity production and export across watershed scales. Further, climate feedbacks associated with RTS activity appear to be scale-dependent. RTSs rapidly generate CO₂, but its outgassing occurs mostly within runoff and comprises a small proportion of watershed-scale fluvial CO₂ efflux (Zolkos et al., 2019). Carbonate alkalinity generated within RTSs represents a much larger positive feedback to climate change, albeit over geological timescales, via carbonate precipitation reactions within the marine carbon cycle (Calmels et al., 2007). Future intensification of RTS activity (Segal et al., 2016) can thus be expected to increase geogenic CO₂ production within headwaters (see also Zolkos et al., 2019) and carbonate alkalinity export across scales (Fig. 7b; Tank et al., 2016). Cross-scale watershed investigations will help to understand these effects across terrains with varying lithologies and permafrost composition and the implications of hillslope thermokarst for climate feedbacks.

5 Conclusions

Climate-driven renewal of geomorphic activity across permafrost-preserved glaciogenic terrain and associated carbonate weathering in the western Canadian Arctic is amplifying aquatic carbon export across scales, despite RTSs disturbing only a fractional proportion of the landscape. Primary consequences include geogenic CO₂ production that is rapid and localized to RTSs and augments soil-respired CO₂ efflux from undisturbed headwater streams. Significant carbonate alkalinity production and export from RTSs project through fluvial networks and likely to Arctic coastal marine environments, forecasting stronger land–fresh water–ocean linkages (Tank et al., 2016) as RTS activity intensifies in glacial margin landscapes across northwestern Canada (Kokelj et al., 2017a). Legacy effects of RTSs on carbon cycling can be expected to persist for millennia, indicating a need for the integration of dynamic-numerical biogeochemical models (Beaulieu et al., 2012) into predictions of weathering–carbon–climate feedbacks (Zolkos et al., 2018) among northern thermokarst terrains (Turetsky et al., 2020).

Appendix A

Figure A1 Estimated vs. measured discharge (Q; p < 0.001, R² = 0.89, F_1,18 = 150) for 20 streams in the Stony Creek watershed. The gray band represents the 95 % confidence interval shown around the regression. Estimates were made using measurements of stream width, Q, and a hydraulic-geometry model (Gordon et al., 2004; see Sect. 2.6). The model (Eq. 1) was used to estimate Q in the Stony Creek tributaries.

Figure A2 Total rainfall in 24 h increments preceding the sampling of each Stony Creek tributary. Rainfall data were obtained from a Government of Northwest Territories weather station on the Peel Plateau located near the RTS FM2. Locations of tributary sampling sites and the weather station are shown in Fig. 1. ^* Indicates no rainfall in the 24 h window.

Table A1 Sampling site characteristics. Retrogressive thaw slump (RTS) FM2 runoff was a tributary to Dempster Creek (confluence upstream of site 2), and Dempster Creek was a tributary to Stony Creek (confluence upstream of site 8). Coordinates reported in decimal degrees (DDs).

^* Not RTS-affected.