Deeper snow increases the net soil organic carbon accrual rate in moist acidic tussock tundra: ²¹⁰Pb evidence from Arctic Alaska

ABSTRACT

The net change in the carbon inventory of arctic tundra remains uncertain as global warming leads to shifts in arctic water and carbon cycles. To better understand the response of arctic tundra carbon to changes in winter precipitation amount, we investigated soil depth profiles of carbon concentration and radionuclide activities (⁷Be, ¹³⁷Cs, ²¹⁰Pb, and ²⁴¹Am) in the active layer of a twenty-two-year winter snow depth manipulation experiment in moist acidic tussock tundra at Toolik Lake, Alaska. Depth correlations of cumulative carbon dry mass (g cm⁻²) vs. unsupported ²¹⁰Pb activity (mBq g⁻¹) were examined using a modified constant rate of supply (CRS) model. Results were best fit by two-slope CRS models indicating an apparent step temporal increase in the accumulation rate of soil organic carbon. Most of the best-fit model chronologies indicated that the increase in carbon accumulation rate apparently began and persisted after snow fence construction in 1994. The inhomogeneous nature of permafrost soils and their relatively low net carbon accumulation rates make it challenging to establish robust chronologic records. Nonetheless, the data obtained in this study support a decadal-scale increase in net soil organic carbon accumulation rate in the active layer of arctic moist acidic tussock tundra under conditions of increased winter precipitation.

KEYWORDS:

• Arctic Alaska
• tundra
• carbon
• ²¹⁰Pb
• snow depth

Introduction

Understanding the arctic carbon cycle today and its likely behavior in the future is essential to recognizing short- and long-term responses and feedbacks to the climate system. The impact of the changing Arctic on the climate system can be substantial because arctic ecosystems harbor as much as 50 percent of the total terrestrial soil carbon pool (Vitharana et al. 2017). Changes in soil organic carbon inventories in arctic regions are expected as a direct response to warmer summers (Chapin and Shaver 1985; Elmendorf et al. 2012; Gallego-Sala et al. 2018) and greater moisture availability (Zhang et al. 2013), including increased winter snow precipitation (Welker, Fahnestock, and Jones 2000; Blanc-Betes et al. 2016; Leffler et al. 2016; Jespersen et al. 2018). However, precipitation predictions are uncertain across the Arctic (some regions are predicted to have more and some are predicted to have less), and this variation is reflected in the annual fluctuations in snow cover and snow water equivalents reported for the past ten years by the annual National Oceanic and Atmospheric Administration Arctic Report Cards (see NOAA 2019).

Increases and decreases in snowfall and winter snowpack can lead to different ecosystem impacts (Bintanja and Selten 2014; Bintanja and Andry 2017). The effects of these changes on soil organic carbon inventories depend on soil carbon cycling processes (Schaeffer et al. 2013; Ricketts et al. 2016; Vitharana et al. 2017) and soil active layer depth (Pattison and Welker 2014; Jespersen et al. 2018). Improved understanding of the mechanisms driving the arctic carbon cycle continues to be at the forefront of research on the role of the “New Arctic” in the global climate system, especially processes that drive winter carbon emissions (Fahnestock et al. 1999; Lupascu et al. 2018; Natali et al. 2019).

Release of stored arctic carbon into the atmosphere, rivers, and ocean is becoming increasingly apparent (Czimczik and Welker 2010; Mishra et al. 2017; Nowinski et al. 2010; Lupascu et al. 2014, 2018; Csank et al. 2019). It is being driven by deepening of the active layer as permafrost thaws and by coastal and riverine erosion causing shoreline retreat and export of carbon-rich sediment and water into the Arctic Ocean (Martens et al. 2018). In contrast, arctic tundra is also becoming more productive as seen by increases in the Normalized Difference Vegetation Index values, due in part to greater leaf area associated with shrub density increases (Tape et al. 2012; Leffler et al. 2016; Mekonnen, Riley, and Grant 2018). The net effect of these accelerated carbon losses and productivity gains on the arctic carbon cycle is difficult to quantify because it requires a full picture of the annual CO₂ and CH₄ emissions (Blanc-Betes et al. 2016) and the export of dissolved organic and inorganic carbon by surface water flow (Lehn et al. 2017; Csank et al. 2019), as well as gains in soil organic carbon and biomass due to increases in net primary productivity (Welker, Fahnestock, and Jones 2000; Matamala et al. 2017). It is difficult to depict changes in carbon inventories in permafrost soils because we lack a baseline for net soil carbon accrual rates and how they may be affected by winter precipitation scenarios.

