Implications of alder shrub growth for alpine tundra soil properties in Interior Alaska

ABSTRACT

The increase in deciduous shrub growth in response to climate change throughout the Arctic tundra has uncertain implications, in part due to a lack of field observations. Here we investigate how increasing alder shrub growth in alpine tundra in Interior Alaska corresponds to active layer thickness and soil physical properties. We documented increased alder growth by combining biomass harvests and dendrochronology with the analysis of remotely sensed Normalized Difference Vegetation Index and fire history. Active layer thickness was measured with a tile probe and carbon and nitrogen pools were assessed via elemental analysis. Shallower organic layers under increasing alder growth indicate that nitrogen-rich, deciduous litter inputs may play a role in accelerating decomposition. Despite the observed reduction in organic carbon stocks, active layer thickness was the same under alder and adjacent graminoid tundra, implying deeper thaw of the underlying mineral soil. This study provides further evidence that the widely observed expansion of deciduous shrubs into graminoid tundra will reduce ecosystem carbon stocks and intensify soil–atmosphere thermal coupling.

PLA IN LANGUAGE SUMMARY

Two consequences of rapid climate warming in the Arctic, where grass-like plants dominate under very cold conditions, are an increased growth and occurrence of shrubs and associated thaw of frozen ground. This exposes organic matter in soils to microbes that can decompose it into carbonaceous greenhouse gases, but some of this carbon loss may be offset by the increased plant growth. Here, we investigate the impacts of greater shrub presence on soil properties at five sites in Alaska. We documented shrub growth by analyzing satellite images, which can help us understand the productivity and/or leaf coverage at each site back in time, and annual growth rings in shrub stems, which show how old the shrubs are and how much they grow each year. We also measured the depth of soil thaw in the field and its organic matter content in a laboratory. Where shrubs were more common, we found a thinner layer of organic matter at the soil surface. Thaw depth remained the same, which may indicate that the presence of shrubs results in deeper thaw of the mineral soil. Our findings support the hypothesis that shrub expansion will further enhance warming-driven increases of greenhouse gas emissions from Arctic landscapes.

KEY POINTS

• Trends in dendrochronology and Normalized Difference Vegetation Index reveal increasing growth of alder shrubs in Interior Alaska.

• More alder cover results in the loss of the soil organic layer and thus soil C and N that is not offset by more shrub biomass.

• Increasing alder growth may promote permafrost thaw not captured by tile probe active layer thickness monitoring.

KEYWORDS:

• Soil organic carbon
• NDVI
• active layer thickness
• shrubification
• permafrost

Introduction

Deciduous shrubs have become more abundant across the Arctic (Sturm et al. 2001; Tape, Sturm, and Racine 2006; Elmendorf et al. 2012; Frost and Epstein 2014; Shevtsova et al. 2020), manifesting in greater shrub height, infilling of open shrub patches, and shrubline advance (Myers-Smith et al. 2011). The causes of this shrub expansion (“shrubification”) are complex (Martin et al. 2017) but include warming (Tape et al. 2012), increased soil moisture (Jespersen et al. 2018; Chen, Hu, and Lara 2021), greater availability of nutrients (Vankoughnett and Grogan 2016), and changes in disturbance frequency (Racine et al. 2004; Frost et al. 2013; Pearson et al. 2013; Plante et al. 2014; Andruko, Danby, and Grogan 2020; Huebner, Buchwal, and Bret-Harte 2022; Liu et al. 2022). Though shrubification has been widely documented, the ramifications of increased shrub growth on ecosystem processes remain less clear (Chapin et al. 2005; Mekonnen et al. 2021). Specifically, vegetation cover and type are known to impact net radiation and surface energy exchange (Jorgenson et al. 2010; Loranty et al. 2018), raising the possibility that large-scale shrubification may trigger important changes in permafrost distribution across the Arctic and the large amounts of organic carbon (C) contained within it (Miner et al. 2022).

