Tundra plant above-ground biomass and shrub dominance mapped across the North Slope of Alaska

We developed empirical models for predicting AGB from Landsat NDVI using plant biomass harvests from a large network of field sites, and then predicted regional AGB by applying these empirical models to Landsat NDVI composite mosaics covering the study area. We quantified uncertainty arising from sampling and sensor calibration errors using a Monte Carlo framework. Independent field data were not available to fully validate our predictions of plant and shrub AGB, so we provide additional lines of evidence supporting these predictions. First, we visually compared our predictions against spatial patterns of plant and shrub AGB that were visible in high resolution satellite and aerial imagery. Second, shrub AGB should be highest in areas with the tallest and most extensive shrub canopies, so we compared our predictions of shrub AGB with (a) field measurements of shrub canopy height and (b) a regional shrub canopy cover map. Third, field studies have shown that plant AGB and shrub dominance are often higher in shrub tundra than graminoid tundra (Shaver and Chapin 1991, Riedel et al 2005), so we tested whether modeled plant AGB and shrub dominance similarly differed among regional plant communities. Finally, we examined how AGB varied along regional air temperature gradients.

We relied primarily on published field measurements of plant and shrub AGB to develop empirical relationships between AGB and Landsat NDVI (Raynolds et al 2002, Walker et al 2003a, Copass et al 2007a, Copass et al 2007b, Walker et al 2012, Vankoughnett and Grogan 2016). These field measurements were derived from clip harvests at tundra sites on the North Slope (n = 14), the Seward Peninsula (n = 7), and in northern Canada (n = 3). With the exception of one Canadian site (Vankoughnett and Grogan 2016), the measurements were made as part of the NSF-sponsored Arctic Transitions in the Land-Atmosphere System project (Walker et al 2003a) and subsequent International Polar Year Greening of the Arctic initiative (Walker et al 2012). Biomass was harvested during the period of peak AGB (mid-June to mid-August). Live and standing-dead vascular plants were clipped at the moss surface in 5–12 quadrats that were distributed on transects or a large plot at each site. Quadrats were typically 20 × 50 cm. Lichen and the green portion of mosses were also harvested. Samples were sorted by plant functional type (PFT; for example, shrub, graminoid, forb, moss, lichen), dried, weighed, and then expressed as dry biomass per unit area (kg m⁻²). We acquired or computed the mean and standard error (SE) of plant and shrub AGB among quadrats at each site ($\overline {{\rm{AGB}}}$, SE_AGB). Plant AGB included the total live and standing dead AGB of all PFTs. Shrub AGB included deciduous and evergreen species with woody stems. Finally, we supplemented the collection of field sites with four unvegetated sites that were situated in mountainous areas around the region. These sites were identified using high resolution satellite, aerial, and oblique imagery available on Google Earth. Together, this data set provided estimates of plant and shrub AGB at 28 sites (figure 1(c)). These sites did not constitute a random or systematic sample of the study area, but did encompass a broad range of tundra plant communities. Additional details are provided in the supplementary material available at stacks.iop.org/ERL/13/035002/mmedia.

We used Landsat data to: (1) create a seamless multi-band composite mosaic of the study area, and (2) determine NDVI at each of the 28 locations for which we had estimates of plant and shrub AGB. Specifically, we used Landsat 5 and 7 surface reflectance data sets from the United States Geological Survey that had been atmospherically corrected using the Landsat Ecosystem Disturbance Adaptive Processing System (Masek et al 2006). We accessed and post-processed the Landsat data using Google Earth Engine (Gorelick et al 2017).