In this study, we explored a ²¹⁰Pb-based chronometric approach to obtain an integrative measure of the temporal change in net soil carbon accumulation rate in the active layer organic horizon of moist acidic tussock tundra in response to the construction of a snow fence near Toolik Lake, Alaska, at 68°38′ N, 149°38′ W; elevation 760 m (Figure 1). We measured and compared profiles of soil organic carbon and ²¹⁰Pb (half-life = 22.3 a) along with ⁷Be (half-life = 53.3 d), ¹³⁷Cs (half-life = 30.2 a), and ²⁴¹Am (half-life = 432.2 a) in a set of active layer cores from within and outside the area influenced by a twenty-two-year snow fence experiment. This was a reconnaissance study conducted with a relatively small set of cores to assess whether this approach could give meaningful results over a decadal timescale, given the known complications of ²¹⁰Pb-based soil chronology (Dörr and Münnich 1989, 1991; Kaste, Heimsath, and Bostick 2007; Kaste et al. 2011; Landis, Renshaw, and Kaste 2016; Olid et al. 2016; Gonzalez-Meler et al. 2018). Our objective was to ascertain the impact of deeper winter snow accumulation (caused by the snow fence) on the net rate of soil organic carbon accumulation. We assumed that all soil organic carbon in the active layer profile was produced from plant and microbial inputs and that ²¹⁰Pb deposited from the atmosphere could be used as a proxy to estimate the net carbon accumulation rate in the organic-rich horizon of the active layer. Key assumptions in this approach were (1) atmospheric ²¹⁰Pb deposited onto the tundra surface was fully retained within the active layer, (2) ²¹⁰Pb was strongly sorbed to organic matter where it has limited mobility (Landis, Renshaw, and Kaste 2014), and (3) soil organic matter decomposition rate was low in the cold arctic environment. Thus, the profile of ²¹⁰Pb deposited onto the organic carbon could allow a quantitative model of the net accumulation rate of soil organic carbon on the timescale of ²¹⁰Pb decay. We used ⁷Be as an indicator of near-surface ²¹⁰Pb mobility (Landis, Renshaw, and Kaste 2016) and made simple modifications to our CRS models to accommodate this effect and also to account for changes in net soil organic carbon accumulation rate and ²¹⁰Pb depositional flux. We compared our ²¹⁰Pb-based chronologies with the depth profiles of bomb-era fallout radionuclides ¹³⁷Cs and ²⁴¹Am.

Figure 1. Location map of the Toolik Field Station in northern Alaska. The inset shows the study area in relation to the circumpolar arctic region (graphic courtesy of Toolik Field Station)

Applications of fallout radionuclides in studies of sediment, soil, and peat accumulation (⁷Be, ¹³⁷Cs, ²¹⁰Pb, ²⁴¹Am) have been explored by a number of authors (Robbins and Edgington 1975; Dörr and Münnich 1989, 1991; Urban et al. 1990; Kaste, Heimsath, and Bostick 2007; Kaste et al. 2011; Abril and Gharbi 2012; Landis, Renshaw, and Kaste 2014, 2016; Shotyk et al. 2015; Olid et al. 2016; Davies et al. 2018; Gonzalez-Meler et al. 2018). The prior studies most relevant to the current study are those dealing with accumulation of soil and peat. There is a growing consensus reflected in the most recent studies that post-depositional downward mobility of ²¹⁰Pb into the accumulating organic horizon of soil and peat generally results in a nonexponential profile of ²¹⁰Pb activity with depth (e.g., Landis, Renshaw, and Kaste 2016; Olid et al. 2016). The downward motion of ²¹⁰Pb is attributed to advective transport of colloidal organic particles with adsorbed ²¹⁰Pb, and this generally results in maximum ²¹⁰Pb activity at a depth of several centimeters to tens of centimeters below the surface vegetation and litter layer (Landis, Renshaw, and Kaste 2016). The depth profile of ⁷Be is a good indicator of the depth of ²¹⁰Pb transport because these two radionuclides behave similarly in such systems (Landis, Renshaw, and Kaste 2016). Various numerical and analytical models have been proposed to accommodate near-surface ²¹⁰Pb mobility so that chronometric information can be obtained regarding mass accumulation rates and calendar dates as a function of mass depth. Most such models assume steady-state behavior with diffusion-like transport of ²¹⁰Pb in the near-surface (typically over a centimeter to decimeter depth range), and derived accumulation rates and calendar dates may be verified with the occurrence of distinct time horizons such as the peak activities of the bomb-pulse radionuclides ¹³⁷Cs and ²⁴¹Am that occurred in 1963–1964 (Kaste, Heimsath, and Bostick 2007; Landis, Renshaw, and Kaste 2016; Olid et al. 2016). Concordant chronologic results may be obtained between ²¹⁰Pb and other fallout radionuclides below the depth of maximum ²¹⁰Pb activity in soil profiles (Landis, Renshaw, and Kaste 2016).

Prior to this study, we anticipated that the consequences of deeper snow in winter could lead to lower soil carbon accumulation, due in part to higher rates of winter CO₂ emissions and accelerated summer decomposition (Schimel, Bilbrough, and Welker 2004), and that any changes in vegetation would be insufficient to counter these losses with soil carbon inputs via litter deposition, root growth and exudate inputs, and microbial growth. In contrast, simple models of our ²¹⁰Pb data from soil cores beneath artificially deepened winter snow are consistent with a substantial increase in the net rate of soil carbon accumulation relative to cores from ambient snow depth conditions outside the influence area of the snow fence. Our data also indicated a rapid response of soil organic carbon accumulation to increased snow depth, because the increase in net soil organic carbon accumulation rate in cores for which it could be reasonably quantified coincided with the calendar date of snow fence construction.

Methods

Field site

As part of the ITEX (International Tundra Experiment; Henry and Molau 1997; M. D. Walker et al. 1999), a 2.8 m-high × 60 m-long snow fence was built at the study location in 1994, perpendicular to the prevailing winter southerly winds, to artificially increase leeward snow accumulation to simulate the anticipated increase in winter precipitation with winter climate warming. Within the snow fence experimental plot, leeward snow drift accumulation zones were defined relative to ambient annual snow depth of ~0.5 m as deep (2–3 m), intermediate (0.5–2 m), and low (<0.5 m) (Figure 2; see also Supplementary Figure S1). The entire experimental plot has heterogeneous microtopography and an ~5° eastward slope. The active layer thickness of moist acidic tussock tundra around Toolik Lake averages ~30 cm during the three- to four-month growing season and reaches a maximum thaw depth of 45 to 50 cm by late August. The moist acidic tussock tundra is an ombotrophic bog, and the active layer consists of a 5- to 20-cm-thick organic layer overlying a mineral horizon derived from weathered glacial drift. Total annual precipitation is near 300 mm (~200 mm summer and ~100 mm winter), with snow cover lasting from September to May.