Shrubification is generally expected to promote deepening of the active layer—the layer of soil above permafrost which seasonally thaws—ultimately resulting in soil carbon loss due to greater microbial activity, especially during the cold season (Sturm et al. 2005; Kropp et al. 2020; Way and Lapalme 2021). Tall shrubs trap snow in winter, which insulates soils from subfreezing air temperatures, and also protrude from the snow surface during spring, which lowers surface albedo and accelerates snow melt (Lawrence and Swenson 2011; Vankoughnett and Grogan 2016; Jafarov et al. 2018; Wilcox et al. 2019). However, cold season shrub–snow interactions may change with climate warming if more frequent snowmelt events occur in autumn, because the resulting ice layers can decrease winter snow drifting and facilitate soil cooling (Barrere et al. 2018). Growing season processes further complicate the picture, because greater evapotranspiration and shading from tall shrubs may lead to cooler and drier soils (Nauta et al. 2015; Vankoughnett and Grogan 2016). Shrub-induced soil drying may even lead to cooler soil temperatures year-round (Göckede et al. 2017); however, the long-term implications of shrubification and surface hydrology on permafrost remain uncertain (Loranty et al. 2018; Heijmans et al. 2022). Manipulation experiments have found evidence of increased soil temperatures caused by shrub presence (Aguirre, Benhumea, and McLaren 2021) but also by shrub removal (Blok et al. 2010), resulting in opposite effects on nutrient cycling and active layer thickness (ALT). Thus, it remains unclear what role increased shrub growth in situ plays in modulating soil physical properties.

Adding to the complexity of assessing the consequences of Arctic shrubification are species effects, which cannot be evaluated at the current resolution of satellite imagery. Shrubification primarily involves three different shrubs, alder (Alnus spp.), birch (Betula spp.), and willow (Salix spp.; Tape, Sturm, and Racine 2006). These genera differ in their tissue chemistry (i.e., nitrogen [N] and tannin content; Christie et al. 2015) as well as their C allocation strategy (Vankoughnett and Grogan 2016; Lynch et al. 2018; Street et al. 2018), both of which impact shrubs’ competitive abilities, herbivore interactions (Olofsson et al. 2009; Cahoon et al. 2012; Christie et al. 2015), and effects on soil properties. Alder, which has N₂-fixing symbionts, is expanding across parts of Alaska (Tape, Sturm, and Racine 2006). Alder presence increases accumulation and bioavailability of N in soil (Hart, Binkley, and Perry 1997; Rhoades et al. 2001), and a recent study showed that alder stands mediate local diversity by encouraging growth of birch and evergreen shrubs and limiting understory tundra species (Wallace and Baltzer 2020). However, additional field studies are needed to explore the relative effects of alder on permafrost soil properties.

Despite the ubiquity of Arctic shrubification and uncertainty surrounding its ecosystem consequences, few field studies have quantified the implications of increased shrub growth on ecosystem and soil C sequestration (Mekonnen et al. 2021). Some have proposed that shrubification has the potential to offset part of the C emissions associated with permafrost thaw. For instance, one modeling study predicted that a shift from graminoid and nonvascular vegetation to woody vegetation in North America under the RCP8.5 scenario would result in additional C sequestration during the twenty-first century because the greater productivity of shrubs, lower decomposition rates of shrub litter, and shading of the soil surface during the growing season compensated for warming-driven increases in heterotrophic respiration (Mekonnen et al. 2018). For temperate soils, chronosequence studies found that woody encroachment into grasslands resulted in net ecosystem C gains in drier regions but C losses in wetter regions (Jackson et al. 2002). An analogous assessment of soil C stocks in the circumpolar permafrost region is particularly challenging because differences in microtopography (and associated differences in net radiation, energy partitioning, and surface hydrology) and buried horizons (formed under previous vegetation types and with legacy C) often inhibit direct comparison of nearby soils under graminoid and shrub tundra.

In recent decades, remote sensing has greatly advanced our understanding of large-scale vegetation change patterns (Myers-Smith et al. 2020). However, it remains limited in its efficacy in detecting soil properties (Odebiri, Odindi, and Mutanga 2021). Therefore, combining high-resolution remote sensing data sets with ground observations is critical to understanding the implications of vegetation change in Arctic landscapes. Here we combine surveys of alder shrub biomass and soil properties with dendrochronology and satellite observations of plant productivity and fire history at five subarctic alpine graminoid tundra field sites in Interior Alaska. This hybrid experimental approach, which combines retrospective remote sensing observations with in situ field observations, offers opportunities to overcome the limitations imposed by pure remote sensing observations. We sought to target belowground and aboveground measurements at sites that were actively responding to drivers associated with vegetation change in the tundra. We asked: What are the implications of increasing alder biomass and/or cover for surface soil properties and active layer thickness?

Materials and Methods

Site description

