Site description

The snow-scoured alpine tundra study site was located on Niwot Ridge in the Colorado Rocky Mountains, USA27. Niwot Ridge begins on the Continental Divide and extends eastward as a narrow arête for about 2 km before it widens into a series of gentle slopes on rounded knobs down to the alpine treeline, approximately 8–9 km east of the Continental Divide30,52. In total, the Niwot Ridge alpine zone covers approximately 20 km2 in area31. Over the period 1982 to 2012, the mean annual air temperature and precipitation near the study site were –2.2 °C and 884 mm (~75% in the form of snow), respectively, but snow cover was sparse due to windy conditions that scoured snow from the soil surface throughout the winter27,53. Continuous meteorological data have been collected on Niwot Ridge since the early 1950s, including the longest high-altitude air temperature record in the conterminous United States54. This makes Niwot Ridge an ideal place to study the effects of climate change and variability on the alpine tundra carbon cycle, and two identical 3-m eddy covariance (EC) towers, believed to be the highest in North America, were established near the T-Van site (40°03′06″ N; 105°35′07″ W; 3480 m asl) in 200726. These towers measured the NEE of CO 2 over a temporally consistent statistical measurement footprint wherein 80% of the cumulative flux originated from within 400 horizontal meters predominantly to the west of each tower along the prevailing wind direction (Fig. 1a)26,55. At each tower, a co-located sonic anemometer (CSAT 3; Campbell Scientific, Logan, UT, USA) and infrared gas analyzer (LI-7500; LI-COR, Lincoln, NE, USA) were used to quantify the vertical wind fluctuations and the density of atmospheric CO 2 , respectively26,27. The NEE was calculated as the covariance between instantaneous (10 Hz) deviations from the 30-min mean of the vertical wind speed and the scalar density of CO 2 . Mann-Kendall and least squares regression analyses were used to identify temporal trends in the cumulative NEE, and statistical correlation between meteorology and the carbon cycle, respectively.

Eddy covariance measurements

Previous studies have focused extensively on the data quality from the T-Van study area including the potential for systematic uncertainty resulting from advective flows, sensor heating, and/or beneath-tower CO 2 storage, which is generally reduced due to the short measurement height, gently sloping terrain, prevailing windy conditions, and persistent westerly wind direction at the study site26,27. Post-processing of the EC data consisted of standard coordinate rotation and the Webb adjustment56. Missing values in the eddy covariance data were replaced by the average value under similar meteorological conditions57, and the GPP and R eco were modeled following ref. 58. We primarily used data from the west eddy covariance tower since it was located ~50 m closer to the solifluction lobes that may correspond to permafrost. However, during multi-day periods of consecutively missing west tower data (5.5% of all data), empirical between-tower relationships specific to daytime and nighttime growing season NEE, and NEE during the rest of the year, were used for infilling purposes. These significant relationships (p ≪ 0.001) were developed using 1 year of simultaneous east and west tower data collection (data collected between November 2007 and October 2009), and validated using east tower measurements to predict the west tower data between 1 August 2013 and 31 December 2013. Using this approach, the modeled west tower cumulative NEE was 3.3% less than the observed west tower cumulative NEE during that time. There was no significant trend in the percentage of infilled values during the study period.

We used the average of three different methods to determine the growing season start and end dates and the variability around the resultant growing season length39. First, a carbon uptake period59 was calculated in which the growing season length corresponds to the total number of days in which the 24-h mean NEE was less than zero. Second, 5-day and 10-day running averages of NEE were evaluated and the growing season start and end dates were identified as the days in which the NEE crossed from positive to negative values and vice versa60. Third, variable length regressions were calculated on 3-day to 9-day subsets of total daily NEE during the spring and fall. For this method, the day with the steepest negative slope during the spring defined the beginning of the growing season and the day with the shallowest positive slope in autumn defined the end of the growing season61.

Radiocarbon measurements

Radiocarbon data are expressed as the parts per thousand (per mil; ‰) deviation from a standard (NBS Oxalic Acid I) of fixed isotopic composition:

$$\Delta {{}^{14}}{\mathrm{C}} = [{\mathrm{FM}} - 1]1000,$$ (1)

where

$${\mathrm{FM}} = \frac{\left(\frac{{\, }^{14}{\mathrm{C}}}{{\ }^{12}{\mathrm{C}}}\right){\mathrm{sample}}}{\left(\frac{{\ }^{14}{\mathrm{C}}}{{\ }^{12}{\mathrm{C}}}\right){\mathrm{standard}}},$$ (2)

and FM is the fraction modern. To correct for mass-dependent fractionation, all samples were corrected to a common δ13C value of −25‰ and also for radioactive decay between the time of sampling and measurement62. The δ13C was measured offline from CO 2 splits at the Institute of Arctic and Alpine Research (INSTAAR) Stable Isotope Laboratory at the University of Colorado Boulder63.

To collect soil respiration samples for radiocarbon analysis during the growing season, semi-permanent 25-cm diameter × 10-cm deep soil collars64 were buried to a depth of 5 cm at four alpine tundra sites including one dry meadow, one moist meadow, and two wet meadow locations on solifluction lobes. Soil respiration samples were collected by attaching a closed dynamic chamber system to the soil collar65; vegetation inside the soil collars was not clipped. Atmospheric (and soil) air from within the chamber was initially removed by circulating air (at 0.9 l min−1; rate controlled by an adjustable flow meter) from within the chamber headspace through a desiccant (Mg(ClO 4 ) 2 ) and then through a column filled with soda lime. The CO 2 was scrubbed in this way until the equivalent of approximately three chamber volumes had passed over the soda lime, in order to maximize the removal of the atmospheric radiocarbon signature. Sample-containing air was then diverted around the soda lime trap to pass through a 2-l glass flask and a CO 2 analyzer (LI-820; LI-COR, Lincoln, NE, USA) before being exhausted back into the chamber headspace (flow rate was reduced to 0.5 l min−1 during this step). The sample was isolated and collected by closing the valves on the flask when sufficient CO 2 had accumulated in the flask (~550 µmol mol−1) to prepare a graphite target for radiocarbon analysis.

Due to the combination of persistent high winds and small magnitude respiratory fluxes, this chamber system was unable to capture sufficient carbon for radiocarbon analysis during the winter. As a result, two types of gas wells were deployed to each of the four sampling locations for wintertime radiocarbon analysis. The first type of gas well was constructed of a length of 0.64-cm diameter stainless steel tubing buried vertically in the soil to 15 cm depth (with ~10 cm of tubing extending above the soil surface to allow for sampling), then capped with a 0.64-cm Swagelok union fitted with a Teflon septum66. The second type of gas well was constructed of a 7-cm length of 3.8-cm diameter polyvinyl chloride (PVC) pipe capped with a two-hole #10 rubber stopper and sealed with silicone gel34,67. For sampling, two lengths of 0.64-cm diameter Teflon tubing extended from the interior of this gas well (through the two-hole rubber stopper) to a height of approximately 5 cm above the soil surface. Samples were collected from each of these gas wells by connecting a 470 cc pre-evacuated stainless-steel canister to the aboveground tubing (the second piece of tubing from PVC gas wells remained tightly capped) via a length of custom-crimped stainless-steel tubing that allowed the canister to fill over the period of approximately 1 h, in order to minimize disturbance to the soil CO 2 profile. The average difference between simultaneous measurements of 14CO 2 collected during the growing season via the two gas well methods was 0.5‰ (n = 6). In winter, the 14CO 2 collected from the stainless-steel gas wells was an average of 6‰ higher than the 14CO 2 from PVC gas wells at the wet meadow site 17, but 19‰ lower than the PVC gas wells at the other wet meadow location (site 18) (n = 12). All air samples were purified and graphitized by the INSTAAR Laboratory for AMS Radiocarbon Preparation and Research at the University of Colorado Boulder, then shipped to the Keck Carbon Cycle AMS Lab at the University of California, Irvine, for analysis.