Figure 2. Experimental configuration showing cross section of snowdrift distribution leeward of the snow fence, with locations of core samples (deep, intermediate, low) and corresponding winter snow depths (m) in this study (modified from M. D. Walker et al. 1999)

Soil surface temperatures in the snow fence experimental plot reach 10°C to 15°C during summer months and decline to about −7°C in late winter under the deepest part of the snow drift, whereas temperatures as low as −30°C occur under ambient late winter conditions (M. D. Walker et al. 1999; Pattison and Welker 2014). During the early winter, before substantial snow accumulation, small amounts of leaves and windblown organic debris are deposited leeward within 25 m of the snow fence (deep snow zone; Figure 2; M. D. Walker et al. 1999). There is typically a two- to three-week delay in the timing of completion of the spring snowmelt in the deeper drift area leeward of the snow fence (M. D. Walker et al. 1999). Soils are acidic (pH < 5), coarse-loamy, and poorly drained. Vegetation is dominated by tussock-forming cottongrass sedge (Eriophorum vaginatum) and intertussock moss (Sphagnum sp.), with wet sedge (Carex bigelowii) and occasional shrubs including dwarf birch (Betula nana) and tealeaf willow (Salix pulchra), in a landscape that is typical of moist acidic tussock tundra across the Alaskan North Slope (M. D. Walker et al. 1999; Leffler et al. 2016). Vegetation cover across the snow fence experimental plot was generally homogeneous, except in the deep snow accumulation zone where the abundance of sphagnum moss and wet sedge was much higher and the abundance of cotton-grass sedge and dwarf birch was much lower (Miller 2014). Standing water is frequently present in the deep snow accumulation zone because it has a slight topographic depression.

Sample collection

Cores were collected in the intertussock areas dominated by sphagnum moss growth at the surface. Soil cores (4.8 cm diameter × 30 cm depth) were taken between 10 a.m. and 3 p.m. in early August 2016, including five cores from within the snow fence plot (one core from deep snow depth zone ~10 m leeward of the snow fence, three cores from intermediate snow depth zone ~35 m leeward of the snow fence, and one core from low snow depth zone ~65 m leeward of the snow fence) and three cores from outside the snow fence plot (intended as controls representing ambient snow conditions; Figure 3; see Supplementary Figure S2 for photographs of typical core sampling sites). Cores were immediately wrapped in Al foil, sealed in plastic bags, and kept frozen until prepared for analysis.

Figure 3. Location map showing collection sites within snow fence experimental plot for core samples measured in this study. Open circles indicate cores reflecting ambient winter snow conditions; filled circles indicate cores for which snow depth was affected by snow fence. Snow fence is approximately 600 m south of Toolik Lake at ~50 m elevation above lake level. Prevailing wind direction is southerly

Sample preparation

Frozen cores were sectioned by using a band saw in a −20°C cold room into 1-cm depth increments that were weighed before drying and then air-dried at 60°C to a constant weight for the direct determination of soil water content and bulk dry density (g cm⁻³) at the University of Illinois at Chicago. Core segments were homogenized after drying by using a high-speed rotary blender.

Gamma spectrometry

Aliquots (~2 to 3 g) of the homogenized dry samples were analyzed by gamma spectrometry at the University of Delaware’s Environmental Isotope Geochemistry Laboratory. Two high-purity Ge gamma spectrometers (model GWL-170-15-LB-AWT with 15-mm well diameter, EG&G Ortec, Ametek, Inc.) were used to measure gamma emissions of ⁷Be (477.2 keV), ¹³⁷Cs (661.7 keV), ²¹⁰Pb (46.5 keV), and ²⁴¹Am (59.5 keV).

Gamma spectrometer efficiencies were calibrated using measurements of certified radionuclide standards, including DL-1A (CANMET 2017) for ²¹⁰Pb and other U- and Th-series nuclides and NIST SRM-4357 (NIST 1997) for ¹³⁷Cs. We assumed that the efficiency of the ²⁴¹Am peak at 59.5 keV was equal to that determined for the ²³⁴Th peak at 63.3 keV, and we estimated the detector efficiency at the ⁷Be energy by linear interpolation between the nearest large peaks of the U- and Th-series nuclides in the DL-1A reference spectra. Specific activities and one-sigma counting errors were calculated using standard statistical techniques (Mook and DeVries 2001) and are reported in millibecquerels per gram dry mass (mBq g⁻¹). Reported activities were corrected for detector background and decay from the time of sampling to the start time of the gamma spectrometry measurement. Limits of quantification for radionuclides were defined as three times the standard deviation of the background under the peak used for the activity quantification, and values less than the limit of quantification are not reported.

Carbon concentrations

Total carbon contents (weight percentage, dry basis) were determined by using a zero-blank elemental analyzer (Costech Analytical, EC4010; Valencia, CA) at the University of Illinois at Chicago. Total carbon concentrations are assumed to equal total soil organic carbon, because inorganic carbon solid phases are unstable at the low pH of the moist acidic tundra active layer. No further acidification was performed during sample preparation. Precision of carbon concentrations was ±1 to 2 weight percentage, based on replicate analyses of samples and standards.

Calculation of carbon accumulation rates using ²¹⁰Pb data

Profiles of unsupported ²¹⁰Pb activity in soil cores were analyzed with a modification of the constant rate of supply (CRS) model (Appleby and Oldfield 1978; Robbins 1978), to obtain estimates of apparent net carbon accumulation rates (mg cm⁻² yr⁻¹) from best-fit linear regressions of cumulative carbon dry mass (mg cm⁻²) vs. ln(A₀/A_Z). Because the 1994 construction of the snow fence caused substantial increases in the annual accumulation of winter precipitation across the experimental site, we considered that ²¹⁰Pb deposition rates also could have increased leeward of the snow fence after 1994.