Shrub biomass data and soil samples were collected at five alpine tundra sites in Interior Alaska along the James W. Dalton Highway (Alaska Route 11), between milepost (MP) 95 and 240, in August 2019 (Table 1), where vegetation is a mosaic of boreal forest, wetlands, and (at higher elevation) alpine tundra.

Table 1. Overview of study sites in Interior Alaska, average (standard deviation), n = 21/site.

	1	2	3	4^a	5
Latitude (°N)	66.34032	66.47886	66.38084	66.59785	68.01826
Longitude (°W)	−150.45989	−150.7016	−150.51224	−150.71356	−149.71002
Elevation (m)	592 (6.2)	605 (2.8)	440 (2.2)	499 (4.0)	789 (20.6)
Aspect	E/NE	E/NE	E/NE	W/NW	W
Parent material	Granite	Granite	Granite	Phyllite	Phyllite
Tundra community^b	Open low mesic shrub/birch–ericaceous shrub	Vaccinium dwarf shrub	Vaccinium dwarf shrub	Open tall alder shrub	Vaccinium dwarf shrub/open tall alder shrub
Mean active layer thickness (cm)	66.0 (25.8)	66.1 (10.7)	67.3 (20.6)	56.2 (21.4)	64.3 (16.8)
Mean Aldermax cover (%)	18 (14)	4 (4)	5 (4)	25 (17)	15 (16)
Mean alder AGB (kg m^−2)	1.1 (2.1)	0.3 (0.7)	0.7 (1.1)	3.8 (6.1)	3.2 (6.9)

^aSite 4 was underlain by shallow bedrock at 25.2 ± 6.2 cm depth (52 percent of survey points). ALT was measured at 48 percent of survey points.

^bCommunity composition as in Viereck (1992).

Sites were selected based on geomorphology, vegetation type (alpine tundra; i.e., absence of trees) and abundance of alder shrubs while also trying to minimize differences (a) in topography between areas with more or less alder at each site and (b) in tundra community composition, aspect, and parent material between sites. Sites 1 to 3 were located within about 20 km of Finger Mountain (MP 98), where granite tors represent remnants of bedrock hills that formed from a Cretaceous pluton and remained unglaciated during the last ice age (D. A. Walker et al. 2009). Site 4 was located near Connection Rock (MP 121) and site 5 in the southern foothills of the Brooks Range (MP 237).

At each site, we established three parallel 30-m transects incorporating areas of higher and lower shrub abundance (Figure 1). Transects were separated by 10 m. Every 5 m along each transect, we recorded the Global Positioning System location and elevation (etrex 10, Garmin) and thaw depth (cm) using a tile probe (Mackay 1977) and quantified alder biomass using a sampling plot (see section “Vegetation surveys and alder biomass and growth history”). Because measurements were conducted at the end of the growing season (22–28 August 2019), we assume that our thaw depth measurements represent the ALT at each site in 2019. Thickening active layers can be indicative of vertical degradation of permafrost following some disturbance (e.g., climate change, fire, vegetation change), depending on how the permafrost formed (Shur and Jorgenson 2007). We interpreted ALT as a permafrost stability proxy in the field, assuming the presence/absence of alder in the gradient transect design the main predictor variable. In addition, alder cover was estimated across 10-m intervals of each transect, and alder disks and soil samples (O horizons) were collected to further assess vegetation history and soil properties, respectively (see sections “Vegetation surveys and alder biomass and growth history” and “Soil collection and analysis”).

Figure 1. Map and photographs of alpine tundra with varying degree of alder abundance. “C” and “D” denote locations where soil cores and alder disks were sampled, respectively. Horizontal gray bars with arrow caps indicate 10-m intervals across which vegetation cover was estimated. All sites with alder >2 m tall also contained shorter alders. Sites 1 and 3 were affected by the Dall City Fire (2004) and site 4 was affected by the North Bonanza Fire (2005). ABoVE standard Alber’s equal area projection. Source: U.S. Census Bureau TIGER geographic database (U.S. Census Bureau 2022).

Vegetation surveys and alder biomass and growth history

