Multi-decadal patterns of vegetation succession after tundra fire on the Yukon-Kuskokwim Delta, Alaska

Most fire activity on the YKD has occurred in the Izaviknek and Kingaglia uplands (IKU), a 3870 km² region of eolian surficial deposits that supports extensive mesic vegetation (Péwé 1975, Wolfe et al 2009, Jorgenson and Roth 2010; figure 1). The IKU composes 5.0% of the YKD ecoregion (Nowacki et al 2003); terrain consists of gentle, rolling hills that range 10–65 m ASL. Unburned vegetation is dominated by lichens, shrubs (mainly <0.5 m height), and sedges corresponding to bioclimate subzone E of the Circumpolar Arctic Vegetation Map (CAVM Team 2003, Raynolds et al 2019). The nearest available long-term instrumental record, approximately 75 km southeast of the study area at Bethel Airport, indicates a mean annual air temperature of −0.7 °C, mean summer (June–August) temperature of 12.2 °C, and mean annual precipitation of 378 mm (1981–2010 base period; NCDC 2019). The IKU is underlain by discontinuous permafrost (Jorgenson et al 2008); permafrost is widespread in upland tundra but is absent from drainages. Approximately 20% (789 km²) of the IKU has experienced fire since 1953, with about half of all fire occurring in 2015 (AICC 2019; figure 2).

Zoom In Zoom Out Reset image size

We established field plots spanning gradients of fire history during 20–31 July 2017 (24 plots) and 8–16 July 2018 (18 plots) (figure 3, table 1). Field data were collected using a vegetation point-intercept (VPI) sampling approach (Karl et al 2017). We supplemented our dataset with 16 VPI plots sampled during 15 June–17 August 2016 (Ludwig et al 2018). The three field campaigns used broadly similar methods, with minor differences to plot layout prompted by the need to address multiple research and land management objectives concurrently in a remote study area.

Zoom In Zoom Out Reset image size

For plot allocation, we used a systematic grid system adapted from the Alaska Landbird Monitoring Survey (Handel and Cady 2004). We assigned fire history categories (years since last fire, or no known fire) to ALMS points within 90 km of the Bethel, AK airport using historical fire perimeters (AICC 2019). Following review of late 1980s Landsat and recent high-resolution imagery, we also mapped and included for plot allocation several 1985 fire perimeters that were represented only as point locations in the historical database. We selected points for field visits using a stratified-random approach, with priority given to (1) points within NASA's Arctic Boreal Vulnerability Experiment remote-sensing flightlines, and (2) points in 2015 burn scars (i.e. sampled 2–3 years post-fire) in order to support assessment of burn severity for related projects.

At most plots, we sampled vegetation along three 30 m linear transects radiating at 120° intervals from the plot center, according to the US Bureau of Land Management's Assessment, Inventory, and Monitoring Program protocol (AIM; Toevs et al 2011). At locations sampled 2–3 years post-fire, most plots instead consisted of four parallel 20 m transects oriented from east to west and spaced 5–10 m apart, following guidance from the Fire Effects Monitoring and Inventory System protocol used to assess the severity of recent fire (Lutes et al 2006). The 16 supplemental plots were sampled using either the AIM layout (6 plots), or along a single randomly oriented 100 m transect (10 plots). At all transects, vegetation was sampled at 0.5 m or 1 m intervals by recording all species 'hit' by the beam from a vertically mounted laser pointer. At most plots, we recorded all vascular plant hits if multiple layers were present. At the supplemental plots, only one hit per point was recorded for each species (6 plots) or only the first hit (i.e. top of canopy) was recorded (10 plots). At all plots measured, only single hits per point were recorded for bryophytes and lichens. We also measured soil thaw depth (42 plots) and maximum shrub height (35 plots) at 5 m intervals along transects; we omitted these measurements at a few plots due to time constraints. Finally, we excavated a soil pit to 40 cm depth near the center of each plot and recorded the thickness of the surface organic layer (no measurement at supplemental plots).

We calculated three species cover metrics from the VPI data, for all plots where the field sampling methods allowed:

• 1.  
  Top-hit cover (THC) is the percentage of points where a species was the first hit. THC represents the top-of-canopy and could be calculated for all 56 plots; the THC of all species sums to 100% per plot. We also calculated THC for litter (dead vegetation and detached leaves), water, and bare ground.
• 2.  
  Any-hit cover (AHC) is the percentage of points where a species occurred. The total AHC for a plot can exceed 100% (i.e. in multi-layered vegetation); we calculated AHC for live vegetation for 46 plots.
• 3.  
  Multi-hit cover (MHC) is the total number of hits for a species, divided by the number of sample points. MHC is a metric of canopy density and can include multiple hits of a species at a single point. MHC thus can exceed 100%; we calculated MHC for live vegetation at 38 plots.

We also calculated these metrics for nine plant functional types (PFTs): deciduous shrub, evergreen shrub, graminoid, forb, Sphagnum, other acrocarpous moss, pleurocarpous moss, liverwort, and lichen. We aggregated these PFTs into broader classes (e.g. all bryophytes) for some applications. Acrocarpous (hereafter, unbranched) mosses are primarily early-successional pioneers, while pleurocarpous (branched) mosses tend to occupy undisturbed sites. We did not distinguish shrub height classes because shrubs rarely exceeded 50 cm height.

Using the AHC summaries, we developed a plot ordination using nonmetric multidimensional scaling (NMDS) for 38 plots with sufficiently detailed species data for nonvascular plants. We aggregated several nonvascular taxa to genus (e.g. Sphagnum, Cladonia) or species-complex (Cladina arbuscula/mitis) because we could not differentiate all species in the field. Taxa that occurred at only one plot were omitted. We used the Bray–Curtis dissimilarity index as the distance measure and restricted the number of ordination dimensions to the minimum necessary to achieve stress values <0.2. We then tested for differences among fire treatments by overlaying 95th percentile confidence-interval ellipses around treatment centroids within the ordination space. We produced a secondary matrix containing plot-level means of soil thaw depth, surface organic layer thickness, and MHC for each PFT and used the envfit function (Oksanen et al 2016) in R (R Core Team 2018) to fit bi-plot vectors for each variable within the ordination. We used MHC because it is sensitive to both the composition and structure of the canopy.

We also calculated the converse Simpson's diversity index, a mathematical measure that characterizes species diversity in a community by the number and abundance of species present (Simpson 1949, Smith and Wilson 1996, Goswami et al 2017). We calculated the Simpson's index using AHC values recorded for each taxon at each plot. Plots sampled 2–3, 10–11, and 45–46 years post-fire were aggregated in the ordination and diversity analyses. We used the AHC metric for the ordination and Simpson's index because it includes all taxa present, but does not overly weight multi-layered PFTs (e.g. shrubs).

We tested for significant differences in the AHC of PFTs, Simpson's diversity index, and soil variables among fire history classes using one-way analysis of variance (ANOVA). Where ANOVA revealed significant differences (p < 0.05), we performed adjusted post hoc multiple comparisons using Tukey's honestly significant difference tests to identify significant differences (adjusted p < 0.05) by class.

We used the VPI data and 60 spectral predictors derived from Landsat data to develop continuous-cover maps for four simplified tundra PFTs at 30 m resolution: shrub, herbaceous, bryophyte, and lichen. We performed modeling in Google Earth Engine (Gorelick et al 2017) using a random forest data-mining approach (Breiman 2001), following the framework used by Macander et al (2017) in northern Alaska.

Spectral predictors consisted of percentile-based surface reflectance composites for 6 bands and 6 indices (e.g. Normalized Difference Vegetation Index) derived from Tier 1 Landsat 7 Enhanced Thematic Mapper + and Landsat 8 Operational Land Imager data from June to September, 2016–2018 (supplemental table S1 is available online at stacks.iop.org/ERL/15/025003/mmedia). We developed surface reflectance composites on a per-pixel basis for five percentiles: 10%, 25%, 50% (median), 75%, and 90%. Although seasonal timing of imagery was not explicitly considered, the percentile-based approach effectively captures spectral conditions during different phenological periods while minimizing spatial artifacts. We masked pixels with a water occurrence of ≥50% in the Joint Research Commission global water product (Pekel et al 2016).