Sediment or soil mass accumulation rates can be calculated using the CRS model, assuming a constant flux of ²¹⁰Pb to the sediment, and this model allows changes in sedimentation rate to be identified graphically. In this case, the integrated activity of unsupported ²¹⁰Pb in the profile is expressed as

(1) ${\mathrm{A}}_{\mathrm{Z}}=\hspace{0.17em}{\mathrm{A}}_0{{\mathrm{e}}^{-}}^{\mathrm{l}\mathrm{t}}$(1)

where A_Z is the cumulative unsupported ²¹⁰Pb activity (Bq cm⁻²) below sediment depth Z

(2) ${\mathrm{A}}_{\mathrm{Z}}={\int }_{\mathrm{Z}}^{\mathrm{\infty }}\mathrm{C}\mathrm{z}\times \mathrm{p}\mathrm{z}$(2)

and A₀ is the cumulative excess ²¹⁰Pb (Bq cm⁻²) for the entire profile

(3) ${\mathrm{A}}_0={\int }_0^{\mathrm{\infty }}\mathrm{C}z\times \mathrm{p}\mathrm{z}$(3)

where λ is the decay constant of ²¹⁰Pb (0.03108 yr⁻¹), C_Z is the unsupported ²¹⁰Pb activity (mBq g⁻¹) in the core section at depth z, and ρz is the dry bulk density (g cm⁻³) of the core section at depth z. Unsupported ²¹⁰Pb activity was obtained by subtracting the mean value of supported ²¹⁰Pb as determined from deeper core sections exhibiting approximately constant total ²¹⁰Pb vs. depth. Supported ²¹⁰Pb activities in seven of the cores were calculated from four to ten deeper sections and ranged from 25 ± 4 to 37 ± 8 mBq g⁻¹. An anomalously high value of 54 ± 7 mBq g⁻¹ was determined for core INT-3 in which only two deeper sections were available for this calculation.

The approximate age t (years) of material at a given depth is given by

(4) $\mathrm{t}=1/\lambda \times ln\left({\mathrm{A}}_0/{\mathrm{A}}_{\mathrm{Z}}\right)$(4)

and the net mass accumulation rate (r) (mg cm⁻² yr⁻¹) is given by

(5) $r=\lambda M$(5)

where λ is the decay constant of ²¹⁰Pb (0.03108 yr⁻¹) and M is the slope of the line of cumulative dry mass (g cm⁻²) vs. ln(A₀/A_z). For the purpose of this study, we transformed cumulative dry mass to cumulative carbon mass by taking the product of the dry bulk density (g cm⁻³) and the weight percentage carbon for each 1-cm core segment and then integrating downward from top to bottom of the core profile. Simple modifications of the CRS model to account for step changes in carbon accumulation and ²¹⁰Pb flux are presented and discussed below with reference to the data obtained in this study.

Calculation of unsupported ²¹⁰Pb, ¹³⁷Cs, and ⁷Be inventories

The inventories of unsupported ²¹⁰Pb, ¹³⁷Cs, and ⁷Be in each core were calculated as follows:

(6) $I=\mathrm{\Sigma }{\rho }_i{A}_i\mathrm{\Delta }{z}_i$(6)

where I is the total inventory of each radionuclide (mBq cm⁻²), ρ_i is the bulk dry density of core segment i, A_i is the specific activity of each radionuclide in core segment i (mBq g⁻¹), and Δz_i is the thickness of core segment i (cm), and the sum is over all core segments.

Calculation of ²¹⁰Pb depositional flux

The depositional flux of ²¹⁰Pb for each core location is calculated as follows:

(7) $Q=\lambda I$(7)

where Q is the depositional flux of ²¹⁰Pb (mBq cm⁻² yr⁻¹), λ is the decay constant of ²¹⁰Pb (0.03108 yr⁻¹), and I is the unsupported ²¹⁰Pb inventory (mBq cm⁻²) as calculated by Equation (6).

Results and discussion

The analytical results of this study are compiled in Supplementary Information Table S1.

Depth profiles for specific activities (mBq g⁻¹) of unsupported ²¹⁰Pb, ⁷Be, ¹³⁷Cs, and ²⁴¹Am are shown in Figure 4. Soil organic C concentrations in cores from the study site ranged from a maximum of 40 to 47 weight percentage in the organic horizon (uppermost 4–15 cm of core). Vegetation (mainly sphagnum moss) and a minor amount of plant litter comprise the top 2 to 4 cm of most cores (Supplementary Figure S3), with the remaining organic layer consisting of variably compacted, partially humified plant material with minimal mineral content. Profiles of dry bulk density and weight percentage carbon vs. depth show that the organic layer has dry bulk density in the range of 0.05 to 0.2 g cm⁻³ and >30 weight percentage carbon (Supplementary Figures S4 and S5). Organic carbon decreases sharply below the organic layer in a 2- to 4-cm-thick transitional layer to about 2 to 5 weight percentage in the underlying fine-grained mineral horizon that has dry bulk density in the range from 0.8 to 1.4 g cm⁻³.

Figure 4. Activities (mBq g⁻¹) of unsupported ²¹⁰Pb, ⁷Be, ¹³⁷Cs, and ²⁴¹Am vs. depth (cm) in active layer soil core samples. Note different x-axis range of core profile DEEP-4