Vegetation cover and composition were determined for each transect along 10-m intervals (three per transect, nine per site). For each transect, these data are reported at the 0-, 10-, 20-, and 30-m points, whereby the two cover estimates for the 10- and 20-m points were averaged. Using three levels (upper canopy, >2 m; lower canopy, ≤2 m; and ground) and a transect width of 0.5 m, we quantified the relative cover of alder (Alnus viridis (Chaix) DC.) and other tall shrub or tree species and of ground cover, separated into dwarf shrub, graminoid, moss, lichen, and bare soil. Additionally, we estimated the maximum alder coverage (Alder_max) at each location as the greater proportional alder coverage of the upper and lower canopy. Where two surveys overlapped, the Alder_max value was averaged.

Alder biomass was quantified by recording the number of alder stems within a 0.5 m × 0.5 m plot location every 5 m along each transect (twenty-one per site) and measuring the basal diameter of each stem using calipers. We also harvested five two six representative stems (including foliage) at sites 2 to 4 and measured their basal diameter. These samples were oven-dried to constant mass and weighed.

We used an allometric fit for the harvested alder stems that was used to predict alder biomass for all plot locations (Equation (1)(1) $\mathit{AG}{B}_{\mathit{model}}=\mathit{aB}{D}^b,$(1) ). Previous studies have found that the basal diameter of tall-stature shrubs is an accurate predictor of aboveground biomass (Berner et al. 2015), but few studies have analyzed this relationship for shrubs growing on permafrost. The fit was calculated as a nonlinear least squares regression using the stats package in R (R Core Team 2021) as in Berner et al. (2015):

(1) $\mathit{AG}{B}_{\mathit{model}}=\mathit{aB}{D}^b,$(1)

where AGB_model is the estimated aboveground dry alder biomass (including leaves, cones, and woody material) per stem, BD is the measured basal diameter of each stem, and a and b are coefficients of the regression representing linear and nonlinear factors, respectively. To correct for nonrandom residuals (such as overprediction of small BD shrubs) we weighted the regression by AGB_model⁻¹. We determined model error predictions from a Monte Carlo simulation of confidence intervals (R predictNLS v1.0–4), which constructs more realistic estimates of nonlinear models than Taylor expansion. We found that the best-fit allometric equation to relate BD to AGB_model had coefficients of a = 10.24 and b = 2.994 (R² = 0.95, n = 16, 1σ = 4.6 g stem⁻¹), which falls within one standard error of previously published coefficient estimates from boreal Alaska for the same species (Berner et al. 2015).

To establish the growth history of alder at our sites, basal stem disks were collected from multiple individuals in each transect (n = 40 total, Figure 1). We collected a disk of the alder with the largest basal diameter (five to eleven per site) within the transects at all sites and an additional two to six disks of large alders outside the transects at sites 2 to 5. Although study designs which sample selectively based on growth can impart biases to dendrochronological data sets (Nehrbass-Ahles et al. 2014), our approach provides an estimate of the timing of alder establishment (maximum age) and growth history at each site.

Disks were air-dried, sanded to 1,500 grit, and measured with a sliding bench micrometer. Shrubs ranged in age from 2 to 81 years, with a mean age of 21.8 years. Ring-width series were examined for dating errors in COFECHA (Holmes 1983). Age-related growth effects were removed from the shrub ring chronologies by detrending with age-based regional curve standardization. Signal-free chronologies were then produced using repeated divisions of the raw ring widths by the chronology indices (ten iterations). Regional curve standardization detrending and signal-free chronology development were performed in CRUST (Melvin and Briffa 2014) both at the site level and using a composite data set of all sampled shrubs.

Trends in NDVI

To evaluate productivity and fire history at each site, images from Collection 1 (Tiers 1 and 2), Landsat 5, 7, and 8 were extracted for Julian days 152 to 243 of years 2000 to 2020 using the Google Earth Engine Python API (van Rossum and Drake 2009; Gorelick et al. 2017; Berner et al. 2020) after field sampling. The study period 2000 to 2020 was chosen due to the limited number of Landsat observations in Alaska before 2000. Images were extracted using the central geographic location for each site with a 50-m circular buffer. These included images with WRS paths 71 to 74 and rows 12 to 14. Images were filtered for clouds, snow, water, band saturation, and geometric uncertainty following Berner et al. (2020). Additionally, land cover classification was analyzed within a 50-m circular buffer around the central geographical coordinate of each site using the ABoVE Landsat-derived Annual Land Cover Across ABoVE Core Domain, 1984–2014 (Wang et al. 2019).