We reserved 20% of the training data (12 plots) as an independent validation dataset for each PFT model, and ran models separately for the validation dataset and the remaining training data. Finally, we ran models for the full study area and evaluated model performance separately using independent predictions for the validation dataset and 'out-of-bag' predictions from the training data. The predictions were made in Google Earth Engine and error metrics were computed in R.

Two NMDS ordination dimensions were sufficient to meet the stress threshold (stress score = 0.184); axis 1 represents a species composition gradient from early to late succession, and axis 2 captures a structural gradient from nonvascular- to shrub-dominated vegetation (figure 4). Field plots formed three nonoverlapping clusters—unburned, 2–3 years post-fire (hereafter, recent burns), and ≥10 years post-fire (hereafter, historical burns)—indicating strong patterns of correspondence between species cover and fire history. Plots sampled 2–3 and 45–46 years post-fire were least and most similar to unburned plots, respectively. Unburned plots were very similar to each other, while burned plots displayed much more variability. The Simpson's diversity index was lowest in recent burns and increased during post-fire succession (figure 5).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Vascular plant diversity was generally low at all plots, and many species were present with high frequency regardless of fire history (table 2, supplementary table S2). For example, the deciduous shrubs Betula nana and Vaccinium uliginosum, and the evergreen shrubs Empetrum nigrum, Rhododendron tomentosum, and V. vitis-idaea accounted for most of the shrub cover regardless of fire history. Graminoid cover consisted almost entirely of sedges and tended to be highest in historical burns; grasses were rare except in historical burns. The tussock-forming sedge Eriophorum vaginatum was present at most plots, but it was never abundant. Forbs were uncommon regardless of fire history. Nonvascular species composition, however, varied strongly among plots with different fire histories.

In unburned tundra, fruticose lichens (e.g. Cladina rangiferina, Flavocetraria cucullata) formed extensive, thick, mixed-species mats (median AHC 83.3%) that formed the canopy at nearly half of sample points (figure 6). Vascular cover was dominated by evergreen shrubs. Bryophyte cover in unburned tundra was dominated by Sphagnum; Sphagnum remained common throughout post-fire succession and reached its highest cover values in the oldest burns.

Zoom In Zoom Out Reset image size

Live vegetation remained discontinuous for several years after 2015 fires; median THC of live vegetation 1, 2, and 3 years post-fire was 24.7%, 28.8%, and 51.4%, respectively. Moribund tufts of Cladina were common in recently burned areas; their bleached color suggested the algal photobiont had been killed by heat. Two yr post-fire, the median AHC of bryophytes was about half that recorded in unburned tundra and was mostly composed of Sphagnum. The unbranched mosses Polytrichum juniperinum and Ceratodon purpureus appeared 2–3 years post-fire; Ceratodon was restricted to early successional sites, but Polytrichum remained common throughout post-fire succession. By 3 years post-fire, lichen cover remained very low but deciduous shrubs, evergreen shrubs, and bryophytes became much more abundant.

In plots sampled 10–11 years post-fire, the combined abundance of all PFTs (median MHC 116.7%) was much higher than in recent burns and was similar to unburned tundra. The dominant PFTs were bryophytes, evergreen shrubs, and deciduous shrubs. Most nonvascular cover consisted of unbranched mosses and Sphagnum.

Plots sampled 32 years post-fire were characterized by a multi-layered plant canopy, and the combined abundance of all PFTs (median MHC 183.9%) was the highest recorded for any class. Evergreen shrubs were the most abundant vascular PFT, followed by deciduous shrubs and graminoids. Nonvascular vegetation was dominated by bryophytes, primarily unbranched mosses. Cover of branched mosses was generally low except in burns ≥32 years old. Lichen abundance was similar to that observed 10–11 years post-fire.

In plots sampled 45–46 years after fire, a multi-layered canopy was present at most points, but median MHC (151.6%) was somewhat lower than in 32 years old burns. The most abundant vascular PFTs were evergreen and deciduous shrubs. Nonvascular vegetation was dominated by bryophytes, primarily Sphagnum. Lichen cover was much higher (median AHC 29.4%) than in younger burns.