The activity of ⁷Be was measured only in five of the cores (those counted within a few months of collection). and in these cores it was found to have declining activity with depth in the upper 2 to 4 cm (Figure 4). The disappearance of ⁷Be with depth generally corresponded to the observed maxima of unsupported ²¹⁰Pb that ranged from 399 to 537 mBq g⁻¹ in the control (CTL) cores, 488 mBq g⁻¹ in the reduced snow depth (LOW) core, 476 to 681 mBq g⁻¹ in the intermediate snow depth (INT) cores, and 1,305 mBq g⁻¹ in the deep snow depth (DEEP) core. Unsupported ²¹⁰Pb exceeded supported ²¹⁰Pb in all cores by a factor of ten or more at the depth of maximum ²¹⁰Pb activity. The activities of ¹³⁷Cs were measured in all cores and each core had a well-defined maximum activity. Maximum activities occurred over a range of core segment depths from 0 to 1 cm up to 9 to 10 cm. Activities of ²⁴¹Am were detectable in all but one core profile.

Carbon accumulation rates determined from ²¹⁰Pb profiles

To obtain a net rate of soil organic carbon accumulation for each core profile, we evaluated the linear correlations of cumulative organic carbon (g cm⁻²) with ln(A₀/A_Z). We refer to this type of data analysis henceforth as a CRS plot. The slope of a best-fit linear regression between these two parameters in a typical CRS plot yields the net organic carbon accumulation rate. Changes in slope on a CRS plot may correspond to changes in net carbon accumulation rate. Carbon accumulation rates (mg cm⁻² yr⁻¹) derived from linear regression slopes are listed in Table 1 with errors, along with the total inventories of unsupported ²¹⁰Pb, ¹³⁷Cs, and ⁷Be and the steady-state equivalent depositional fluxes of ²¹⁰Pb.

Table 1. Carbon accumulation rates, radionuclide inventories, and ²¹⁰Pb fluxes.^a

	C accumulation rate^b(mg cm^−2 yr^−1)
Core	Pre-1994	Post-1994	Unsupported ^210Pbinventory (mBq cm^−2)	Steady-state equivalent ^210Pb flux(mBq cm^−2 a^−1)	^137Cs inventory(mBq cm^−2)	^7Be inventory(mBq cm^−2)
CTL-2	1.5 ± 0.1		105 ± 3	3.3 ± 0.1	83 ± 1	—
CTL-4	1.7 ± 0.1		79 ± 2	2.5 ± 0.1	57 ± 1	—
CTL-5	0.7 ± 0.1	2.8 ± 0.2	122 ± 3	3.8 ± 0.1	104 ± 1	0.8 ± 0.1
LOW-2	1.8 ± 0.3	4.6 ± 0.3	145 ± 3	4.5 ± 0.1	73 ± 1	7.9 ± 1.0
INT-1	1.1 ± 0.1	3.2 ± 0.1	146 ± 2	4.5 ± 0.1	70 ± 1	—
INT-3	—	—	135 ± 4	4.2 ± 0.1	45 ± 1	5.6 ± 0.1
INT-5	1.1^c	3.2 ± 0.1	138 ± 3	4.3 ± 0.1	54 ± 1	1.3 ± 0.1
DEEP-4	1.2^c	—	481 ± 5	15.0 ± 0.2	82 ± 1	16.2 ± 0.2

Notes. ^aErrors on inventories and flux are ±1 standard deviation based on counting statistics.

^bErrors on C accumulation rates are 95 percent confidence limits based on linear regressions. Range of R² values for regressions including more than two points is 0.973–0.999.

^cSlope is based on two-point regression; error in C accumulation rate not determined.

In order to assist our interpretation of the observed changes in slope observed in CRS regression lines for these cores, we did a simple numerical experiment (described in Supplementary Information Addendum and Table S2) to demonstrate the effects of step changes in arbitrary values of the net carbon accumulation rate and the ²¹⁰Pb depositional flux on the chronologic records obtained from CRS plots. We evaluated hypothetical changes in these parameters corresponding with the time of construction of the snow fence in 1994; these are portrayed in a demonstrative plot (Figure 5). Model A is a typical CRS model with a constant ²¹⁰Pb depositional flux and a constant carbon accumulation rate. Also shown in model A are the calendar dates 1994 and 1964, corresponding to the time of snow fence construction and the peak time of bomb-pulse radionuclide fallout, respectively. Model B shows the effect of doubling the ²¹⁰Pb flux after 1994 while maintaining a constant carbon accumulation rate; the slope of the line remains constant as in model A, but the calendar dates 1994 and 1964 occur at higher values of ln(A₀/A_Z). In model C, the ²¹⁰Pb flux is held constant while the carbon accumulation rate is doubled after 1994. In this case there is a doubling of slope after 1994 while the calendar dates remain fixed at the same values of ln(A₀/A_Z) as found in model A. In model D, both the ²¹⁰Pb flux and carbon accumulation rate are doubled after 1994, and the effects of these changes are a doubling of slope after 1994 as well as an increase in the values of ln(A₀/A_Z) corresponding to the calendar dates 1994 and 1964. Finally, in model E, both the ²¹⁰Pb flux and carbon accumulation rate tripled after 1994, resulting in a further increase in slope and further increases in the values of ln(A₀/A_Z) corresponding to the calendar dates 1994 and 1964, in comparison with model C.

Figure 5. Hypothetical modifications of CRS models for simulated core profiles showing cumulative carbon (g cm⁻²) vs. ln(A₀/A_Z). Carbon accumulation rates (CAR) and ²¹⁰Pb depositional fluxes were plotted with different arbitrary Y-intercept values to allow different models to be easily distinguished graphically, using equal values of CAR prior to 1994 in all models. Model A is a standard CRS model with ²¹⁰Pb flux and CAR held constant. In model B, ²¹⁰Pb flux is doubled after 1994 while CAR remains constant. In model C ²¹⁰Pb flux remains constant and CAR is doubled after 1994. In model D both ²¹⁰Pb flux and CAR are doubled after 1994. In model E both ²¹⁰Pb flux and CAR are tripled after 1994