Landsat Normalized Difference Vegetation Index (NDVI) was calculated for each image using the red and infrared bands from each sensor (bands 3 and 4 for Landsat 5 and 7; bands 4 and 5 for Landsat 8). NDVI_max was taken as the maximum observation per year per site. We note that NDVI_max is sensitive to the number of quality images in a summer, because the likelihood that true peak greenness was measured increases as the number of observations increases. The difference between phenologically corrected NDVI_max and Landsat-measured NDVI_max has been shown to be insignificant when at least seven observations are available for a particular site and year (Berner et al. 2020). At least seven observations were available for all sites and years except for site 2 in 2012 (n = 6) and for site 4 in 2003 (n = 5) and 2006 (n = 4). No quality images were available in 2003 for sites 1 or 2.

Trends in NDVI_max over time were analyzed for significance using the Mann-Kendall trend test and Theil-Sen slope using the zyp package (Bronaugh & Consortium, 2019) in R (R Core Team 2023), which accounts for autocorrelation in temporal data using Yue et al.’s method (Yue, Pilon, and Cavadias 2002). Similar methods for detecting NDVI trends have been used for Arctic and boreal sites (Berner et al. 2013, 2020; Guay et al. 2014). Wildfires that may have affected study sites in the time period 2000 to 2020 were explored using the MODIS burned area product (MCD64A1, monthly, 500 m) with a 3 km × 3 km buffer (Giglio et al. 2018). Major fires were identified using the Alaska Wildland Fire Information Map Series from the U.S. Bureau of Land Management (Alaska Fire Service (AFS) Alaska Wildland Fire Information Map Series n.d.).

Soil collection and analysis

To assess the relationship among alder biomass, cover, and the properties of the soil organic layer, a total of twenty-eight soil cores (three to eight per site, Figure 1) were collected with an electric auger with a scalloped edge and an inner diameter of 2.54 cm (Nalder and Wein 1998) from the surface (dead moss or other O horizon) to a depth of 6 to 36 cm below the surface, determined by the distance to frozen or saturated ground, saprolite, or auger length. Thus, only complete organic layers are analyzed here. Sampling locations were chosen to capture end members of alder abundance and a representative range between. Cores were stored intact in capped plastic liners and frozen (−20°C) within 24 hours until analysis. Frozen soils were sectioned by horizon, oven-dried, and homogenized to <40 μm diameter or ground to powder, respectively. Aliquots (1–6 mg of organic and 8–25 mg of mineral soil) were weighed into tin capsules (5 mm × 9 mm) and analyzed for their total C and N content with an elemental analyzer. Percent C was used to determine whether soil horizons had been correctly classified as organic (≥20 percent C) or mineral (<20 percent). The total amount of C or N in each complete organic layer was calculated using bulk density, depth, and elemental data.

Statistical methods

To examine the influence of variable state factors and burn history on response variables (C and N stocks, organic layer thickness, ALT), Welch’s two-sample t tests were performed on soil measurements with <10 percent Alder_max coverage using the stats package in R (R Core Team 2023). Because burn sites (1, 3, and 4) burned around the same time (2005 and 2006), burn history was treated as binary (yes or no). t Tests failed to reject the null hypothesis that the mean response variable was not significantly different according to parent material, aspect, or burn history. We therefore proceed with analysis assuming that these variables do not significantly alter relationships among alder and response variables.

To examine the relationship between alder (biomass and coverage) and response variables, linear mixed effect (LME) models were created using the lme4 package in R (Bates et al. 2015; R Core Team 2023). All response variables but ALT were square root transformed prior to regression analysis to meet the assumption of normality (square root transformation was chosen to account for zero values). Independence of residuals was examined using the Durbin-Watson test for autocorrelated errors, normal distribution of residuals was examined using the Shapiro-Wilks test, and homoskedasticity was examined using the score test for nonconstant error variance using the car and stats packages in R (Fox and Weisberg 2019; R Core Team 2023). Models of the relationship between ALT and alder biomass or cover did not meet assumptions. LME models were chosen to account for the random variance created by site and transect in the study design. Models that treated both site and transect (as well as site and transect independently) as random effect intercepts were considered. The best fit of those models was chosen using Akaike’s information criterion (Bolker 2022); in all cases, models that treated only site as a random effect were chosen.

Results and Discussion

Alder growth history