Maximum shrub height showed little variation across the fire history gradient. The maximum shrub height recorded was only 48 cm, and plot-level mean maximum heights ranged 6.8–19.7 cm (figure 7). Maximum shrub heights tended to be highest in 10–11 years old and 32 years old burns.

Zoom In Zoom Out Reset image size

Across all plots, thaw depth ranged 40.8–62.1 cm (figure 8). Thaw depth was shallowest in unburned tundra (median 41.6 cm), and deepest 10–11 years post-fire (52.3 cm). Thaw depth was highly variable in plots measured 2–3 years post-fire, and least variable in unburned plots.

Zoom In Zoom Out Reset image size

Surface organic thickness ranged 7–45 cm and was most variable in recent burns (figure 8). There were, however, no clear patterns of correspondence between organic thickness and burn history. We observed exposed mineral soil only at a handful of convex microsites that were severely burned in 2015.

Based on the independent validation dataset and 'out-of-bag' error metrics, model performance was best for lichens, followed by shrubs and bryophytes (figure 9). Model performance was poor for herbaceous, likely because live herbaceous cover was generally low and was similar among all fire history classes. Fire history was not an input in the modeling, but modeled PFT abundance displayed relationships with fire history similar to those seen in the field data. For example, shrub cover was highest in 1985 (32 years old) burns and declined somewhat in older burns.

Zoom In Zoom Out Reset image size

The mapped cover of PFTs displays strong patterns of correspondence with fire perimeters (figure 10). Shrub cover decreases abruptly in the footprint of 2015 fire scars, and increases abruptly in burns 32 years and older. A similar pattern is evident for bryophyte AHC, except that bryophyte AHC is also high in unburned, poorly drained terrain in the southeastern fringe of the study area. The lichen model displays the most striking patterns, with very high cover in unburned uplands and low cover in both recent and historical burns.

Zoom In Zoom Out Reset image size

4.1.1. Vegetation

4.1.2. Soils and permafrost

Investigations of tundra fire are frequently hampered by the low resolution and data gaps for historical fire perimeters (Kasischke et al 2002) and the lack of reliable fire severity data. For example, the Landsat-based Monitoring Trends in Burn Severity (MTBS) product (Eidenshink et al 2007) only exists for fires since 2002 on the IKU. MTBS did, however, capture the extensive 2015 fire footprint, of which >91% was classified as low or moderate severity. The lichen-dominated fuels and the scarcity of exposed mineral soils observed in recent firescars suggest that high-severity fire has been limited in extent on the IKU. We also know that fire perimeters invariably contain unburned patches, and inadvertent field sampling of such patches could confound analysis of post-fire succession. However, the high mortality of lichen after fire, and the persistence of Polytrichum in historical burns, provided obvious field characters by which we avoided unburned inclusions.

Uncertainties related to our chronosequence-based sampling could arise from local differences in pre-fire vegetation and soil conditions, raising the question of whether fire activity was concentrated in areas where antecedent conditions differed from unburned areas. We focused our study on upland tundra and avoided sampling wetlands and streambanks. We consider it unlikely that variability within upland tundra in the IKU was sufficient to influence broad-scale fire behavior. The quantitative cover maps presented here, coupled with high-resolution topography and soil moisture metrics, could inform future modeling efforts to characterize burn probability across the landscape.

The field dataset was relatively small due to logistical constraints, but each plot comprised 90 point measurements of vegetation spread across a 700 m² area. In addition, our spatial models of PFT abundance show patterns of correspondence with fire history that corroborate the field observations, even though fire history was not a predictor in the models. Important sources of variability in post-fire succession, however, likely arise from differences in the landscape setting and severity of specific fires. While successional patterns likely vary according to landscape position and fire severity, our field observations are representative of the contiguous upland areas where most fire activity has occurred. The emergence of quantitative vegetation mapping techniques (He et al 2019) and fire severity metrics appropriate for tundra regions (Loboda et al 2013) provide a way forward for assessing tundra recovery after fire.