Patterns and persistence of hydrologic carbon and nutrient export from collapsing upland permafrost

.

upland thermokarst may be a dominant linkage transferring carbon and nutrients from terrestrial to aquatic ecosystems as the Arctic warms.

Introduction
Arctic tundra and boreal forest have accumulated a vast pool of organic carbon, twice as large as the atmospheric carbon pool and three times as large as the carbon contained by all living things (Schuur et al., 2015;Hugelius et al., 2014;Tarnocai et al., 2009).Climate change is simultaneously causing widespread permafrost degradation (Slater and Lawrence, 2013) and altering high-latitude hydrology (Peterson et al., 2006;Rawlins et al., 2010), exposing carbon and other elements previously protected in permafrost to transport and processing in Arctic rivers, lakes, and estuaries.Fluxes of dissolved organic carbon (DOC), nutrients, and other ions are changing across the permafrost region and the rate of change is projected to accelerate (Frey and McClelland, 2009;Jones et al., 2005;Laudon et al., 2012;McClelland et al., 2007;McClelland et al., 2014;O'Donnell et al., 2012;Petrone et al., 2006;Rawlins et al., 2010;Striegl et al., 2005;Tank et al., 2012).The interaction between changing hydrology and degrading permafrost is one of the key uncertainties in predicting the response of aquatic ecosystems to high-latitude climate change (Abbott et al., 2015;Koch et al., 2013b;McClelland et al., 2008;Rawlins et al., 2010;Vonk and Gustafsson, 2013).
Permafrost degradation follows two basic trajectories.In permafrost with little ground ice, the soil profile can thaw from the top down without disturbing the surface, gradually exposing organic matter and solutes to hydrologic export as the seasonally thawed active layer deepens (Koch et al., 2013a;Petrone et al., 2006;Striegl et al., 2005).Alternatively, in permafrost where ground ice volume exceeds soil pore space, thaw may cause surface subsidence or collapse, termed thermokarst (Fig. 1; Kokelj and Jorgenson, 2013).When thermokarst occurs on hillslopes it can abruptly mobilize sediment, organic matter, and solutes from meters below the surface, impacting kilometers of stream reach or entire lakes (Bowden et al., 2008;Kokelj et al., 2005;Kokelj et al., 2013;Vonk et al., 2012).
The term thermokarst includes a suite of thermo-erosional features with different morphologies determined primarily by ice content, substrate type, landscape position, and slope (Osterkamp et al., 2009).In upland landscapes, the three most common thermokarst morphologies are retrogressive thaw slumps, active-layer detachment slides, and thermoerosion gullies (Fig. 1; Jorgenson and Osterkamp, 2005;Kokelj and Jorgenson, 2013;Krieger, 2012).Retrogressive thaw slumps (hereafter slumps) often form on lakeshores, have a retreating headwall, and can be caused by a variety of ground ice types including glacial ice, ice wedges, and cave ice.Active-layer detachment slides (hereafter slides) form when the seasonally thawed surface layer of vegetation and soil slips downhill over an ice-rich transition zone (Lewkowicz et al., 2007).Thermo-erosion gullies (hereafter gullies) form due to melting of ice wedges, growing with a generally linear or dendritic pattern, and are often associated with water tracks or headwater streams.These three morphologies currently impact ca.1.5 % of the landscape in the western foothills of the Brooks Range (Krieger, 2012) and could affect 20-50 % of uplands in the continuous permafrost region by the end of the century based on projected thaw and estimates of ground ice distribution (Slater and Lawrence, 2013;Zhang et al., 2000), though circumarctic prevalence and development of upland thermokarst are poorly constrained (Jorgenson et al., 2006;Lantz and Kokelj, 2008;Yoshikawa et al., 2002).
These changes in sediment delivery, light penetration, and nutrients can alter aquatic food webs in receiving ecosystems (Mesquita et al., 2010;Thienpont et al., 2013;Thompson et al., 2012).
Despite a growing understanding of the potential effects of upland thermokarst on aquatic biogeochemical cycles, there is conflicting evidence on the temporal duration of these effects and their overall importance to ecological functioning, precluding conceptualization of patterns of thermokarst impacts and their incorporation into coupled climate models.If thermokarst disturbance is hydrologically connected to aquatic ecosystems, substantial loading of sediment, carbon, and nutrients can occur (Bowden et al., 2008;Kokelj et al., 2005;Kokelj et al., 2013;Shirokova et al., 2013;Thienpont et al., 2013;Vonk et al., 2013), though not all features connected to surface waters result in enhanced carbon and nutrient export (Thompson et al., 2012).Conversely, if thermokarst is hydrologically isolated from surface waters, such as when failures occur high on hillslopes, even dramatic disturbance can have little or no impact on aquatic chemistry and elemental budgets (Lafreniere and Lamoureux, 2013;Lewis et al., 2012).The duration of carbon and nutrient release, and the persistence of biogeochemical disturbance in affected ecosystems after feature stabilization is largely unknown, with altered surface water chemistry lasting for decades in some cases of nutrient loading or surface disturbance (Kokelj et al., 2005;Thienpont et al., 2013), or fading after less than a year in others (Lafreniere and Lamoureux, 2013).
To address these knowledge gaps, we sampled surface outflow from thermokarst features in various stages of development across a broad portion of the North Slope of Alaska.
We focused on two questions.First, how does thermokarst formation alter hydrologic release of carbon and nutrients, and second, can the type and duration of hydrologic release be predicted based on feature morphology or landscape characteristics?We hypothesized that upland thermokarst would initially stimulate nutrient release due to disruption of soil aggregates, accelerated organic matter mineralization in impacted soils, decreased plant uptake, and direct release from melting ground-ice.However, following nutrient retention theory (Vitousek and Reiners, 1975), we hypothesized that this pulse of nutrients would be followed by a period of elemental retention due to enhanced nutrient uptake by recovering vegetation and diminished pools of organic matter and nutrients following disturbance.We hypothesized that DOC export would depend on the balance between DOC production from soil disruption and DOC removal via adsorption by exposed mineral soil as well as enhanced processing of DOC within features due to abundant nutrients and biodegradable DOC from permafrost.In regards to feature morphology, we hypothesized that fundamental differences in formation and functioning of slides, gullies, and slumps, such as the amount of organic and mineral soil displaced, type of ground ice, location on the landscape, and duration of disturbance would result in systematic differences in carbon and nutrient release.We predicted that slumps would have the largest and longest impact, slides would have a large but short-lived impact, and gullies would have a muted impact of intermediate duration.

Study sites
We tested our hypotheses about thermokarst carbon and nutrient export with observations from 83 slides, gullies, and slumps on the North Slope of Alaska (Fig. 2).Features were identified by aerial surveys, satellite imagery, and previous studies (Abbott et al., 2014;Bowden et al., 2008;Gooseff et al., 2009)  The Toolik Field Station is located 254 km north of the Arctic Circle and 180 km south of the Arctic Ocean.The mean annual temperature is -10°C with mean monthly temperatures ranging from -25°C in January to 11.5°C in July.The region receives 320 mm of precipitation annually with 200 mm falling between June and August (Toolik Environmental Data Center Team, 2014).Feniak Lake is located 360 km west of the Toolik Field Station in the central Brooks Range at the northeast boarder of the Noatak National Preserve.The mean annual temperature is -7°C (Jorgenson et al., 2008) and mean precipitation is 450 mm (WRCC, 2011).The Kelly River Ranger Station is located on the western boarder of the Noatak National Preserve, 170 km west of Feniak Lake.Average annual temperature is -5.4°C and the area receives a mean of 300 mm of precipitation, a third of which falls during the growing season (Stottlemyer, 2001).
Vegetation is typical of Arctic tundra across the study region and includes moist acid tundra characterized by the tussock-forming sedge Eriophorum vaginatum, moist non-acidic tundra, and shrub tundra (Bhatt et al., 2010;Walker et al., 1998), with isolated stands of white spruce (Picea glauca) near the Kelly River Ranger Station (Sullivan and Sveinbjornsson, 2010).All three areas occur in bioclimate subzone E, the warmest region in the continuous permafrost zone (Walker et al., 2010).The foothills of the Brooks Range have been affected by multiple glaciations starting in the late Tertiary and continuing to 11 ka B.P. (Hamilton, 2003).Repeated rounds of glacial advance and retreat have resulted in a patchwork of glacial till, bedrock, and loess parent materials of various ages (Hamilton, 2010).Time since last glaciation can be associated with ecosystem properties including pH, organic layer depth, nutrient pools, vegetation community, and biogeochemical rates (Epstein et al., 2004;Hobbie et al., 2002;Lee et al., 2011;Walker et al., 1998).

Experimental design and sampling
To test our hypotheses concerning the intensity and duration of thermokarst impacts on aquatic chemistry, we sampled thermokarst features in all stages of development across landscape ages and vegetation types.We collected water from 83 thermokarst outflows and 61 adjacent undisturbed water bodies such as water tracks and first-order streams (22 locations did not have a suitable paired reference site).To quantify the evolution and duration of thermokarst effects through time, we classified features on a 0 -3 index (Fig. 3) following the development of a hypothetical feature from before initiation (0) to after stabilization (3).Development stages were defined as follows: 0. no apparent present or past thermo-degradation, 1. active thermodegradation (>25 % of headwall is actively expanding) with completely turbid outflow, 2. moderate thermo-degradation (<25 % of headwall is expanding) with somewhat turbid outflow, and 3. stabilized or limited thermo-degradation with complete or partial revegetation and clear outflow.Ideally we would have tested for trends in elemental export based on absolute feature age rather than a development stage proxy.However, identifying a reliable time since formation requires high-resolution remote sensing or radiocarbon dating (Kreiger 2012, Balser and Jones 2014, Pizano et al., 2014) which was beyond the scope of our study given the large number of features sampled.Features were classified in the field prior to any chemical analyses, precluding the possibility of bias in classification based on chemical signature.We also performed a sensitivity analysis, randomly excluding a third of the thermokarst features, to explore robustness of the classification, which did not substantively change the results or interpretation.
Vegetation class was determined in the field and cross-referenced with published vegetation maps when available (Walker et al., 2005).Glacial geology and surface age were based on recent maps of the study region (Hamilton, 2010(Hamilton, , 2003;;Kanevskiy et al., 2011).Most site ages ranged from 10-200 ka, though six sites occurred on surfaces unglaciated for more than 1000 ka.We classified sites on surfaces younger than 25 ka as young, and sites over 50 ka as old, corresponding to the split between the Itkillik I and II advances (Hamilton, 2003).
Samples for carbon and nutrient analysis were filtered in the field (0.7 µm effective pore size, Advanctec GF-75) into 60 ml high density polyethylene (HDPE) bottles, except when excess sediment required settling overnight when samples were filtered within 24 hours.After filtration, samples were frozen until analysis.We measured DOC and dissolved inorganic carbon (DIC) with a Shimadzu TOC-5000 connected to an Antek 7050 chemiluminescent detector to quantify total dissolved nitrogen after combustion to NOx.We analyzed major ions (NO3 -, NO2 -, SO4 2-, Cl -, NH4 + , Ca 2+ , Na + , Mg 2+ , and K + ) on a Dionex DX-320 ion chromatograph.We calculated dissolved organic nitrogen (DON) by subtracting dissolved inorganic nitrogen (DIN = NO3 -+ NH4 + + NO2 -) from total dissolved nitrogen, and we calculated the DOC to DON ratio (C:N) of dissolved organic matter, an indicator of organic matter source and degree of prior processing (Amon et al., 2012).To determine the percentage of thermokarst outflow coming from ground ice, we analyzed δ D and δ 18 O on a Picarro L1102i via cavity ringdown spectroscopy.
Because lateral fluxes of dissolved gas can constitute a considerable portion of Arctic carbon budgets (Kling et al., 1992;Striegl et al., 2012), we measured dissolved CO2, CH4, and N2O in feature outflows and reference water.At each site we collected a 30 ml sample of bubble-free water in a 60 ml gas-tight syringe accompanied by an ambient atmospheric sample in a 15 ml evacuated gas vial.Upon return to the lab or camp we added 30 ml of atmosphere to the syringe and shook vigorously for two minutes to facilitate equilibration of dissolved gases with the introduced headspace, and then injected a sample of the headspace into an evacuated gas vial for storage until analysis.We determined CO2, CH4, and N2O concentration of the headspace sample on a Varian 3300 gas chromatograph with a flame ionization detector and methanizer for carbon species and an electron capture detector for N2O.We calculated the proportion of total gas dissolved in solution and in the headspace using Henry's constants adjusted for extraction temperature (Wilhelm et al., 1977), and subtracted ambient gas introduced during extraction to determine initial concentration.We calculated saturation as the percent of equilibrium water concentration based on atmospheric partial pressure and water temperatures at the time of sampling and extraction.
To determine the direct contribution of carbon and nutrients from ground ice, we sampled exposed headwall ice at 24 sites.Because gravel and cobbles prevented motorized coring, we collected ice scrapings with a hand corer into Ziplock™ bags, which we filtered and analyzed after melt as previously described.We compared concentrations of carbon and nutrients in ground ice to feature outflows to determine whether solute concentrations were changing as water flowed through the feature.At these sites we used the difference between the δ 18 O of ground ice and adjacent reference water (stream or water track) to determine the proportion of outflow contributed by ground ice.We calculated the proportion from ground ice with a simple two end-member model: Fraction from ground ice = (δ 18 Oout − δ 18 Osw) / ( δ 18 Oice − δ 18 Osw) (Equation 1) where out is feature outflow, sw is undisturbed surface water, and ice is headwall ice.
To convert carbon and nitrogen concentrations into elemental loads and areal yields we measured discharge at the outflow of 26 thermokarst features using salt-dilution gauging (Figs. 1 and 3).We logged electrical conductivity with a YSI Professional Plus conductivity meter and added 10-100g of dissolved NaCl upstream of the probe by 10-20 m, depending on the size of the outflow.Discharge was determined by total dilution of the tracer as it passed by the probe (Wlostowski et al., 2013).We mapped feature perimeters with a commercial-grade, handheld GPS, except for four sites around Toolik, which were mapped by the Toolik Field Station GIS staff with a survey-grade GPS and base station.We determined areal DOC and DIN daily yields from the gauged sites by multiplying outflow concentration by discharge and dividing by the area of the feature.For sites with surface water flowing into the top of the feature (primarily gullies but also some ALDS and slumps) we subtracted the reference concentration of each solute before calculating yield, so the estimate represented only the contribution from the area disturbed by thermokarst (assuming that any unmeasured lateral water inputs along the feature margin have the same carbon and nitrogen concentrations).
Because determining contributing area for reference sites was not possible due to the low resolution of digital elevation models for the study region, we compared yields from thermokarst features to published yields of DOC and DIN from upland Arctic tundra (Giesler et al., 2014;McClelland et al., 2007;McClelland et al., 2014;Olefeldt et al., 2013;Peterson et al., 1993;Peterson et al., 1986;Townsend-Small et al., 2011).

Statistical analyses
We used a linear mixed-effects model to test for effects of thermokarst development stage, feature type, vegetation, and landscape age on water chemistry while accounting for spatial and temporal non-independence in the data.For each water chemistry parameter we used a mixedeffects analysis of variance (ANOVA) with development stage crossed with feature type and vegetation and landscape age as fixed effects.We included site as a random effect to pair thermokarst outflows with their adjacent reference water.The models included seasonal and inter-annual variability both across and within sites.We visually inspected residual plots for deviations from normality and homoscedasticity, and transformed response and predictor variables when necessary.We simplified the full model by automated backwards elimination, using restricted maximum likelihood to evaluate fixed effects and likelihood ratio tests for random effects.To test for differences between groups, we performed post-hoc Tukey honest significant difference tests on the least squares means using Satterthwaite approximation to estimate denominator degrees of freedom.We used Pearson product-moment correlation to test for associations between water chemistry parameters and development stage, which we recoded low to high and treated as a continuous variable of disturbance intensity.A decision criterion of α = 0.05 was used for all tests.
All analyses were performed in R 3.0.2(R Core Team, 2013) with the lme4 and lmerTest packages (Bates et al., 2013;Kuznetsova et al., 2014).The complete dataset is available through the Advanced Cooperative Arctic Data and Information Service at www.aoncadis.org/dataset/soil_water_gas_alaskan_thermokarst.html.

Thermokarst distribution and characteristics
Feature types were not distributed equally among vegetation classes with most active-layer detachment slides occurring on non-acidic tundra, most thermo-erosion gullies occurring on acidic tundra, and thaw slumps distributed among tundra types (Table 1).Feature types were also unevenly distributed between development stages with over half of slumps classified as stage 1 (very active) compared to approximately 30 % of slides and gullies.Over 90 % of all features were associated with, or intersected a water body (Table 1).Slides and gullies occurred primarily on or next to water tracks or headwater streams and the majority of thaw slumps were on lakeshores.Slides tended to occur in the highest topographic positions, slumps were distributed across high and low gradient surfaces, and gullies were most common on foot slopes or valley bottoms.
Discharge from thermokarst features varied widely by feature type and individual features in the study, from no flow at some stabilized slumps and slides to 9.4 L sec -1 at one slide (Table 1).Mean discharge was highest for slides and lowest for slumps.For sites where we estimated the proportion of outflow derived from ground ice, the ice contribution varied from 0-97 %.
Slumps had the highest average ground ice contribution and slides had the lowest, though these values are not representative of all features, since they are only based on sites with exposed ground ice.Generally sites with high discharge (> 2 L sec -1 ) had little contribution from ground ice, except several large slumps with very active headwall retreat.

Effects of development stage and morphology on water chemistry
Thermokarst significantly altered concentrations of carbon, nitrogen, and other solutes but the magnitude and duration of these effects differed by feature type (Figs. 4,5,and 6).For most parameters, effects were largest at the most active features, with differences tapering off as activity decreased.However, DOC in slide outflows as well as DIC, Mg 2+ , Ca + , and dissolved N2O concentrations in gully outflows were highest in stabilized features.Slumps tended to have the largest effect on solute concentrations.For example, SO4 2-concentration was 30-fold higher than reference in stage-1 outflows, compared to 3.3-and 1.5-fold higher for gullies and slides, respectively.Gully reference and outflow chemistry was generally distinct from slides and slumps, with higher dissolved gas concentrations and DOC:DON, but lower concentrations of ions and DIC.
Thaw slumps caused the greatest increase in dissolved organic matter concentration, with DOC and DON 2.6-and 4.0-fold greater in stage-1 features, compared with 1.6-and 1.4-fold increases in slides, and 2.2-and 1.6-fold increases in gullies of DOC and DON, respectively (Fig. 4).Thermokarst had a much larger impact on inorganic nitrogen, with mean NH4 + and NO3 -concentrations 9-to 27-fold greater in stage-1 features (Fig. 5).Consequently, the relative proportion of DIN, which made up less than 10 % of total nitrogen in reference waters, constituted 26 to 38 % of total nitrogen in stage-1 features and 48 % of total nitrogen in stage-2 gullies (Fig. 7).NH4 + was the dominant form of DIN for all feature types and development stages except stage-3 (stabilized) slides where NO3 -made up 70 % of DIN.Elevated DIN persisted through stage 2 for slumps and through stabilization for gullies.
Dissolved CH4 concentration was 92 % and 89 % lower than reference for stage-1 gullies and slumps, respectively (Fig. 4).However, there were no significant differences by development stage for dissolved CO2 and dissolved N2O was only significantly elevated in stabilized gullies.
Across all development stages and feature types, 93 % and 97 % of all samples were supersaturated with CO2 and CH4, respectively, whereas 51 % of samples were supersatured with N2O.
Specific yields of DOC and DIN from stage-1 thermokarst features were 30-and 57-fold higher than literature values for undisturbed tundra, respectively (Fig. 8).The geometric mean yield for level-1 features was 0.45 g C m -2 day -1 for DOC and 3.8 mg N m -2 day -1 for DIN, though there was considerable variability between individual sites within activity levels.DOC and DIN yields from stabilized features were within the range of literature values for undisturbed tundra.Yields varied more strongly by activity level than by feature type, with similar yields from the most active ALDS, gullies, and slumps.
For the five sites with repeated measures of thermokarst outflow chemistry, solute concentrations were variable between samplings but did not show systematic seasonal or interannual trends, except for DIC concentration and δ 18 O, which both increased through the growing season (Figs.S1 -S8).

Ground-ice, vegetation, and landscape age
Permafrost ice was high in dissolved carbon, nitrogen, and solutes and had a depleted δ 18 O signature relative to reference waters (Table 3).Average concentrations of DIC, NH4 + , and K + were higher in ground ice than feature outflow, indicating uptake or dilution during transport from the feature headwall to outflow.However, all other solutes, notably DOC, NO3 -, and SO4 2- , were higher in outflows than in ground ice, indicating net production or contribution from soils or more concentrated flowpaths during transit.
Landscape age modulated the effect of upland thermokarst on water chemistry, with much larger differences between impacted and undisturbed concentrations of DOC, NH4 + , Cl -and SO4 2-at sites occurring on surfaces older than 50 ka (Fig. 9).Vegetation had a smaller effect on fewer parameters with only DOC, Ca + , and Cl -differing significantly by vegetation community independent of development stage, feature type, and landscape age, with different patterns between vegetation communities for each solute (Fig. 10).

Discussion
There is conflicting evidence of the impacts of upland thermokarst on concentrations and fluxes of DOC, nutrients, and other solutes (Bowden et al., 2008;Thompson et al., 2012), as well as the intensity and duration of these effects (Kokelj et al., 2005;Lafreniere and Lamoureux, 2013;Thienpont et al., 2013).Our spatially extensive sampling of active and stabilized features revealed that upland thermokarst consistently increases DOC and other solute concentrations with a particularly large effect on inorganic nitrogen.Magnitude and duration of thermokarst effects on water chemistry differed by feature type and secondarily by landscape age.Most solutes returned to undisturbed concentrations after feature stabilization, but elevated inorganic nitrogen and several other parameters persisted in gully and slump outflows, suggesting these feature types could have long-lasting impacts on aquatic nutrient dynamics.

Patterns of carbon and nitrogen release from upland thermokarst
We hypothesized that thermokarst would increase or decrease DOC concentration in surface waters depending on the balance of DOC production and removal processes active during feature formation.Despite large organic layer losses and abundant exposed mineral soil (Pizano et al., 2014), upland thermokarst significantly increased average DOC concentration and yield for all feature types.Additionally, DOC from active thermokarst features is three to four times more bio-and photo-degradable than active-layer-derived DOC (Abbott et al., 2014;Cory et al., 2013) changing the implications of this release at different spatial scales.DOC mobilized by thermokarst is likely to be mineralized rapidly in receiving soils, streams, and lakes, accelerating transfer of permafrost carbon to the atmosphere (Vonk et al., 2013) but reducing the impact of this disturbance on estuaries of the Arctic Ocean (McClelland et al., 2012;Striegl et al., 2005).
Upland thermokarst had a relatively larger effect on aquatic nitrogen than carbon concentrations, reducing the C:N ratio of dissolved organic matter and causing substantial and long-lasting release of inorganic nitrogen.Phosphorus, not nitrogen, is typically the most limiting nutrient in Arctic freshwater systems (O'Brien et al., 2005;Slavik et al., 2004), however, nitrogen and silica limit productivity in Arctic estuaries and the Arctic Ocean (McClelland et al., 2012;Vancoppenolle et al., 2013).If thermokarst nitrogen release is accompanied by bioavailable phosphorus, more nitrogen will be retained in inland aquatic ecosystems, whereas if thermokarst outflows have relatively little phosphorus, a larger proportion of liberated nitrogen will reach the ocean.Thermokarst can increase phosphorus loading (Bowden et al., 2008;Hobbie et al., 1999), but the relative impact of upland thermokarst on nutrient stoichiometry remains an important unknown.
Along with changes in solute concentrations and characteristics, upland thermokarst may affect the seasonality of solute flux.For most aquatic ecosystems in the Arctic, the majority of annual carbon and nutrient load occurs during snowmelt or early spring (Holmes et al., 2012).While carbon and nitrogen concentrations in thermokarst outflow do not appear to vary systematically over the season, thermokarst discharge, which depends primarily on air temperature and net radiation, peaks in mid to late summer (Kokelj and Jorgenson, 2013;Lantuit and Pollard, 2005;Lantz and Kokelj, 2008).Late-season delivery of carbon and nitrogen would have a larger relative impact on surface water concentrations, further modifying functioning of Arctic rivers and lakes.This shift could also affect Arctic estuaries, where nutrients and carbon are taken up quickly during open-water season but transported to the Arctic Ocean during ice cover (Townsend-Small et al., 2011).
Feature morphology strongly influenced magnitude and duration of thermokarst effects on water chemistry, with slides having a smaller and shorter impact than gullies or slumps.This could be due to differences in feature depth and duration of feature growth.In permafrost soil, leachable solutes are typically highest below the transition layer at the top of the permafrost table (Keller et al., 2007;Kokelj and Burn, 2003;Malone et al., 2013) and the age and characteristics of soil carbon differ strongly with depth (Guo et al., 2007;Neff et al., 2006;Nowinski et al., 2010;Schuur et al., 2009).Shallow slides are less likely to expose deeper, solute-rich soils to hydrologic export than slumps and gullies, which cut meters into permafrost.
However, slides caused a similar magnitude of increase as gullies and slumps for inorganic nitrogen concentration, suggesting that altered dynamics at the surface rather than depth of disturbance may determine nitrogen available for export.For all feature types, effects on carbon, nitrogen, and other solutes were largely limited to the period of active feature formation, meaning that the influence of upland thermokarst is directly related to period of active growth.In this regard slides, gullies, and slumps are dramatically different.Slides typically form suddenly, over a period of weeks, days, or even hours (Lewkowicz, 2007) and stabilize the same season they appear (Lafreniere and Lamoureux, 2013).In contrast, large thaw slumps commonly remain active for 12-50 years (Burn, 2000;Kokelj et al., 2013;Lewkowicz, 1987) though small slumps stabilize in less than ten years (Kokelj et al., 2009).
Less is known about gully longevity, but based on average feature size and rates of headwall retreat, they remain active for five to ten years (Jorgenson and Osterkamp, 2005), with large features lasting over a decade (Godin and Fortier, 2012).Differences in outflow chemistry between feature types agree with findings from high Arctic systems suggesting that slide formation may have relatively limited impact on water chemistry (Lewis et al., 2012), and suggest that gullies and slumps, with their long active periods and influential position in hydrologic networks (Krieger, 2012), are likely to have a persistent and widespread effect on aquatic ecosystems.

Decrease in dissolved methane
There are several possible mechanisms behind the unexpected 90 % decrease of dissolved CH4 in gully and slump outflows.Greater thaw-depth within features could facilitate infiltration, creating a larger aerated zone where CH4 oxidation can occur (Schuur et al., 2009).
Slides may have had no effect on dissolved CH4 because they do not affect thaw depth as profoundly as gullies and slumps.However soils affected by slides, gullies, and slumps have partial pressures of CH4 higher or equal to reference tundra (Abbott and Jones, 2015), suggesting that low CH4 in thermokarst outflows is due to changes in production or consumption in the water column, rather than in soils.For slumps this decrease may be due to high concentrations of SO4 2-released during thermokarst formation.SO4 2-is an energetically favorable electron acceptor compared to the low molecular weight organic compounds or CO2 used by methanogens (Dar et al., 2008), and sulfate-reducing bacteria can inhibit methane production through competition for molecular substrates (Muyzer and Stams, 2008).SO4 2- concentration was negatively associated with dissolved CH4 across site types and development stages, further supporting this hypothesis.However, SO4 2-release does not explain decreased dissolved CH4 in gully outflows since we observed no change in gully SO4 2-.One possibility is that high inorganic nitrogen concentration is stimulating CH4 consumption in gully and slump outflows.While elevated DIN can suppress high-affinity methanotrophs responsible for CH4 oxidation in low-CH4 environments, DIN can stimulate consumption by low-affinity methanotrophs that dominate consumption in high CH4 environments and are commonly nitrogen-limited (Bodelier and Laanbroek, 2004).This would explain the large CH4 decrease in gully outflows where CH4 concentration was high, and the lack of response in slide outflows where CH4 was 10-fold lower despite similar changes in DIN concentration.
Similar concentrations of SO4 2-have been observed in outflows of thaw slumps in the Mackenzie delta (Kokelj et al., 2005;Malone et al., 2013) and there is evidence of enhanced sulfur availability in lakes throughout the Arctic (Drevnick et al., 2010).The widespread release of SO4 2-from upland thermokarst may have important implications for carbon cycling as the permafrost region thaws.Increases in freshwater SO4 2-could accelerate anaerobic decomposition of organic carbon liberated from permafrost (Einsele et al., 2001) and suppress CH4 production after permafrost thaw, modulating one of the key feedbacks from the permafrost system on global climate (Walter et al., 2006).

4.3
Where is thermokarst nitrogen coming from?
Though primary production in high-latitude terrestrial ecosystems tends to be limited by nitrogen, suggesting that bioavailable forms of nitrogen should be retained (Vitousek and Reiners, 1975), there are numerous reports of inorganic nitrogen loss from landscapes affected by permafrost degradation (Jones et al., 2005;Mack et al., 2004;McClelland et al., 2007).
Contrary to our hypothesis that high demand for nutrients by re-establishing plants would decrease nutrient concentrations in thermokarst outflows during recovery, NH4 + concentration was elevated in stabilized gullies and in no case was DIN significantly lower in recovering features than in undisturbed tundra.This suggests that either nitrogen is not limiting plant growth during revegetation or pathways of nitrogen loss bypass locations of high uptake (e.g.preferential flowpaths below plant rooting zones).
Microenvironments in thermokarst can favor deciduous shrub establishment including nitrogen-fixing species (Lantz et al., 2009), a potential source for thermokarst nitrogen.
However, even in the absence of nitrogen-fixing species, surface soils in recovering thermokarst features re-accumulate nitrogen rapidly (Pizano et al., 2014).Upland thermokarst can warm wintertime soil temperature by up to 6°C due to conductive heat flux to soils during summer and added insulation in winter from deeper snow (Burn, 2000).If nitrogen mineralization continues through the fall and winter in thawed soils beneath thermokarst scars, hydrologic activity in the spring or deep shrub roots could transport inorganic nitrogen to the surface, fueling productivity and hydrologic export.The isotopic signature of NO3 -draining a high Arctic catchment impacted by upland thermokarst suggests DIN from thermokarst is derived from the heterotrophic decomposition of organic matter found in the mineral soil (Louiseize et al., 2014), supporting this hypothesis.Additionally or alternatively, a portion of inorganic nitrogen in upland thermokarst outflow may come from mineralization of labile dissolved organic matter in the water column or soil solution.This would explain the strong correlation between DIN concentration and DOC biodegradability observed in several Arctic and boreal ecosystems (Abbott et al., 2014;Balcarczyk et al., 2009;Wickland et al., 2012).

Shifts in landscape-scale water chemistry
As high latitudes warm, ecosystems are experiencing widespread shifts in aquatic chemistry including an increase in DOC flux in areas with peat and thick organic soil (Frey and McClelland, 2009), a decrease in DOC where organic soil is shallow (McClelland et al., 2007;Petrone et al., 2006;Striegl et al., 2005), increases in major ion concentrations (Frey and McClelland, 2009;Giesler et al., 2014;Keller et al., 2010), and increased inorganic nutrient flux (Jones et al., 2005;McClelland et al., 2007;Petrone et al., 2006).These changes in catchment-scale solute fluxes have primarily been attributed to mechanisms associated with gradual thaw such as deepening of surface flowpaths and changes in residence time.However, thermokarst may also be contributing to these shifts in catchment-scale chemistry (Frey and McClelland, 2009).The chemical signature of dissolved organic matter from thermokarst closely matches biodegradable DOC recently detected in boreal streams and rivers (Abbott et al., 2014;Balcarczyk et al., 2009;Wickland et al., 2012) and increases of DIN and solutes from thermokarst match circumpolar changes attributed to a shift towards greater ground-water inputs (Frey and McClelland, 2009;Frey et al., 2007).
Currently a scarcity of observations of the spatial extent and distribution of upland thermokarst features and the annual elemental yields for different feature and landscape types limits our ability to evaluate the relative importance of gradual thaw and thermokarst in determining the evolution of high-latitude biogeochemistry.Though our estimates of DOC and DIN daily yield are based on individual measurements from a relatively small set of features, if elemental yields from upland thermokarst are similar to the range observed here, this spatially limited disturbance may have a large influence on landscape-level carbon and nitrogen fluxes.A simple scaling exercise based on projections of permafrost degradation, average feature lifetimes, and daily yields measured here suggests that though upland thermokarst is only expected to directly impact 3% of the total circumarctic watershed by 2100, it may cause a 2.7 -23% increase in annual circumarctic DOC flux and a 2.2 -19% increase in dissolved inorganic nitrogen averaged over 2050 -2100 (see Table S1 for assumptions).While these fluxes are highly speculative, they underline the potential of this spatially limited disturbance to influence the rate of carbon and nitrogen release from thawing permafrost.

Conclusions
Upland thermokarst across the foothills of the Brooks Range caused substantial increases of inorganic nitrogen, DOC, and other solute concentrations.Thaw slumps and thermo-erosion gullies had larger impacts on solute concentrations and are likely more important than slides to surface water chemistry because they can remain active for multiple years.The delivery of labile carbon and nutrients such as SO4 2-and inorganic nitrogen to downstream or downslope ecosystems could have important consequences for offsite carbon cycling, accelerating decomposition of organic matter in anoxic environments and priming the decomposition of recalcitrant organic matter.The fact that individual features can impact entire lakes or river reaches over multiple years in combination with the large portion of the landscape underlain by ice-rich permafrost suggest that upland thermokarst may be the dominant disturbance affecting aquatic ecosystems as the Arctic warms.

Author contributions
Abbott and Jones designed the experiment and worked closely on the manuscript written by Abbott.Abbott, Godsey, and Larouche carried out sample collection, preparation, and analysis.All authors helped refine the experimental design and provided input on the manuscript.(Giesler et al., 2014;McClelland et al., 2007;McClelland et al., 2014;Olefeldt et al., 2013;Peterson et al., 1993;Peterson et al., 1986;Townsend-Small et al., 2011).See Table 1 for complete definition of development stages: 0=reference, 1=most active, and 3=stabilized and were located in three areas of upland tundra underlain by continuous permafrost in the foothills of the Brooks Range.We collected samples during the growing season (June-August) of 2009-2012 and May of 2011 in the region surrounding the Toolik Field Station, with additional sampling in the Noatak National Preserve near the Kelly River Ranger Station in 2010 and Feniak Lake in 2011.Most sites were sampled a single time over the course of the study, except for the five most accessible features near the Toolik Field Station, which we sampled 2 -4 times each summer from 2009 -2011.

Figure 2 .
Figure 2. Map of study area.Features near the Kelly River Ranger Station were sampled in July of 2010, near Feniak Lake field camp in July of 2011, and near the Toolik Field Station in the summers of 2009-2012.

Figure 3 .
Figure 3. Aerial view of three thaw slumps in various development stages.Samples were collected from feature outflows and adjacent water bodies such as unimpacted water tracks, streams, and lakes to assess the impact of thermokarst on water chemistry.

Figure 4 .
Figure 4. Dissolved carbon species and characteristics in outflow from 22 active-layer detachment slides, 19 thermo-erosion gullies, 42 thaw slumps, and 61 reference features in upland tundra on the North Slope of Alaska.Open circles signify statistical difference from stage-0 undisturbed sites, α = 0.05.Different letters above panels represent significant differences between feature types.Error bars represent SE estimated by mixed-effectsANOVA after accounting for between-site variability.See Table1for complete definition of

Figure 9 .
Figure 9. Mean (±95 % CI) of parameters that varied significantly by surface age.Impacted and reference concentrations are shown independently when the interaction between thermokarst impact and surface age was significant, otherwise, results are combined.

Figure 10 .
Figure10.Mean (±95 % CI) of DOC for thermokarst outflow and reference water from sites occurring on moist acidic (MAT), non-acidic (MNAT), and shrub tundra.Different letters above panels represent significant differences between feature types.DOC, Ca 2+ , and Cl - were the only parameters for which vegetation was a significant predictor in the mixedeffects ANOVA when accounting for development stage, feature type, and landscape age.

Table 1
Figure 7.The relative proportion of carbon and nitrogen species in thermokarst outflow by feature type and development stage.See Figs. 2 and 3 for estimates of error and statistical tests for each parameter andTable 1 for complete definition of development stages: +. Symbology the same as Fig.2.

Table 1 .
Characteristics of upland thermokarst features in study.SE) characteristics of upland thermokarst on the North Slope of Alaska.*Sample size for discharge and ground ice contribution measurements, **Sample size for landscape position and development stage.Development stages were defined as follows: 0. No apparent present or past thermo-degradation, 1. active thermo-degradation (>25 % of headwall is actively expanding) with completely turbid outflow, 2. moderate thermo-degradation (<25 % of headwall is expanding) with somewhat turbid outflow, 3. stabilized or limited thermodegradation with complete or partial revegetation and clear outflow.

Table 2 .
Correlations between water chemistry parameters for 83 thermokarst features and 61 reference water tracks and first order streams.

37 0.35 0.63 0.42
Strength of relationships was determined by Pearson product-moment correlation.Significant correlations (p < 0.05) are in bold.Relationships were visually inspected and transformed when necessary to meet the assumption of linearity (log and exponential transformations noted in the variable column on the left).All units are μM except dissolved gases (CO2, CH4, and N2O), which are ppmv, DOC:DON, which is a unitless ratio, and activity, which is recoded development stage 1-4 (low to high) treated as a non-parametric continuous variable.

Table 3 .
Water chemistry for ground ice, thermokarst outflows and reference waters.Mean (SE) water chemistry from ground ice, outflow, and reference water for the 5 slides, 3 1 gullies, and 16 slumps where we sampled ground ice exposed by thermokarst formation.2 *DOC:DON is a unitless ratio and δ 18 O is ‰. 3