Our results suggest that alder biomass and cover increased in recent decades at our alpine tundra sites, though with varying disturbance histories. Alder was consistently the dominant tall deciduous shrub species at all sites (aspen, birch, and willow accounted for <5 percent), and its growth chronologies suggest an increase in growth in the past four decades (Figure 2a). Among alder individuals sampled, which were the largest individuals per site, age varied by a factor of 3.7 across sites, ranging from a mean age of nine years at site 1 to thirty-four years at site 2; 43 percent of sampled stems established before 2000, and 53 percent established between 2000 and 2010. The cross-site growth chronology showed steadily increasing annual alder growth from the mid-1970s to about 2010, such that average radial growth post-2010 was approximately double the norms from the late 1970s. Site-specific growth chronologies displayed variable trends, with sites 2, 3, and 4 showing increasing recent growth and sites 1 and 5 displaying neutral or negative trends, but generally lacked sufficient within-site replication for reliable trend analysis (only one site generated a chronology with expressed population signal >0.85). Previous population studies of A. viridis report variable stand age structure (Lantz, Gergel, and Henry 2010), with one recent study reporting the mean age of a stand on Alaska’s North Slope to be 77 ± 15 years (Drew et al. 2023). Together, the low average age and increasing growth of large alder shrubs at these study sites point toward recent response to change, which may include recovery from fire disturbance and/or response to more favorable climatic conditions.

Figure 2. Records of (a) alder ring width index (1937–2018) and (b) maximum yearly NDVI (2000–2020) at five alpine tundra sites. Vertical bands indicate the following burn events and the sites that they affected: Dall City Fire of 2004 (sites 1 and 3) and North Bonanza Fire of 2005 (site 4), verified by MODIS MCD64A1. Black lines show intersite local regressions (span = 0.3), weighted by sample size in (a); gray band shows standard error.

NDVI_max was obtained from Landsat for the entirety of each study site as a metric for site productivity through time and disturbance history, because Landsat resolution is too coarse to capture the fine-scale vegetation gradient targeted for field sampling. Though there was a significantly positive trend in average NDVI_max for all sites over the period 2000 to 2020 (tau = 0.495, Theil-Sen slope = 0.0044, p = .003; Figure 2b), this trend was stronger among the burned sites (1, 3, and 4) than the unburned sites.

Time series of NDVI_max support an increase in shrub cover and/or productivity at all sites, though evidence suggests variable site history. NDVI is known to be linked to leaf area index in both shrubland and graminoid tundra ecosystems (Gamon et al. 1995; Filella et al. 2004; Riedel, Epstein, and Walker 2005). Though NDVI is also known to be linked to biomass (Boelman et al. 2003), this correlation may be weakened in shrublands due to the abundance of woody biomass (Gamon et al. 1995; Cunliffe et al. 2020); NDVI at other phenological points (e.g., NDVI pre leaf out) may be better for estimating total aboveground biomass (Boelman et al. 2011). The mean NDVI_max for all sites between 2000 and 2020 ranged from 0.6 ± 0.06 to 0.8 ± 0.04 (mean ± SD).This range, along with the positive trend observed over two decades, is consistent with shrubland tundra summer NDVI with high shrub cover (Beck and Goetz 2011), suggesting that an increase in shrub cover is the main driver of this trend, rather than an increase in productivity of graminoid species. However, the fire response evident at sites 1, 3, and 4 is consistent with those previously observed in boreal forest succession (Beck and Goetz 2011). Specifically, the NDVI_max recovered and then exceeded prefire levels, likely due to increased deciduous recruitment (Johnstone et al. 2010; Beck and Goetz 2011; Beck et al. 2011), though this postfire recruitment may be linked to prefire deciduous species density (Johnstone et al. 2004). Though site selection targeted alpine tundra using vegetation characteristics described for the area (Hewitt et al. 2017), it seems possible that especially those sites at which alder growth was not observed prior to burn (sites 1 and 4) are in successional stages. However, analysis of percentage landcover classification for each burned site suggests that these sites were not classified as forest prior to burn (Wang et al. 2019). In 1984 (earliest year of data availability), site 1 was classified as 100 percent herbaceous; site 3 was 63 percent tussock tundra, 25 percent low shrub, and 13 percent herbaceous; and site 4 was 83 percent low shrub, 33 percent tall shrub, and 17 percent tussock tundra. Nevertheless, though more extensive studies would be needed to precisely determine the drivers of shrub growth at these study sites, our observations and data support the assumption that alder biomass and cover are increasing adjacent to alpine tundra.

Alder biomass and cover impacts on soils