For our measured core profiles, we estimated carbon accumulation rates from the slopes of linear regressions and considered the ²¹⁰Pb depositional flux as a fitting parameter that would allow the observed changes in slope to begin at a calendar date of 1995. This resulted in a fairly unique solution in terms of both carbon accumulation rate and the apparent change in ²¹⁰Pb flux. In order to account for the expected near-surface downward mobility of ²¹⁰Pb (Landis, Renshaw, and Kaste 2016), we excluded from our linear regressions any core segment in which ⁷Be was detected. In cores where ⁷Be was not measured, we excluded the top two core segments (top 2 cm) from the CRS plot regression analyses on the assumption that these were most likely affected by downward ²¹⁰Pb mobility. The chronologies implied by the resulting modified CRS models are considered accurate only below the depth limit of significant downward ²¹⁰Pb mobility (Landis, Renshaw, and Kaste 2016).

Individual modified CRS plots prepared for each of the measured core profiles are shown in Figure 6. Three “control” cores (CTL prefix in sample ID) were collected. These sites were designated in the original experimental design to represent ambient snow conditions; that is, outside the leeward snowdrift influence of the snow fence. Two of these cores (CTL-2 and CTL-4) showed good linear correlations over the entire depth range in which unsupported ²¹⁰Pb was present (Figure 6), consistent with a constant ²¹⁰Pb depositional flux and constant net organic carbon accumulation rates of 1.5 ± 0.1 to 1.7 ± 0.1 mg cm⁻² yr⁻¹ over the past ~100 years (Table 1). In these two cores, ²⁴¹Am was detected in the core segment that would be expected to contain bomb-pulse radionuclides based on its estimated calendar date, supporting the validity of the chronologic record estimated from the bulk values of each core segment. However, the ¹³⁷Cs peaks in these cores were found in segments well above the depth of ²⁴¹Am, indicating likely mobility of ¹³⁷Cs in the low-pH pore waters of the active layer. In core CTL-4, ²⁴¹Am was also detected in segments having calendar dates after the bomb-pulse peak. The third control core CTL-5 was, contrary to expectation, apparently influenced by the snow fence because its modified CRS plot showed a substantial increase in slope after 1994 (Figure 6). This core profile indicated that the carbon accumulation rate increased substantially after the construction of the snow fence (Table 1). The location of CTL-5 was ~65 m NNW of the western edge of the snow fence. Wind directions during winter at Toolik are generally southerly (katabatic winds), with occasional southeasterly or westerly trends (S. Sturm and Stuefer 2013). Occasional southeasterly winds during winter may thus have caused consistent snowdrift pile-up at CTL-5 following the snow fence construction. The unsupported ²¹⁰Pb inventory in CTL-5 is higher than that in both CTL-1 and CTL-2 (Table 1). This further supports our interpretation of CTL-5 as being influenced by enhanced snow accumulation that was caused by the snow fence, because additional snow accumulation would be likely to cause an increase in the ²¹⁰Pb depositional flux.

Figure 6. Modified CRS plots for measured core profiles showing cumulative organic carbon (g cm⁻²) vs. ln(A₀/A_Z). Note differences in ranges of x- and y-axis values. Dashed lines are best-fit linear regressions. Numerals adjacent to data points correspond to calendar years calculated from ln(A₀/A_Z). Net carbon accumulation rates (mg cm⁻² yr⁻¹) were obtained from slopes of regression lines. Data points having ⁷Be (detected or assumed) are excluded from regressions

Most of the cores affected by enhanced snow accumulation leeward of the snow fence (i.e., CTL-5, LOW-2, INT-1, and INT-5) had modified CRS plots best described by two-slope models having a distinct change in slope at depth, with concomitant increases in net carbon accumulation rate and ²¹⁰Pb depositional flux (Figure 6). Two of the cores (INT-3 and DEEP-4) were difficult to fit with modified CRS models because they had insufficient numbers of core segments below the depth of measured ⁷Be activity with which to define values for two different slopes. For the cores having successful two-slope modified CRS models, the carbon accumulation rates were given by the slopes of the linear regressions, and the ²¹⁰Pb fluxes were treated as a fit parameter that was adjusted to obtain a calendar date of 1995 for the change in slope. All of the two-slope models required increases in ²¹⁰Pb flux (ranging from +82 to +310 percent) to fit a calendar date of 1995 for the change in slope. This implies that the increase in winter snow accumulation leeward of the snow fence was accompanied by a roughly proportional increase in ²¹⁰Pb deposition or possibly that there had been an overall increase in ²¹⁰Pb flux deposited by summer precipitation after 1994. A gradual increase in summer ²¹⁰Pb flux in the tundra region may be likely as the thaw depth increases to allow more release of ²²²Rn from the active layer.

The portions of the successful two-slope modified CRS models corresponding to the period before the change in slope indicated net carbon accumulation rates ranging from 0.7 ± 0.1 to 1.8 ± 0.3 mg cm⁻² yr⁻¹, which bracket those obtained from the CRS models for the two control cores CTL-2 (1.5 ± 0.1 mg cm⁻² yr⁻¹) and CTL-4 (1.7 ± 0.1 mg cm⁻² yr⁻¹). The slopes of the two-slope modified CRS models beginning in 1995 (i.e., after snow fence construction) yielded uniformly larger carbon accumulation rates ranging from 2.8 ± 0.2 to 4.6 ± 0.3 mg cm⁻² yr⁻¹ (Table 1), roughly two to three times as much as those of the control cores.

Inventories of unsupported ²¹⁰Pb, ¹³⁷Cs, and ⁷Be