Creating a seamless composite mosaic of the study area first involved collecting all available scenes from Landsat 5 (n = 361) and Landsat 7 (n = 1360) that were acquired over the region from late June through late August between 2007 and 2016. Next, we masked pixels that had been flagged as water, snow, clouds, or cloud shadows using the function of mask algorithm (Zhu and Woodcock 2012, Zhu et al 2015) and took additional steps to mask residual water, clouds, and cloud shadows (Ju and Masek 2016, Pekel et al 2016). After extensive masking, we: (1) computed the NDVI for each scene, (2) bias-corrected the NDVI from Landsat 5 to match Landsat 7 (figure S2; Ju and Masek 2016), (3) computed the 80th percentile NDVI values for each pixel from all available scenes, and then (4) identified the scene with the NDVI closest to this percentile. Finally, we created a multi-band Landsat composite mosaic of the study area by extracting spectral reflectance information from these 'peak greenness' scenes on a per pixel basis.

In addition to generating the composite mosaic, we also computed Landsat peak summer NDVI for each field site. We used the same masking and compositing approach, but instead of considering the entire 2007–2016 period, we selected Landsat scenes that spanned a four-year window around when a cluster of sites was sampled. A four-year window balanced data availability against the need to ascertain NDVI concurrent with field sampling. We computed the mean and standard error of Landsat peak summer NDVI for all unmasked pixels within a 50 m radius buffer around each field site ($\overline {{\rm{NDVI}}}$, SE_NDVI). We tested multiple buffer sizes (50 m, 100 m, or 200 m) and found that they had minimal effect on mean NDVI at these sites. After generating the Landsat composite mosaic and estimating peak summer NDVI at each field site, we mapped plant and shrub AGB across the study area.

We modeled plant AGB (kg m⁻²), shrub AGB (kg m⁻²), and shrub dominance (%) over the study area by first developing empirical AGB–NDVI relationships at the field sites and then applying these relationships to regional Landsat NDVI mosaics. We tracked uncertainty in these predictions using a Monte Carlo approach (Harmon et al 2007, Berner et al 2017). Specifically, we focused on: (1) uncertainty in the AGB–NDVI models resulting from sampling errors in both variables within each site, and (2) uncertainty in the regional NDVI given uncertainty in the sensor calibration (Markham and Helder 2012). Our approach involved repeatedly (n = 1000) permuting $\overline {{\rm{AGB}}}$ and $\overline {{\rm{NDVI}}}$ at each site, fitting nonlinear (exponential) models to the permuted field data, and then predicting regional AGB by combining each nonlinear model with a permutation of the regional NDVI mosaic. We fitted these models using the nonlinear least squares function in R (R Core Team 2015) and, after predicting AGB, computed shrub dominance for each permutation. Finally, for each of these plant attributes, we computed the pixel-wise median, 2.5th percentile, and 97.5th percentile among the suite of permutations. This yielded both a best estimate (median) and 95% confidence interval (CI) for each plant attribute at 30 m resolution for the North Slope. These CIs incorporate several important sources of uncertainty, but do not provide a full account of potential errors (for example, atmospheric correction) and thus are a conservative estimate of uncertainty. We implemented this modeling routine on a high-performance computing cluster at Northern Arizona University. Furthermore, we pre-processed, analyzed and visualized these data using a combination of ArcGIS 10.5 (ESRI, Redlands, CA), ArcScene 10.5, and the R packages dplyr (Wickham and Francois 2015) and raster (Hijmans and van Etten 2013). Additional details are provided in the supplementary material available at stacks.iop.org/ERL/13/035002/mmedia.

We visually compared our predictions of plant and shrub AGB with spatial patterns visible in a high resolution GeoEye-1 satellite image (1.8 m) and an oblique aerial image. The GeoEye-1 image was acquired during summer 2010 and obtained from the University of Minnesota Polar Geospatial Center. The oblique aerial image was acquired from a helicopter during the summer of 2001 as part of an earlier study (Tape et al 2006).