The site random effect accounted for very little of the variance explained by the total models, which is clear in the R² values (method described in Nakagawa and Schielzeth 2013) shown in Figure 3a–f (Bartoń 2023). Only in the case of organic thickness as a function of aboveground biomass did site explain any notable variance (0.0146 random effect and 0.8392 residual). Thus, the marginal fit (which considers fixed effects—i.e., the predictor variable—but not random effects; i.e., site) is shown in Figure 3a–f.

Figure 3. Relationships between modeled alder aboveground biomass (AGB_model) or maximum alder cover (Alder_max) in alpine tundra and square root transformed organic layer (a), (b) thickness, (c), (d) C stock, (e), (f) N stock, or (g), (h) ALT. Soil response variables (organic layer thickness, C stock, N stock) are square root transformed. The marginal fit taken from LME models is shown in (a) to (f). Marginal R² (R² m) and conditional R² (R² c) are shown in inset plots. Ninety-five percent confidence interval bands were estimated by nonparametric bootstrapping (Bates et al. 2015; Wickham 2016; Bartoń 2023; R Core Team 2023).

Organic layer thickness was negatively related to both AGB_model and Alder_max (Figure 3a,b; R² m = 0.25 and 0.43 for AGB_model and Alder_max, respectively; n = 28), leading to lower organic C stocks under alder (Figure 3c,d; R² m = 0.16 and 0.31, respectively; n = 28) and lower organic N stocks under alder (Figure 3e–f; R² m = 0.14 and 0.31, respectively; n = 28). Alder_max was a stronger predictor of organic layer thickness, C stocks, and N stocks than AGB_model.

Results indicate that organic layer properties may be responding to changes in alder litter quantity and/or quality over changes in alder biomass. Deciduous shrub leaf litter is relatively labile compared to other tundra plant functional types, which also promotes faster decomposition (Weintraub and Schimel 2005). Moreover, alder increases bioavailable N through a symbiotic relationship with N-fixing Frankia bacteria (Franche, Lindström, and Elmerich 2009), which may increase rates of decomposition and alter C balance (Salmon et al. 2019; Heslop et al. 2021; Pérez et al. 2021). Additionally, an increase in alder cover may affect organic layer C stocks through crowding of shade intolerant ground cover species, like mosses and lichen (Cornelissen et al. 2001; M. D. Walker et al. 2006; Pajunen, Oksanen, and Virtanen 2011; Elmendorf et al. 2012; Fraser et al. 2014; Mekonnen et al. 2021). The thick moss layer in the tundra promotes C accumulation due to slow decomposition rates (Hobbie et al. 2000) and reduces soil temperature fluctuations in the growing season (Turetsky et al. 2012; Fisher et al. 2016). A recent manipulative study found that shrub presence, rather than litter inputs, changed nutrient cycles and microbial processes indirectly though changes to soil temperature (Aguirre, Benhumea, and McLaren 2021). Also, the relatively low predictability of the linear models presented here may be attributed to heterogeneity in shrub canopy and soil properties in alder shrub tundra, which is not fully captured with the chosen sample design.

Contrary to expectations, no difference in ALT was observed with varying alder abundance (Figure 3g–h). Measured ALT was 63 ± 19 cm (mean ± SD) beneath alder and 66 ± 19 cm beneath tussock tundra. Previous studies have found both positive and negative relationships between ALT and shrub growth (Blok et al. 2010; Tape et al. 2012; Wilcox et al. 2019), because the expansion of shrubs on the tundra introduces complex feedback mechanisms which may produce net cooling (canopy shading, increased evapotranspiration) or net warming (increased albedo, increased snow accumulation, increased decomposition, increased fire occurrence) of the soil (Sturm et al. 2005; Lantz, Marsh, and Kokelj 2013; Loranty et al. 2018; Kropp et al. 2020). We found no evidence of an effect in either direction of shrub presence on ALT and thus permafrost stability. We therefore found that the observed organic horizon shallowing beneath alder had no effect on ALT. Previous work in the boreal region has shown that a greater (canopy and understory) leaf area index and organic layer thickness (particularly moss) are associated with a shallower ALT (O’Donnell et al. 2009; Fisher et al. 2016; Matyshak et al. 2017). It is perplexing then that the observed organic horizon shallowing beneath alder had no effect on ALT, especially because moss cover and thickness tended to be lower under alder. Taken together, the combination of a lack of trend in thaw depth and a negative trend in organic horizon thickness indicates that ALT could be deepening into mineral soil, albeit with high variability. However, this proposed type of active layer deepening is not typically observed with vegetation change (Fisher et al. 2016) and could still be influenced by fire history.