Total inventories and annual fluxes of unsupported ²¹⁰Pb, ¹³⁷Cs, and ⁷Be in the cores were calculated as described in the Methods section and are listed in Table 1. The total inventories of unsupported ²¹⁰Pb range from 79 ± 2 to 481 ± 5 mBq cm⁻². The values for control cores CTL-2 and CTL-4 (105 ± 3 and 79 ± 2 mBq cm⁻², respectively) are in good agreement with unsupported ²¹⁰Pb inventories obtained previously for sediments from Toolik Lake (75 mBq cm⁻²; Cornwell 1985) and for intertussock and ridgetop tundra soils (95 and 104 mBq cm⁻², respectively) from the Imnavait Creek watershed, ~16 km east of Toolik Lake (Cooper et al. 1991). Mean annual ²¹⁰Pb depositional fluxes obtained from our Toolik area cores ranged from 2.5 to 15.0 mBq cm⁻² yr⁻¹. Inventories and fluxes of ²¹⁰Pb are lowest in control cores CTL-2 and CTL-4 (3.3 ± 0.1 and 2.5 ± 0.1 mBq cm⁻² yr⁻¹, respectively) and moderately higher in most cores affected by enhanced snow accumulation (3.8 ± 0.1 to 4.5 ± 0.1 mBq cm⁻² yr⁻¹), with an anomalously high value measured in core DEEP-4 (15.0 ± 0.2 mBq cm⁻² yr⁻¹), which is further discussed below. The general trend of these values is consistent with that expected from increased snow depth, which appears to be accompanied by a higher depositional flux of atmospheric ²¹⁰Pb. However, because snow fence construction occurred in 1994, only one ²¹⁰Pb half-life (approximately twenty-two years) had elapsed between snow fence construction in 1994 and our core sampling in 2016, so the observed total inventories and fluxes mainly reflect prevailing conditions prior to 1994. Using the modified CRS models that we generated for this study, we verified that the best-fit post-1994 increases in ²¹⁰Pb depositional flux reproduced the observed increases in unsupported ²¹⁰Pb inventories relative to those of the control cores. Calculated pre-1994 ²¹⁰Pb depositional fluxes obtained from our best-fit modified CRS model parameters ranged from 2.1 to 2.9 mBq cm⁻² a⁻¹, with a mean value of 2.6 ± 0.4, which is in reasonable agreement with the ²¹⁰Pb depositional fluxes obtained from our measurements of unsupported ²¹⁰Pb inventories for the control cores (2.5 ± 0.1 and 3.3 ± 0.1 mBq cm⁻² a⁻¹).

The ²¹⁰Pb depositional fluxes that we derived from our measurements (Table 1) are in good agreement with a range of values from other high-latitude studies. A study of fifty-five sediment cores from five small arctic lakes near our Toolik field site yielded fluxes ranging from 2.9 to 7.7 mBq cm⁻² yr⁻¹ (Fitzgerald et al. 2005). A global compilation of ²¹⁰Pb depositional fluxes showed that the mean value at nine sites in 60° to 70° N latitude was 3.7 ± 3.5 mBq cm⁻² yr⁻¹ (Baskaran 2011). The ²¹⁰Pb flux measured in Greenland ice cores Summit and Dye 3 ranged from 0.4 to 2.5 mBq cm⁻² yr⁻¹ (Dibb 1992).

The inventories of ¹³⁷Cs in our cores ranged from 45 ± 1 to 104 ± 1 mBq cm⁻². Values obtained for the control cores CTL-2 and CTL-4 (83 ± 1 and 57 ± 1 mBq cm⁻², respectively) were within the range of those affected by the snow fence (45 ± 1 to 104 ± 1 mBq cm⁻²). These values agree well with ¹³⁷Cs inventories measured in tundra samples collected in 1990 from the Imnavait Creek watershed (Cooper et al. 1995) that had a range of 42 to 97 mBq cm⁻² (reported range of 77 to 177 mBq cm⁻², corrected for twenty-six years of decay from 1990 to 2016). Cooper et al. (1995) concluded that there was no evidence for widespread downslope ¹³⁷Cs transport in the Imnavait Creek watershed, although they noted that the elevated ¹³⁷Cs inventories in some of the Toolik Lake sediment cores near lake inlets indicated some apparent mobility. In several of our cores (INT-1, INT-5, CTL-5) ¹³⁷Cs peaks were concordant with ²¹⁰Pb chronologies, but most of the other cores showed apparent evidence for ¹³⁷Cs mobility generally toward the surface. This may be attributable to the low pH of pore waters (<5) along with possible advective transport of dissolved Cs through pore water flow or surface runoff during the spring snowmelt season and possible uptake by plant roots. Further study of ¹³⁷Cs activities in tundra pore waters and streams around Toolik Lake may provide useful constraints on the possible mechanism(s) and extent of ¹³⁷Cs mobility in this system.

The inventories of ⁷Be in the cores affected by the snow fence ranged from 0.8 ± 0.1 to 16.2 ± 0.2 mBq cm⁻², in remarkably good agreement with the range of 0.6 to 15.5 mBq cm⁻² measured in tundra samples of the Imnavait Creek watershed (Cooper et al. 1991). The range of ⁷Be inventories in snowpack cores collected from that watershed in 1989 and 1990 was 1.0 to 8.1 mBq cm⁻², similar to that measured in tundra samples. These inventories represent a subannual time period because of the relatively short (53 d) half-life of ⁷Be.

The anomalously high ⁷Be inventory of 16.2 ± 0.2 mBq cm⁻² in core DEEP-4 is correlated with its anomalously high unsupported ²¹⁰Pb inventory of 481 ± 5 mBq cm⁻² (Table 1). This may indicate that short-term lateral advective water flow carrying particles having adsorbed ⁷Be and ²¹⁰Pb into the local topographic lows (i.e., incipient thermokarst) that developed in the area of deepest snow accumulation. An additional factor that could help explain these anomalous values in core DEEP-4 is the observation of windblown leaves and surficial organic debris deposited leeward within 25 m of the snow fence in the deep snow zone (Figure 2; M. D. Walker et al. 1999).