We also compared our predictions of shrub AGB against field measurements of shrub canopy height made on large plots (250 × 250 m) at 25 sites that were widely distributed along the Chandler and lower Colville Rivers as well as the Dalton Highway (figure 1(a); Duchesne et al 2016). These field measurements focused on shrubs > 50 cm tall and were made during 2010 and 2011. We computed the mean and SE of shrub AGB among pixels within a 100 m radius buffer around each site and then assessed the strength of the association between mean shrub AGB and mean shrub canopy height using Spearman's correlations (r_s).

We additionally compared our predictions of shrub AGB against a new regional map of shrub canopy cover (Macander et al 2017) in which shrub canopy cover (%) was mapped across ~12.5 Mha of the North Slope at a 30 m spatial resolution by linking point-intercept field measurements with Landsat imagery and other geospatial datasets. We performed a pixel-wise comparison between modeled shrub AGB and canopy cover for the area of overlap between the data sets.

We summarized plant biomass by plant community type using the Ecosystems of Northern Alaska data set (Jorgenson and Heiner 2003), which was derived using Landsat and depicts the spatial extent of 36 ecosystem classes on the North Slope. We selected and consolidated ecosystem classes to those representing shrub, tussock, sedge, and mixed tundra plant communities and then computed the mean and total plant and shrub AGB within each plant community. We tested for differences in mean and total AGB among plant communities using a boot-strapping approach based on the Monte Carlo simulations. Additional details are provided in the supplementary material, as are tabular summaries for each consolidated and non-consolidated ecosystem class (tables S3, S4, S5).

We examined how plant biomass changed along regional gradients in both monthly mean air temperatures (T_avg) and the summer warmth index (SWI). The WorldClim-2 data set provided estimates of monthly T_avg (1970–2000) at 1 km resolution (Fick and Hijmans 2017), which we resampled to 30 m resolution. We computed the SWI as cumulative monthly T_avg exceeding 0 °C (Jia et al 2003, Walker et al 2003a). We then calculated mean plant AGB, shrub AGB, and shrub dominance at 0.1 °C increments along the regional temperature gradients and assessed the association between each mean plant attribute and regional air temperatures using Spearman's correlations.

We estimated the total plant AGB of the North Slope as 112 (91, 135) Tg. Shrub AGB accounted for 53 (32, 75) Tg, or 43 (28, 53)% of the regional plant AGB. Plant AGB averaged 0.74 (0.60, 0.89) kg m⁻² over the study area and shrub AGB averaged 0.35 (0.21, 0.49) kg m⁻² (table S1). Plant and shrub AGB respectively exceeded 1.09 and 0.62 kg m⁻² across 10% of the region (figure S3). The biomass maps revealed a general decline in plant AGB, shrub AGB, and shrub dominance from the Brooks Foothills to the Beaufort Coastal Plain, as well as substantial variability in AGB at finer spatial scales (figures 4 and 5). For instance, the biomass maps resolved shrub-dominant areas with high AGB along streams and rivers, as well as nearby upland areas and bogs with lower AGB and less influence from shrubs. These verdant riparian corridors and sparser upland areas are clearly visible in the high resolution satellite and aerial images (figures 4 and 5).

Areas with taller and more extensive shrub canopies tended to have higher predicted shrub AGB than areas with lower, sparser shrub canopies (figure 6). There was a strong positive relationship between modeled shrub AGB and mean shrub canopy height measured at 25 field sites (r_s = 0.88, p < 0.001). Modeled shrub AGB appeared to increase more rapidly with mean shrub canopy height at sites with the shortest canopies; however, this likely reflects that canopy height was only measured on shrubs > 50 cm tall, whereas modeled shrub AGB includes shrubs of all sizes (figure 6(a)). In addition to the relationship with shrub height, there was also a strong positive relationship between modeled shrub AGB and the regional shrub canopy cover map (r_s = 0.76, p < 0.001, figure 6(b)).

Zoom In Zoom Out Reset image size