The role of wildfire should also be considered in the interpretation of both organic horizon and ALT results. Fire inhibits shrub growth immediately after occurrence but promotes shrub growth in the long term on the order of a few decades (Racine et al. 2004). Because shrubs can resprout following fire and have resilient mycorrhizal symbionts, they may facilitate vegetation recovery (Lantz, Gergel, and Henry 2010; Hewitt et al. 2013; Heim et al. 2021; Mekonnen et al. 2021). It is therefore possible that the observed organic horizon properties are predictive of alder presence rather than vice versa, because exposed mineral soil after fire provides area for seed germination (Racine et al. 2004). However, in this case we would expect to see additional evidence of active fire recovery at burned sites. We do not observe any significant difference in organic layer thickness, C or N stocks, or ALT in non-alder tundra between burned and unburned sites. We also do not see evidence of increased ALT corresponding with either shallow organic layer or burn occurrence. In a boreal forest study, thicker organic layers resulted in decreased ALTs, and fire occurrence resulted in increased ALTs through direct effects on canopy, soil moisture, and organic layer thickness (Fisher et al. 2016). Tundra vegetation usually recovers within five to ten years in small tundra fires, but compositional shifts may persist for decades after severe fires and/or result in regime shifts (Bret-Harte et al. 2013). It is evident that tundra fires have a role in promoting shrub growth and altering organic layer properties in the short term, but based on the results of this study, it seems evident that shrub growth facilitates the persistence of these effects. Increased burn area and frequency of large fires were observed in the boreal region in the twentieth century (Kasischke and Turetsky 2006; Kelly et al. 2013), and tundra fires are projected to increase in frequency and severity (Hu et al. 2010); thus, fire will have an increasingly present role in moderating tundra shrub growth.

Current Earth system models suggest that the increase in net primary production due to warming will offset the loss of C from Arctic soils (Todd-Brown et al., 2013), but our data suggest the opposite may be the case, considering that the organic layer is high in C content (Tarnocai et al. 2009). In the set of locations for which a soil core was harvested (n = 28), we found that the maximum stock of soil C in the organic horizon was 20 kg soil C m⁻² (associated with 1.7 kg alder C m⁻² at the same plot location) and the maximum stock of alder biomass C was 5.7 kg alder C m⁻² (associated with 2.6 kg soil C m⁻²). Thus, at the observed maximum alder biomass, there was a loss of about 62 percent combined biomass and organic layer C stock. Though limited in scope, our study reveals that the accelerated growth of alder in alpine tundra systems is resulting in a loss of the organic soil horizons containing large quantities of C and N.

Conclusions

This study investigated the impact of alder shrub biomass and cover on soil organic layer properties in graminoid, alpine tundra in Interior Alaska. Though limited in site selection by terrain and infrastructure, this study sought to fill a need for in situ measurements of the impact of alder shrub growth on the landscape. We found that increasing alder biomass and cover—at five sites that have exhibited increasing growth over the past four decades—was associated with a reduction in organic layer thickness, C, and N, resulting in net C loss among organic soil horizons and biomass. However, we found no difference in ALT associated with increasing alder growth. Given the scale and speed of Arctic shrubification, more studies are urgently needed to understand the impacts of shrub growth on C sequestration and permafrost degradation, especially because it may vary for other species of shrubs.

Acknowledgments

We thank S. Sackett (Fairbanks ABoVE field office) and NSF Toolik Field Station for logistical assistance during field campaigns.

Disclosure statement

No potential conflict of interest was reported by the authors.

Data availability statement

The data collected for this study are available at the Oak Ridge National Laboratory Distributed Active Archive Center (https://doi.org/10.3334/ORNLDAAC/2300). Figures and models were made with R v4.2.1: https://cran.r-project.org/bin/windows/base/old/.

Supplementary material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/15230430.2023.2285334.

Additional information

Funding

This work was funded by NASA ABoVE [#NNX15AU16A] (to MLG and CIC). Further, RGJ was supported by NSF-OPP [1836873] (to J. M. Welker), AMW by the Ridge to Reef Research Traineeship [NSF-DGE 1735040], and BM and YK by University of California, Irvine Undergraduate Research Opportunities Program.