Although atmospheric CO 2 was removed from the chambers prior to sampling, we accounted for leakage of atmospheric CO 2 into the chambers during the accumulation period (which is of particular importance when soil respiration fluxes are low) by calculating the amount of ambient CO 2 (f air ) in each sample based on the δ13C ratios of CO 2 in the sample (δ13C obs ), ambient air (δ13C air ), and heterotrophic soil respiration (δ13C resp ):

$$f_{{\mathrm{air}}} = \frac{{\left( {{\rm{\delta}} \,{}^{\!\!13}{\mathrm{C}}_{{\mathrm{obs}}}} \right) - \left( {{\rm{\delta}} \,{}^{\!\!13}{\mathrm{C}}_{{\mathrm{resp}}}} \right)}}{{\left( {{\rm{\delta}} \,{}^{\!\!13}{\mathrm{C}}_{{\mathrm{air}}}} \right) - \left( {{\rm{\delta}} \,{}^{\!\!13}{\mathrm{C}}_{{\mathrm{resp}}}} \right)}},$$ (3)

where δ13C air = −8‰ and δ13C resp = −28‰68. Using this estimate of f air , we then calculated air-corrected 14C signatures (Δ14C cor ) from each observed soil respiration 14C signature (Δ14C obs ):

$$\Delta \,{}^{\!\!14}{\mathrm{C}}_{{\mathrm{cor}}} = \frac{{\Delta \,{}^{\!\!14}{\mathrm{C}}_{{\mathrm{obs}}} - \left( {f_{{\mathrm{air}}}\Delta \,{}^{\!\!14}{\mathrm{C}}_{{\mathrm{air}}}} \right)}}{{\left( {1\! - f_{{\mathrm{air}}}} \right)}}.$$ (4)

Since radiocarbon analysis of respired soil air yields only a bulk ∆14C value, the SOM was separated by physical means prior to radiocarbon analysis. Soils from 5 to 15 cm depth at the dry, moist, and wet meadow sites were air dried, sieved through 2 mm mesh, and then fractionated into inter-aggregate particulate organic matter (free light fraction (fLf)), particulate organic matter occluded within soil aggregates (occluded light fraction (oLf)), and organic matter that is complexed with minerals (heavy fraction; Hf) fractions (based on a density of 1.85 g cm−3) at the Lawrence Berkeley National Laboratory in Berkeley, California, USA25. All three soil fractions were similarly prepared for radiocarbon analysis at the INSTAAR Laboratory for AMS Radiocarbon Preparation and Research and shipped to the Keck Carbon Cycle AMS Lab for radiocarbon content analysis.

Soil carbon model

We used the software R to model the mean turnover time of the Lf (fLf) and Hf* (oLf + Hf) SOM fractions and their seasonal contribution to soil respiration from a dry meadow and two wet meadow sites on solifluction lobes. We chose these sites for modeling based on the availability of supporting information and contrasts in the respired 14CO 2 flux. For the wet meadow site, separate models were fit to the summer and winter respiration fluxes; winter respiration was not detectable at the dry meadow site precluding fitting of a separate model. In total, we assessed the output from three different model parametrizations (dry meadow, wet meadow summer, and wet meadow winter). Differences between the parameter fits of the summer and winter models and their implications are addressed in the manuscript text.

We used the two-pool model in series from the SoilR package as our basic modeling framework, with the Lf and Hf* data from the density separation as constraints on the partitioning of the total carbon (TC) between the two model pools. Each model was initialized at time = 0 with a single depth increment and TC = 12.1 kg C m−2, which was measured from depth-distributed sampling of the wet meadow site. We then ran the non-steady state models with annual time steps for the period from 1901 to 2015. Based on ref. 69, we applied estimated annual carbon inputs of 72 and 124 g C m−2 for dry and wet meadow simulations, respectively. We derived the radiocarbon content of annual inputs from an interpolation of pre- (IntCal13) and post-bomb (Northern Hemisphere zone 2)70 atmospheric values, without a time-lag.