There were notable differences in plant and shrub AGB among plant communities (figure 7 (a)). The mean plant AGB was significantly (p < 0.05) higher in shrub tundra (~0.97 kg m⁻²) than in tussock tundra (~0.71 kg m⁻²) or sedge tundra (~0.50 kg m⁻²), as was the mean shrub AGB (figure 7(b)). Tussock tundra had higher mean plant AGB than sedge tundra, though mean shrub AGB did not differ significantly between these plant communities (p = 0.08). Mean shrub dominance in shrub tundra (~53%) was also higher than in sedge tundra (~34%), but did not differ significantly from tussock tundra (~43%, p = 0.14). In total, mixed tundra accounted for the largest percentage of regional plant AGB (~41%), followed by shrub (~31%), sedge (~16%), and then tussock (~10%) tundra (figure 7(b)). Regional shrub AGB was distributed similarly between mixed tundra (~40%) and shrub (~36%) tundra, with much smaller percentages occurring in sedge (~12%) and tussock (~10%) tundra (p < 0.05).

Zoom In Zoom Out Reset image size

Regional variation in plant and shrub AGB was strongly influenced by summer air temperatures (figure 8). Specifically, mean plant AGB, shrub AGB, and shrub dominance increased persistently with increasing mean monthly air temperatures during June, July, and August (r_s =0.96–0.98, p < 0.01). These attributes also increased notably along the regional gradient in SWI (all r_s=0.97, p < 0.01), though showed considerably weaker relationships with mean monthly air temperatures during spring, fall, and winter (r_s = −0.11–0.82; figure 8(b)). For instance, the mean plant AGB nearly doubles (~0.49 kg m⁻² versus ~0.90 kg m⁻²) between areas with the coldest and warmest mean air temperature during June, while shrub AGB nearly quintuples (~0.10 kg m⁻² versus ~0.47 kg m⁻²) along this thermal gradient (figure 8(c)). Additionally, shrub dominance nearly doubles (~28% versus ~50%) between the coldest and warmest areas (figure 8(c)).

Zoom In Zoom Out Reset image size

The maps of plant and shrub AGB for the North Slope of Alaska that we present here complement the circumpolar Arctic biomass map that was developed by linking biomass harvests with NDVI derived from AVHRR (Raynolds et al 2012). Several associated studies examined trends in plant AGB using AVHRR time series (Jia et al 2003, Epstein et al 2012), yet ours is the only other satellite-derived map of plant AGB that covers the North Slope. Strikingly, regional estimates of the total plant AGB differed by just 4% between our map (~112 Tg) and the circumpolar map (~116 Tg), and spatial patterns of plant AGB were broadly similar after aggregating our map to the coarser spatial resolution of the circumpolar map (r_s = 0.63, p < 0.01, figures S4 and S5). Although our regional maps cover just ~3% of the Arctic tundra biome covered by the circumpolar map, they provide estimates of AGB at ~171 000 times greater detail (30 m versus 12 400 m resolution). As a result, they resolve considerable landscape variation in both plant and shrub AGB, as is evident through comparison with high resolution imagery (figures 4 and 5). Our maps thereby provide a unique, management-relevant resource that can inform ecological monitoring and analysis of tundra plant communities on the North Slope.

Despite our focus on AGB in this study, we recognize that below-ground biomass (BGB) comprises a large but highly variable proportion of total plant biomass in tundra communities. For instance, root-to-shoot ratios range from approximately 1:2 for shrub tundra and 3:1 for tussock tundra to 20:1 for wet sedge-moss tundra (Bliss 2000). Global syntheses report root-to-shoot ratios of ~4.8 (Mokany et al 2006) to ~6.6 (Jackson et al 1996) for the tundra biome. We did not try to map plant BGB in this study given the paucity of the available BGB data and the wide variation in estimates for various vegetation types. However, using this collection of root-to-shoot ratios, we estimate that regional plant BGB on the North Slope roughly totals 536–737 Tg. Combining plant AGB and BGB suggests that plant biomass totals 647–848 Tg, which represents a carbon stock of 323–424 Tg assuming that biomass is 50% carbon. This plant biomass carbon stock is small (< 10%) relative to the amount of carbon stored in permafrost soils to 1 m depth on the North Slope (~6930 Tg C; Hugelius et al 2013). Nevertheless, plant biomass influences near-surface permafrost dynamics by shading and insulating the soil horizon (Walker et al 2003b, Blok et al 2010), as well as by mediating the fire regime (Berner et al 2012, Loranty et al 2016). Future efforts are needed to more accurately quantify and map plant BGB in tundra ecosystems, as well as to evaluate interactions among above-ground and below-ground ecosystem carbon pools (Webb et al 2017).