Summary and conclusions

The principal objective of this study was to evaluate the effect of deeper winter snow accumulation on the net accumulation rate of soil organic carbon in the active layer of moist acidic tussock tundra, following the construction of a snow fence in 1994 as part of the ITEX (Jones et al. 1998; Welker, Fahnestock, and Jones 2000; Welker et al. 2005; Jespersen et al. 2018). We explored the application of ²¹⁰Pb as a chronometer to constrain decadal-scale carbon accumulation rates in the active layer. Based on the values we obtained from our modified CRS models for net soil organic carbon accumulation rates in the control cores and the pre-1994 portions of the other cores, we estimated the average background (i.e., ambient snow conditions) net soil organic carbon accumulation rate in the active layer of the study area at ~14 g C m⁻² yr⁻¹. Increases in winter snow depth over a period of twenty-two years (since construction of the snow fence in 1994) resulted in apparent increases of this rate at this site by ~90 to 310 percent to rates of ~26 to 44 g C m⁻² yr⁻¹. This represents a substantial increase in the soil carbon sink strength, as predicted by a number of climate warming models (e.g., Charman et al. 2013; Jiang et al. 2016; Chaudhary, Miller, and Smith 2017). Similar values of carbon accumulation rates were found in a chronologic study of soil profile TFS2 sampled ~50 m west of our study area (at a location where the dominant vegetation consisted of Sphagnum capillifolium, Aulacomnium turgidum, and Salix reticulata), which showed a substantial increase in carbon accumulation rate from ~14 g C m⁻² yr⁻¹ to ~43 g C m⁻² yr⁻¹ corresponding to an apparent ecosystem change from minerotrophic to oligotrophic conditions after the Little Ice Age at ~1850 CE (Taylor et al. 2019).

These results suggest that though deeper snow leads to elevated winter CO₂ emissions, the positive effects of deeper snow on shrub growth, aboveground biomass and net CO₂ fixation, and presumably greater root biomass collectively overcome the gaseous losses of C, leading to the positive net effect of increasing active layer soil carbon accumulation rates. This positive and negative feedback contest between ecosystem C losses (gaseous emissions) and gains (combinations of C fixation, biomass accumulation, and soil C formation) has been at the center of some of the earliest studies of tundra responses to changes in conditions with ITEX (Welker et al. 1997). The findings in this study suggest that over decades of deeper snow in winter, soil carbon accumulations of tussock tundra will be appreciably greater due in large part to vegetation changes associated with shrub growth and elevated soil carbon inputs (M. Sturm et al. 2005; Tape et al. 2012; Lynch et al. 2018) and lower organic matter decomposition rates as soil water content increases with snow (Blanc-Betes et al. 2016), even though CO₂ emissions in winter are higher.

However, abiotic and biotic interactions may change the magnitude, extent, redistribution, and even longer-term net effect of this increased soil carbon sink. For instance, thawing of previously frozen soils may initially stimulate rates of microbial degradation of organic carbon, and as accessible carbon pools are depleted over time and soil conditions change, microbial communities may shift toward more effective decomposition and metabolic strategies (Ricketts et al. 2016; Melillo et al. 2017). Competing rates of soil methanogenesis versus respiration, which respectively emit CH₄ and CO₂ to the atmosphere, depend in large part on whether the soils are wet or dry (Blanc-Betes et al. 2016). These changes in soil carbon accumulation will not occur in isolation from other climate change–induced effects. Permafrost thaw may change the ecohydrology of the system (Blanc-Betes et al. 2016; Throckmorton et al. 2016; Jespersen et al. 2018), and an even greater dominance of the landscape by shrubs with more leaf area will continue to increase carbon fixation (Leffler et al. 2016; Lynch et al. 2018), leading to greater carbon allocation to roots. Denser shrub canopies may cause shading and soil cooling, countering greater soil carbon emissions (Mekonnen, Riley, and Grant 2018) and amplifying the biocomplexity of the tussock tundra’s response to climate change (D. A. Walker, Epstein, and Welker 2008; Flower et al. 2019).

The multidecadal net effect of the complex ecohydrologic interactions in the arctic tundra are difficult to predict (Welker et al. 1997). The gradual expansion of shrubs along with the effects of warmer winters on the duration and magnitude of winter snowpack may add to the longer-term complexity of carbon source–sink relations in this arctic ecosystem (Fahnestock et al. 1999; Natali, Watts, and Rogers et al. 2019). Nonetheless, it is likely that the rapid climate change now being observed in the Arctic will continue, and the consequences of changing winter snow depth accompanying climate change will likely be one of the most important factors in the future evolution of the arctic ecosystem (Welker, Fahnestock, and Jones 2000; Natali et al. 2019). The approach demonstrated here may provide a new tool for better understanding of the response of the arctic tundra’s soil organic carbon inventory to the changing climate.

Acknowledgments

We thank the staff of the Toolik Field Station and CH2MHill Polar Services for logistical support, and we are grateful to the anonymous reviewers for their insightful and constructive comments on an earlier version of this article.

Disclosure statement

No potential conflict of interest was reported by the authors.

Supplemental Material

Supplemental data for the article can be accessed on the publisher’s website

Additional information

Funding

This research was funded by the U.S. Department of Energy, Terrestrial Ecosystem Science Program (DE-SC0006607) to MAG-M and by the U.S. National Science Foundation, Office of Polar Programs over the 20+-year duration of the experiment with awards from the Arctic System Science and Arctic Observing Network programs, through award numbers 1504141, 1433063, 0119279, 0856728, 0632184, 9617643, 9321730 to JMW. KCD acknowledges receipt of a Grad Student Research Grant from the Geological Society of America.