Our approach to mapping AGB was informed by prior studies demonstrating strong relationships between tundra plant AGB and NDVI measured using hand-held spectrometers (Shippert et al 1995, Boelman et al 2003, Jia et al 2003) and satellite-borne spectrometers, typically AVHRR (Jia et al 2003, Walker et al 2003a, Raynolds et al 2012). The strong nonlinear relationship we observed between plant AGB and Landsat NDVI was consistent with these prior studies, though unique in that we demonstrated the relationship using Landsat and a relatively large network of tundra field sites. We also showed there was a strong relationship between shrub AGB and Landsat NDVI, thus expanding on earlier work using NDVI from hand-held spectrometers (Jia et al 2003). Shrubs are frequently the tallest plants in Low Arctic tundra and their canopy characteristics exert a strong influence over tundra spectral reflectance (Boelman et al 2011, Loranty et al 2011, Bratsch et al 2017). The NDVI is often considered to be an indicator of plant photosynthetic capacity (Tucker 1979) and likely tracks variation in tundra plant AGB given co-variation among plant community photosynthesis, leaf area, and AGB (Boelman et al 2003, Street et al 2007, Walker et al 2012). Our study highlights Landsat NDVI as a useful metric of spatial variation in plant and shrub AGB within the tundra biome. Furthermore, our results indicate that satellite observations of arctic greening likely reflect widespread increases in plant AGB and shrub dominance within this biome (Tape et al 2006, Epstein et al 2012, Frost et al 2014, Frost and Epstein 2014).

Our approach to mapping AGB depended on open access to the Landsat archive and the ability to efficiently process the large volumes of information in this archive. For this study we mined spectral information from over 2000 Landsat scenes using Google Earth Engine (Gorelick et al 2017). As a comparison, Beck et al (2011) mapped shrub cover over roughly the same region using seven Landsat scenes that were carefully mosaiced and inter-calibrated, yet which nevertheless exhibited some spectral irregularities and were missing data where clouds were present. Open access to the Landsat archive (Wulder et al 2012), automated cloud and shadow masking (Zhu and Woodcock 2012, Zhu et al 2015), and high-throughput computing (Gorelick et al 2017) are enabling rigorous new analyses of Earth's valuable land and water resources at spatiotemporal extents not previously achieved (Hansen et al 2013, Ju and Masek 2016, Pekel et al 2016).

We found that areas with taller and more extensive shrub canopies typically had higher predicted shrub AGB. Similarly, shrub canopy height and AGB scale allometrically among individual willow, alder, and birch from northern Alaska and Siberia (figure S6; Berner et al 2015). Regional warming contributed to modeled increases in shrub height and cover over the past half century in our study region, particularly within riparian areas (Tape et al 2006, Naito and Cairns 2011, Tape et al 2016a) where repeat aerial photographs show conspicuous shrub expansion since the late 1940s (Tape et al 2006). Shrub AGB is related to both shrub height and cover, thus the changes noted above were likely accompanied by an increase in regional shrub AGB (Jia et al 2003); however, the magnitude of change remains uncertain and warrants future investigation.

We found that plant communities previously mapped as shrub tundra had higher predicted plant AGB and shrub dominance than communities mapped as tussock or sedge tundra. The differences that we observed among plant communities were similar to several prior studies (Walker et al 2003a, Riedel et al 2005). For instance, we estimated that plant AGB in shrub communities averaged 0.97 (0.81, 1.18) kg m⁻² across the region and that shrubs typically comprised ~53% of this plant AGB. Plant AGB averaged 1.20 (1.00, 1.39) kg m⁻² in a shrub community near Ivotuk, where shrubs similarly comprised 51% of plant AGB (Riedel et al 2005). Plant AGB and shrub dominance near Ivotuk were also highest in the shrub community, followed by tussock and sedge communities where shrubs composed 8%–39% of plant AGB. The differences in shrub dominance that we observed among plant communities was also consistent with differences in shrub cover predicted from Landsat (Beck et al 2011). These observations add credence to our predictions of plant and shrub AGB, while also highlighting the fact that future shrub encroachment into graminoid tundra will likely be accompanied by an increase in plant AGB (Pearson et al 2013).

We found that mean plant AGB nearly doubles between the coldest and warmest parts of the North Slope, while shrub AGB nearly quintuples. Consequently, shrub dominance was highest in areas with the warmest summers. Prior field studies similarly demonstrated that tundra sites with the warmest summers tend to have the highest plant AGB and shrub dominance (Epstein et al 2008, Walker et al 2012), along with the tallest and fastest growing willows (Walker 1987). Remote sensing studies have further shown that tundra NDVI (Raynolds et al 2008, Berner et al 2013) and shrub cover (Beck et al 2011) tend to increase along spatial gradients of increasing summer air temperatures. Our current study corroborates and advances these prior findings by linking field and remote sensing measurements to simultaneously demonstrate the pronounced temperature response of plant AGB, shrub AGB, and shrub dominance over a large region.

Given these observed patterns and recent ecosystem modeling (Pearson et al 2013, Yu et al 2017), we anticipate that continued warming will further increase plant productivity, AGB, and shrub dominance on the North Slope. Nevertheless, long-term warming could drive a net loss of carbon from these tundra ecosystems if the decomposition of soil organic matter (carbon loss) is stimulated to a greater degree than plant productivity (carbon gain; Mack et al 2004, Belshe et al 2013, Abbott et al 2016, Commane et al 2017). The greening of North Slope tundra ecosystems (Guay et al 2014, Ju and Masek 2016) was accompanied by increasing ecosystem respiration, such that the region was a net source of CO₂ to the atmosphere from 2012–2014 (Commane et al 2017). Future efforts are needed to better understand how continued warming might affect climate feedbacks associated with tundra ecosystem carbon dynamics (Mack et al 2004, Abbott et al 2016), as well as with changes in tundra surface albedo and evapotranspiration (Pearson et al 2013).

Our new maps could be a useful tool for regional wildlife habitat modeling. Shrub height and canopy volume are strong predictors of the habitat suitability for moose (Alces alces gigas) and snowshoe hares (Lepus americanus) on the North Slope since these shrub attributes influence potential foraging and cover from predators (Zhou et al 2017). Shrub expansion and enlargement in riparian areas led to more suitable habitat for moose and snowshoe hares over the past century, likely causing these herbivores to expand their ranges into increasingly verdant parts of the North Slope (Tape et al 2016a, Tape et al 2016b). On the other hand, arctic greening has been linked to declining caribou populations across tundra ecosystems in North America (Fauchald et al 2017), potentially due in part to reduced foraging as shrubs outcompete lichen (Joly et al 2009, Fraser et al 2014, Fauchald et al 2017). Changes in shrub dominance can also impact migratory birds (Boelman et al 2014) and arthropod communities (Rich et al 2013, Sweet et al 2015). Our new maps of plant AGB, shrub AGB, and shrub dominance provide insight into the habitat suitability for herbivores, insectivores, and their associated predators (for example, lynx, wolves) in this rapidly changing region.