Patterns of Ecosystem Structure and Wildfire Carbon Combustion Across Six Ecoregions of the North American Boreal Forest

Introduction

The boreal forest is one of the largest biomes on earth, covering almost 1.9 billion hectares and encompassing ∼30% of the global forested area (Brandt et al., 2013). From a global perspective, the most significant and critical function of the boreal forest is its ability to sequester and store carbon (C), as it contains approximately one third of terrestrial C stocks (Pan et al., 2011). This globally important biome is becoming increasingly vulnerable to change as the climate continues to warm.

One of the most rapid pathways through which boreal forests can be altered is with changes to the fire regime. Climate warming and drying has led to an intensification of boreal forest fires (Balshi et al., 2009; Flannigan et al., 2009), with large increases in the average annual area burned over recent decades (Coops et al., 2018). Further increases, up to a factor of five, are expected by the end of the century (Balshi et al., 2009; Boulanger et al., 2014). Increasing fire extent, frequency, and severity could alter boreal forest C storage (Bond-Lamberty et al., 2007; Walker et al., 2019): from the historical net accumulation of C from the atmosphere over multiple fire cycles, to a net loss, which in turn would cause a positive feedback to global climate warming (Oris et al., 2013; Li et al., 2017). Understanding the structure and function of these fire-prone ecoregions is needed to quantify the role of fire in the global C cycle and its feedbacks to climate change.

The accumulation of C in boreal ecosystems increases with stand age or time after fire, but the rate of accumulation depends on vegetation composition and environmental conditions that impact net primary productivity and decomposition (Luyssaert et al., 2008; Jonsson and Wardle, 2010). Black spruce [Picea mariana (Mill.) B.S.P.] is one of the dominant tree species in boreal North America and is generally less productive and slower growing than jack pine (Pinus banksiana Lamb.) or deciduous trees (Subedi and Sharma, 2013; Alexander and Mack, 2016). The majority of organic C sequestered in black spruce boreal forests resides in the soil organic layer (SOL; Ping et al., 2008; Johnson et al., 2011). Combustion of the SOL dominates C loss during fire, typically accounting for 80–90% of total C emitted in boreal forests of Alaska and Canada (Boby et al., 2010; Rogers et al., 2014; Walker et al., 2018b). Despite its importance, the greatest uncertainty in modeling C combustion from wildfire is estimating the SOL component (French et al., 2004), which requires information on SOL bulk density (g cm^–3) and C fraction (%) to model C stocks (g C m^–2). Both the rate of C accumulation and SOL characteristics are likely to be ecoregion-specific as a consequence of differences in long-term climate, geological and biogeographical history, soil development, and parent materials (Houle et al., 2017).

Fire severity can be defined as the loss of above- and belowground organic material that occurs as a direct consequence of fire (Keeley, 2009). Given that the majority of combustion takes place in the SOL (Boby et al., 2010; Rogers et al., 2014; Walker et al., 2018b), estimates of SOL burn depth and post-fire residual SOL depth have been used as metrics of fire severity (Greene et al., 2004, 2007; Johnstone and Kasischke, 2005; Johnstone and Chapin, 2006; Walker et al., 2018a). Total C combusted from fire is also a widely used metric of fire severity (Boby et al., 2010; Walker et al., 2018b), but variations in above- and belowground C combustion, C combusted relative to pre-fire C, and the proportion of belowground C combusted relative to total C combusted are also important metrics for understanding fire severity and its ecological impacts. In this study, we compiled data on pre-fire C pools, C loss from fire, and the potential drivers of these C metrics from 417 recently (2004–2015) burned sites within 18 fires. We also collected depth-wise measurements of soil characteristics from 167 of these sites plus an additional 110 sites (27 recently burned and 83 control sites with no historical record of burning in the last 50–75 years). Sites were established by five independent research projects, sampled one year after fire, and distributed across six ecoregions of North America’s western boreal forests (Figure 1 and Supplementary Table S1). Using this unique and spatially extensive dataset, we addressed the following questions:

FIGURE 1

Figure 1. Distribution of field sites across six ecoregions in the western boreal forests of North America. See Supplementary Table S1 for details on sites used for ecoregion models and sites used for soil property models.

(1) How do pre-fire aboveground (stand age, density, basal area, proportion of black spruce, biomass, and C pools) and belowground (SOL depth and C pools) forest attributes differ among fire-prone ecoregions?

(2) Do aboveground, belowground, and total C pools increase over time at a similar rate among ecoregions?

(3) In the SOL profile, what is the relationship between depth and (i) bulk density, (ii) C fraction, and (iii) C stocks, and does this vary between ecoregions in a way that is dependent on total SOL thickness?

(4) Do fire severity and C combustion (i.e., residual SOL depth, burn depth, above- and belowground C combustion, proportional combustion) vary among ecoregions?

Materials and Methods

Study Areas and Data Acquisition

We collated data from 527 sites established by five independent research projects located in six ecoregions in the western North American boreal forest (Figure 1 and Supplementary Table S1). Sites were established in the northern ecoregions of Alaska Boreal Interior, Boreal Cordillera, Taiga Plains, and Taiga Shield, and in southern ecoregions of Softwood Shield and Boreal Plains. Ecoregions differ in geological history, soils development, parent materials, and mean annual temperatures and precipitation (US EPA, 2015). We collated site-level data on pre-fire stand structure, stand age, topography, and pre- and post-fire above- and belowground C pools from 417 sites across 18 wildfires. These sites were established by four of the five research projects (Supplementary Table S1), sampled one year after fire, and included in this project if they had all of the measured variables in Table 1. We also collated depth-wise measurements of soil characteristics from 167 of the 417 sites plus an additional 110 sites (27 burned, 83 unburned) across 12 wildfires and eight spatially independent unburned areas. Three out of four of the previously mentioned research projects recorded these depth-wise soil measurements and additional sites were included from a fifth research project (Supplementary Table S1). Sites were chosen to be representative of burned or unburned forests within each ecoregion using remote sensing imagery and fire history records and stratifying by a combination of drainage conditions or fire severity. See references in Supplementary Table S1 for details on study-specific methods of site selection.

TABLE 1

Table 1. Ecoregion specific mean ± standard error for each of the variables of interest associated with pre-fire structure or fire severity and combustion.

Across all research projects, calculations largely followed the methods described in Walker et al., 2018b, but see references in Supplementary Table S1 for details. Briefly, each site was assigned a moisture class based on topography-controlled drainage and adjusted for soil texture and presence of permafrost, ranging from xeric to subhygric, using the method outlined by Johnstone et al. (2008). Stand age at the time of fire was based on tree ring counts from five to ten dominant trees per site. Stem counts, diameter at breast height (DBH) measurements, or study-specific allometric equations were used to calculate tree density (stems m^–2), basal area (m² ha^–1), aboveground biomass (g m^–2), and aboveground C content (g C m^–2). Tree combustion estimates were used to quantify aboveground C combustion These were estimated for each tree as either total percent burned or combustion of structural classes (i.e., foliage, fine branches, large branches, bark). Residual SOL depth was measured at five to 25 points per site. Soil organic layer (SOL) burn depth was estimated based on the height of black spruce adventitious roots above the residual SOL measured on 10 trees per site and calibrated using depth of black spruce adventitious roots in unburned sties (as per Boby et al., 2010; Walker et al., 2018a). In 59 out of the 417 sites, black spruce trees were absent and burn depth was estimated by subtracting the residual SOL from the SOL depth in unburned sites of the same moisture-class (as per Walker et al., 2018a). Pre-fire SOL depth was calculated as the sum of the residual SOL and SOL burn depth. We compiled site-level estimates of residual SOL C, pre-fire SOL C, and SOL C combusted. These were estimated by sampling five to ten soil profiles per site, assessing depth-wise bulk density (g cm^–3), C content (%), and C stocks (g C m^–2), and modeling C stock as a function of SOL depth. Based on the above project-specific estimates, we calculated total C combusted as the sum of above- and belowground C combustion, proportion of pre-fire C combusted as total C combusted divided by total pre-fire C, and proportion of total C combusted attributed to the belowground C pool as belowground C combustion divided by total C combusted.

In addition to the site-level variables, we acquired depth-wise measurements of soil characteristics from 167 of these sites plus an additional 110 sites (27 burned, 83 unburned) (Supplementary Table S1). The additional 27 burned sites had no data on aboveground C pools and the additional 83 unburned sites had no data on C combustion and thus could not be included in our analyses of ecoregion comparisons. Five to ten SOL profiles were collected per site, and bulk density (g cm^–3), C content (%), and C stock (g C m^–2) were assessed from either pre-determined depth increments (e.g., 5 cm) or depth measurements of horizons (e.g., brown moss, fibric, humic). To ensure that all measurements were acquired from the SOL and not mineral soil, we excluded all depth-wise SOL samples (n = 3794) with bulk density >0.75 g cm^–3 and C content <20%. This resulted in 2685 measurements from 1075 residual or unburned SOL profiles nested within 283 sites (Supplementary Table S1). We adjusted each residual SOL profile depth measurements by the site-level estimate of burn depth. Given that each study grouped SOL depth-wise measurements uniquely, we binned the data into 10 cm increments. For example, if the depth of the soil sample or horizon in the pre-fire or unburned SOL profile was <10 cm, we gave it a depth of 10 cm, if the depth was ≥10 cm and <20 cm it was given a depth of 20 cm and so on to a maximum depth of 60 cm.

Statistical Analyses

All statistical analyses were performed using R statistical software version 3.5.2 (R Development Core Team, 2018).

To assess if the structure of these fire-prone ecoregions differed (Q1), we fit linear (LMM) or generalized linear mixed effects models (GLMM) in the package “lme4” (Bates et al., 2015) or the package “glmmTMB” (Brooks et al., 2017). Each model included the fixed effect of ecoregion (six levels) and random effects of fire (18 levels) to allow for different intercepts for each fire scar in order to account for the spatial non-independence of sites within fire scars. Note that because ecoregions cannot be replicated or spatially interspersed, our study is not randomized with replication at the ecoregion level but is randomized with replication within each fire in each ecoregion. We tested response variables of moisture class, pre-fire stand age, density, basal area, proportion of black spruce, aboveground biomass, aboveground C pool, SOL depth, belowground C pool, total C pool, and the proportion of the total C contained in the belowground pool. For stand age, we used a LMM with a normal distribution. For moisture class, we used a GLMM with a Poisson distribution. For each of the proportional response variables we first added a constant of 0.000001 to any zero values and subtracted this constant from any values of one in order to fit a GLMM with a beta distribution and a log-link function. We used GLMM for the remaining variables with a gamma distribution and a log-link function. For these and all mixed models that follow, the significance of fixed effects were determined using maximum likelihood ratio tests comparing the full model to a reduced model and confirmed using Akaike information criterion (AIC, ΔAIC < 2.0; Zuur et al., 2009). Because we had multiple response variables with the same experimental units we corrected the p-values obtained from likelihood ratio tests comparing the full model to a reduced model using the false discovery rates (FDR; Benjamini and Hochberg, 1995). When a significant effect of ecoregion was detected based on the adjusted p-values (p-value <0.05) we used Tukey-Kramer post hoc analysis for multiple comparisons in the R package “emmeans” (Lenth et al., 2019) to test for pairwise differences in marginal means between ecoregions. For these and all LMM or GLMM that follow we assessed bias in model fit by visually inspecting residual versus fitted values.

To assess if aboveground, belowground, and total C pools increased over time at a similar rate between ecoregions (Q2), we fit LMM with response variables of pre-fire aboveground C, belowground C, and total C, all of which were natural log transformed to ensure normality. For each model, we included fixed effects of stand age at the time of fire (i.e., time after establishment), ecoregion (six levels), and their first order interaction and a random effect of fire (18 levels). The importance of fixed effects and p-value adjustments proceeded as described above. When a significant interaction with ecoregion was found, we used a Tukey-Kramer post hoc analysis to complete pairwise comparisons of marginal mean intercepts and slopes between ecoregions.

To determine how (i) bulk density (g cm^–3), (ii) C content (%), and (iii) C stocks (g C m^–2) change with depth in the SOL profile and if these vary with total pre-fire or unburned SOL thickness and ecoregion (Q3), we fit LMM for each soil characteristic. Bulk density and C stocks were natural log transformed for normality. We grouped ecoregions into four categories to ensure sufficient sample size within each unique combination of ecoregion and total pre-fire or unburned SOL thickness category. Taiga Plains (n = 290 soil profiles) and Taiga Shield (n = 192 soil profiles) were left as is, but Alaska Boreal Interior and Boreal Cordillera were grouped as “Alaska” (n = 311 profiles) and the Boreal Plains and Softwood Shield were grouped as “Saskatchewan” (n = 282 soil profiles). We assigned each site to a SOL thickness category (total SOL; 0–20 cm, 20–40 cm, and >40 cm) and modeled soil characteristics using fixed effects of depth within soil profile (depth), ecoregion, total pre-fire or unburned SOL thickness, and the first-order interactions between depth and ecoregion and depth and SOL thickness. For these LMM, we used random effects of soil profile (1075 levels), nested within site (283 levels), nested within fire or unburned area (20 levels). The importance of fixed effects and p-value adjustments proceeded as above. Differences in marginal mean intercepts and slopes of bulk density, C content and C stocks as a function of depth were compared between ecoregions within each SOL thickness class and between SOL thickness classes within each ecoregion using Tukey-Kramer post hoc analysis.

To assess if metrics of fire severity and C combustion varied among ecoregions (Q4), we fit LMM or GLMM, with response variables of residual SOL depth, residual SOL C, burn depth, aboveground, belowground, and total C pools, the proportion of total pre-fire C combusted, and the proportional contribution of the belowground C pool to total C combusted. For each model, we included random effects of fire scar (18 levels) and a fixed effect of ecoregion (six levels). For burn depth, we used a LMM with a normal distribution. For each of the proportional response variables, we fit a GLMM with a beta distribution and a log-link function. We used GLMM for the remaining variables with a gamma distribution and a log-link function. The importance of fixed effects and p-value adjustments proceeded as described above. When a significant effect of ecoregion was detected, we used Tukey-Kramer post hoc analysis to test for pairwise differences in marginal means between ecoregions.

Results

In assessing the pre-fire structure of these fire-prone ecosystems (Q1), we found that stand age, proportion of black spruce, stand basal area, aboveground biomass, aboveground C pools, SOL depths, belowground C pools, total C pools, and the proportion of total C stored belowground differed between ecoregions (Figure 1, Table 1 and Supplementary Table S2). Only site drainage class and pre-fire stand density did not differ between ecoregions (Table 1). In general, the northern Taiga and Alaskan ecoregions were more similar to one another than to the southern boreal ecoregions (Boreal Plains and Softwood Shield). Compared to the southern ecoregions, sites in the northern boreal ecoregions burned at an older age (Figure 2A), were more dominated by black spruce (Figure 2B), had lower aboveground C pools (Figure 2C), and stored proportionally more of the total ecosystem C belowground than aboveground (Table 1). In addition to these north versus south differences, we also observed distinct characteristics of the Taiga ecoregions, where pre-fire SOL depth, belowground C pools (Figure 2E), and total pre-fire C pools were higher compared to the Alaskan or southern boreal ecoregions.

FIGURE 2

Figure 2. Ecoregion differences between (A) stand age, (B) pre-fire proportion of black spruce, (C) pre-fire aboveground C pool, (D) pre-fire belowground C pool, (E) burn depth, (F) residual belowground C pool, (G) aboveground C combusted, and (H) the proportion of total C combustion from belowground. Generalized linear mixed effects models were used to assess the effect of ecoregion and letters represent results of Tukey-Kramer post-hoc pairwise comparisons. See Supplementary Table S2 for model results and Table 1 for means.

Ecoregion differences were also apparent in the rates at which C accumulated (Q2). We found significant interactions between stand age and ecoregion in predicting each of pre-fire aboveground, belowground, and total C (Figures 3A–C, Table 2 and Supplementary Table S3). Pre-fire aboveground C increased slightly in the Taiga Plains and Taiga Shield, but we found no significant relationships between aboveground C pools and time of stand establishment in the Alaskan (Alaska Boreal Interior and Boreal Cordillera) or southern boreal (Boreal Plains and Softwood Shield) ecoregions (Figure 3A and Table 2). Pre-fire belowground C pool increased over time in the Boreal Plains but did not change with time after fire in the other ecoregions (Figures 3B,C and Table 2). Total C pools, the majority of which is from belowground (Table 1 and Figure 3C), significantly increased over time in the Taiga Shield and Boreal Plains (Figure 3A and Table 2). The rate of increase in each of aboveground, belowground, and total C pools was much higher, albeit not significantly different due to large standard errors, in the southern boreal ecoregions compared to the northern.

FIGURE 3

Figure 3. Ecoregion specific effects of time after the previous fire, based on stand age at the time of the most recent fire, on natural log transformed pre-fire (A) aboveground C pools, (B) belowground C pools, and (C) total C pools (g C m^–2). Colors represent ecoregions and lines represent model fits of ecoregion specific significant relationships and 95% confidence intervals. Non-significant relationships are not shown. See Supplementary Table S3 for model output and Table 2 for estimated intercepts and slopes.

TABLE 2

Table 2. Marginal mean intercept and slope estimates and standard errors of pre-fire aboveground carbon (C), belowground C, and total C as a function of stand age for each ecoregion.

In examining the relationships between soil characteristics and depth (Q3), we found that bulk density (g cm^–3), C fraction (%), and C stock (g C m^–2) varied with depth and that these relationships differed between ecoregions and total thickness of the pre-fire or unburned SOL profile (Figure 4, Table 3 and Supplementary Table S4). We found that bulk density (g cm^–3) always increased with depth and within each total SOL thickness categories this increase was always greatest for the Taiga Plains. When comparing within each ecoregion, increases with depth were greatest when total SOL thickness was 0–20 cm (Figure 4 and Table 3). Carbon content (%) decreased with depth in all ecoregions and SOL thickness category, but this decrease was greatest for all ecoregions when SOL thickness was <20 cm. Within each total SOL thickness category the decrease in C content with depth was greatest in the Taiga Plains. In all cases, C stocks (g C m^–2) increased with depth and this rate of increase was greatest for shallow SOL (<20 cm) in the Taiga Plains (Figure 4 and Table 3). Despite having high C stocks, the Taiga Shield had the smallest increase in cumulative C stocks with depth across all SOL thickness categories (Figure 4 and Table 3).

FIGURE 4

Figure 4. Soil properties of bulk density (g cm^–3), carbon content (%) and carbon pool (g C m^–2) as a function of sample depth, total plot SOL thickness, and ecoregion. Lines represent significant ecoregion specific relationships between soil properties and depth in the soil profile and shading represents 95% confidence intervals. Note that points are jittered horizontally at each 10 cm depth increment (width = 2 cm). See Supplementary Table S4 for model results and Table 3 for Tukey-Kramer post-hoc pairwise comparisons of marginal mean intercepts and slopes compared across total SOL thickness categories within ecoregions and across ecoregions within total SOL thickness categories.

TABLE 3

Table 3. Linear mixed effect model estimated marginal mean intercepts and slopes (± standard error) for the effect of sample depth (depth) on the response variables of bulk density, carbon, and carbon pool for each combination of ecoregion (four levels) and total soil organic layer thickness (three levels).

We found ecoregion differences in most of our metrics of fire severity (Q4). Specifically, burn depth, residual SOL depth and C, aboveground C combustion, proportion of C combustion from belowground, and proportion of total pre-fire C combusted were different between ecoregions (Figure 2, Table 1 and Supplementary Table S2). Burn depth (Figure 2E) was slightly higher in the Alaskan ecoregions and residual belowground C pools (Figure 2F) were highest in the Taiga ecoregions, but these differences did not result in any significant ecoregion differences in belowground or total C combustion (Table 1). The most notable difference in fire severity between ecoregions, was that aboveground C combustion was higher (Figure 2F) and the proportion of total C combustion coming from belowground was lower (Figure 2G) in the southern compared to the northern ecoregions.

Discussion

Through compiling a spatially extensive dataset spanning six ecoregions in North America’s western boreal forests, we were able to assess how the structural and functional attributes that relate to above- and belowground C stocks of these fire-prone ecoregions differ. Our study sites captured a broad gradient in pre-fire conditions of tree productivity, stand age, and ecosystem C storage as well as burn depth and C combustion from fire both within and among ecoregions. Given the breadth of our study sites, and that they were chosen to be representative of burned forests within each ecoregion, our results should allow future research to estimate C accumulation, storage, and combustion throughout the western North American boreal forest.

Our regional analysis highlighted pronounced ecoregion differences associated with latitude in boreal forest structural and functional attributes relating to C sequestration and storage. In general, southern ecoregions (Boreal Plains and Softwood Shield) were more similar to one another than to the northern ecoregions, suggesting that differences in long-term climate between these regions are more important than ecoregion differences in geological and biogeographical history, soil development, and parent materials, for predicting ecosystem C dynamics. The age distributions of our examined field sites suggest that there is a higher probability of younger stands burning in southern compared to northern boreal ecoregions. Northern ecoregions also had smaller aboveground C pools, larger belowground C pools, and lower rates of C accumulation. Taken together, these differences suggest that sufficient biomass for fire ignition and spread occurs at an earlier stage in stand development in the southern compared to the northern boreal ecoregions. However, we did not account for potential differences in lightning rates between ecoregions, which is correlated to temperature and precipitation and a dominant driver of total burned area in boreal North America (Veraverbeke et al., 2017).

Although C accumulates slowly in northern boreal forests, the relative low severity and frequency of past fires has allowed soil C pools to accumulate and act as a reservoir for long-term C storage (Harden et al., 2000; Chapin et al., 2006). This has resulted in deep SOL and large belowground C pools in the Alaskan and Taiga ecoregions. These soil C sinks are particularly high in the Taiga ecoregions where peatlands, bogs, and fens are common (Tarnocai et al., 2009; Hossain et al., 2015). However, the soil C pools in northern ecoregions may be more likely to combust with a changing fire regime (Walker et al., 2019) or changing land use (Turetsky et al., 2011a) than they have in the past. The higher frequency of fires in southern boreal regions, in combination with storing proportionally more C above- than belowground and having lower belowground C pools, suggests that the long-term C sink of southern boreal forests is weak compared to the northern ecoregions.

Warmer temperatures increase rates of organic soil decomposition, which release nutrients and increase aboveground C production and total net primary productivity (Chapin et al., 2002; Day et al., 2019). This leads to increased litter fall, root turnover, and an overall higher yield in C input per year in boreal forest soils (Deluca and Boisvenue, 2012). As the southern ecoregions have higher mean annual temperatures and a greater number of growing degree days (>5°C) than the northern ecoregions, these long-term climate parameters are likely important determinants of the higher aboveground C storage and accumulation rates we observed in the southern ecoregions. Similar shifts in C distribution from below- to aboveground pools as temperatures increase have been observed in black spruce forests of Interior Alaska (Kane and Vogel, 2009).

Whether stands established post-fire or post-timber harvest also influence boreal C dynamics (Seedre et al., 2014; Dieleman et al., 2020). Boreal forest stands that originate following harvest tend to have higher rates of aboveground C recovery than those that establish after fire due to reforestation practices or because harvesting activities usually occur on more productive sites (Ilisson and Chen, 2009; Seedre et al., 2014). In the southern boreal ecoregions, almost half of the sites (19 of 44) established following harvest rather than fire, which might also explain why these ecoregions had higher aboveground C pools compared to the northern ecoregions (Dieleman et al., 2020).

Differences in species composition also influence C storage and accumulation rates; black spruce are generally less productive and slower growing than jack pine or deciduous trees (Subedi and Sharma, 2013; Alexander and Mack, 2016). Black spruce was dominant in the majority of sites in northern ecoregions, but in the southern boreal ecoregions, jack pine dominated in 18 out of 44 sites. As the frequency and severity of boreal wildfires continue increasing in these northern ecoregions, stands previously dominated by black spruce are likely to transition to jack pine or deciduous dominance (Johnstone et al., 2010; Whitman et al., 2018; Hart et al., 2019). This change in species composition in combination with direct impacts of climate warming on C storage and accumulation rates, suggests that the future C dynamics of northern boreal forests will resemble those of southern boreal forests.

As expected, soil bulk density and C stocks increased and C content decreased with depth in the SOL profile and the rate at which these soil characteristics changed with depth were dependent on ecoregion and total pre-fire or unburned SOL thickness. The differences in these rates were generally greater between total SOL thickness categories than between ecoregions, suggesting that soil development is more closely related to local drainage conditions than to ecoregion differences. Interestingly, within each total SOL thickness category the Taiga Plains ecoregion had the highest bulk density and C content, potentially due to extensive lacustrine deposits from the last major glaciation that are found throughout this ecoregion (Ecosystem Classification Group et al., 2009). These results highlight the importance of using regionally specific estimates of SOL C storage for predicting C combustion.

Fire severity, measured as the proportion of total ecosystem C combusted, was close to 50% in Alaska and approximately 30% in the other ecoregions (Table 1). Lower proportional combustion in the southern ecoregions is likely due to more C stored aboveground in tree stems that is unavailable for combustion. In contrast, we attribute the lower proportional combustion in the Taiga Plains and the Taiga Shield to deep SOL and large belowground C pools (Table 1) in poorly drained areas where deep soils are protected from burning by permafrost or saturated water table conditions (Turetsky et al., 2011b). Burn depth, residual SOL depth and C, and the proportion of total C combustion from belowground were lowest in the southern ecoregions, emphasizing the high C combustion from more rapidly accumulating aboveground biomass and shallow pre-fire SOL depths in these ecoregions. The fact that burn depth differed but belowground C combustion did not highlights the importance of local or regional conditions in modeling C combustion.

As climate continues to warm and boreal wildfire regimes continue to intensify, the structure and function of the boreal forest is likely to change. Notably, the boreal net ecosystem C balance could switch from net C sink to a net C source (Walker et al., 2019). Understanding C combustion from fires is critical for forecasting such changes and the results we present in this study should allow future research to estimate C accumulation, storage, and combustion throughout the western North American boreal forest. Our findings that C storage, accumulation rates, and sources of C combustion from fire vary latitudinally suggests that the C dynamics of northern boreal forests are likely to shift and resemble the dynamics of southern boreal forests in the future. Specifically, as climate continues to warm northern regions are likely to start accumulating C more rapidly, store more C aboveground, have smaller belowground C stocks, burn more frequently, and emit proportionally more C from the aboveground component than from belowground. Understanding these changes to the structure and function of fire-prone boreal ecoregions are needed to quantify the role of fire in the global C cycle and its feedbacks to climate change.

Data Availability Statement

The data used in this manuscript is archived at the Oak Ridge National Laboratory Distributed Active Archive Center (ORNL DAAC), doi: 10.3334/ORNLDAAC/1744 (Walker et al., 2020).

Author Contributions

XW and MM conceived the study. XW compiled datasets, analyzed the data, and led the writing in collaboration with MM. All authors established the sites and contributed to field data, read, edited, and contributed to the writing of this manuscript.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Funding

This writing of this manuscript and synthesis of data was supported by funding the NASA Arctic Boreal and Vulnerability Experiment (ABoVE) Legacy Carbon grant NNX15AT71A awarded to MM. The original field studies were supported by funding in the United States from NSF DEB RAPID grant #1542150 to MM, NASA Rapid Response grant NNX15AD58G and NASA ABoVE grant NNX15AT83A to LB-C, NASA ABoVE grant NNX15AU56A to BR, SV, and MT, Joint Fire Science Program grant 05-1-2-06 to JJ, NSF grant 0445458 to MM, NSF support to the Bonanza Creek LTER (DEB-0423442); and in Canada from NSERC Discovery Grant funding to JJ and MT; Government of the Northwest Territories Cumulative Impacts Monitoring Program Funding project #170 to JB; NSERC PDFs to ND and CD; GNWT logistical and financial

support through the Laurier-GNWT Partnership Agreement; Polar Knowledge Canada’s Northern Science Training Program funding awarded to Canadian field assistants; SV acknowledges Vidi grant support from the Netherlands Organisation for Scientific Research (NWO).

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/ffgc.2020.00087/full#supplementary-material

References

Alexander, H. D., and Mack, M. C. (2016). A canopy shift in interior alaskan boreal forests: consequences for above- and belowground carbon and nitrogen pools during post-fire succession. Ecosystems 19, 98–114. doi: 10.1007/s10021-015-9920-7

CrossRef Full Text | Google Scholar

Balshi, M. S., Mcguire, A. D., Duffy, P., Flannigan, M., Kicklighter, D. W., and Melillo, J. (2009). Vulnerability of carbon storage in North American boreal forests to wildfires during the 21st century. Glob. Change Biol. 15, 1491–1510. doi: 10.1111/j.1365-2486.2009.01877.x

CrossRef Full Text | Google Scholar

Bates, D., Mächler, M., Bolker, B., and Walker, S. (2015). Fitting linear mixed-effects models using lme4. J. Stat. Softw. 67, 1–48.

Google Scholar

Benjamini, Y., and Hochberg, Y. (1995). Controlling the false discovery rate: a practical and powerful approach to multiple testing. J. R. Stat. Soc. Series B 57, 289–300. doi: 10.1111/j.2517-6161.1995.tb02031.x

CrossRef Full Text | Google Scholar

Boby, L. A., Schuur, E. A., Mack, M. C., Verbyla, D., and Johnstone, J. F. (2010). Quantifying fire severity, carbon, and nitrogen emissions in Alaska’s boreal forest. Ecol. Appl. 20, 1633–1647. doi: 10.1890/08-2295.1

CrossRef Full Text | Google Scholar

Bond-Lamberty, B., Peckham, S. D., Ahl, D. E., and Gower, S. T. (2007). Fire as the dominant driver of central Canadian boreal forest carbon balance. Nature 450, 89–92. doi: 10.1038/nature06272

PubMed Abstract | CrossRef Full Text | Google Scholar

Boulanger, Y., Gauthier, S., and Burton, P. J. (2014). A refinement of models projecting future Canadian fire regimes using homogeneous fire regime zones. Can. J. Forest Res. 44, 365–376. doi: 10.1139/cjfr-2013-0372

CrossRef Full Text | Google Scholar

Brandt, J. P., Flannigan, M. D., Maynard, D. G. I, Thompson, D., and Volney, W. J. A. (2013). An introduction to Canada’s boreal zone: ecosystem processes, health, sustainability, and environmental issues. Environ. Rev. 21, 207–226. doi: 10.1139/er-2013-0040

CrossRef Full Text | Google Scholar

Brooks, M. E., Kristensen, K., van Benthem, K. J., Magnusson, A., Berg, C. W., Nielsen, A., et al. (2017). glmmTMB balances speed and flexibility among packages for zero-inflated generalized linear mixed modeling. R J. 9, 378–400.

Google Scholar

Chapin, F. S., Matson, P. A., and Mooney, H. A. (2002). Principles of Terrestrial Ecosystem Ecology. New York, NY: Springer-Verlag.

Google Scholar

Chapin, F. S. I., Woodwell, G. M., Randerson, J. T., Rastetter, E. B., Lovett, G. M., Baldocchi, D. D., et al. (2006). Reconciling carbon-cycle concepts, terminology, and methods. Ecosystems 9, 1041–1050.

Google Scholar

Coops, N. C., Hermosilla, T., Wulder, M. A., White, J. C., and Bolton, D. K. (2018). A thirty year, fine-scale, characterization of area burned in Canadian forests shows evidence of regionally increasing trends in the last decade. PLoS One 13:e0197218. doi: 10.1371/journal.pone.0197218

PubMed Abstract | CrossRef Full Text | Google Scholar

Day, N. J., Dunfield, K. E., Johnstone, J. F., Mack, M. C., Turetsky, M. R., Walker, X. J., et al. (2019). Wildfire severity reduces richness and alters composition of soil fungal communities in boreal forests of western Canada. Glob. Change Biol. 25, 2310–2324.

Google Scholar

Deluca, T. H., and Boisvenue, C. (2012). Boreal forest soil carbon: distribution, function and modelling. Forest. Int. J. Forest Res. 85, 161–184. doi: 10.1093/forestry/cps003

CrossRef Full Text | Google Scholar

Dieleman, C., Rogers, B. M., Veraverbeke, S., Johnstone, J. F., Laflamme, J., Solvik, K., et al. (2020). Drivers of carbon stocks and combustion in response to wildfire events in the southern boreal forest. Glob. Change Biol. 1–18. doi: 10.1111/gcb.15158

PubMed Abstract | CrossRef Full Text | Google Scholar

Ecosystem Classification Group, Northwest Territories, and Department of Environment and Natural Resources (2009). Ecological Regions of the Northwest Territories: Taiga Plains. Yellowknife, NWT: Dept. of Environment and Natural Resources. Govt. of the Northwest Territories.

Google Scholar

Flannigan, M. D., Krawchuk, M. A., de Groot, W. J., Wotton, B. M., and Gowman, L. M. (2009). Implications of changing climate for global wildland fire. Int. J. Wildland Fire 18, 483–507.

Google Scholar

French, N. H. F., Goovaerts, P., and Kasischke, E. S. (2004). Uncertainty in estimating carbon emissions from boreal forest fires. J. Geophys. Res. Atmos. 109:D14S08.

Google Scholar

Greene, D. F., Macdonald, S. E., Haeussler, S., Domenicano, S., Noel, J., Jayen, K., et al. (2007). The reduction of organic-layer depth by wildfire in the North American boreal forest and its effect on tree recruitment by seed. Can. J. Forest Res. Rev. Can. Rech. Forest. 37, 1012–1023. doi: 10.1139/x06-245

CrossRef Full Text | Google Scholar

Greene, D. F., Noel, J., Bergeron, Y., Rousseau, M., and Gauthier, S. (2004). Recruitment of Picea mariana, Pinus banksiana, and Populus tremuloides across a burn severity gradient following wildfire in the southern boreal forest of Quebec. Can. J. Forest Res. Rev. Can. Rech. Forest. 34, 1845–1857. doi: 10.1139/x04-059

CrossRef Full Text | Google Scholar

Harden, J. W., Trumbore, S. E., Stocks, B. J., Hirsch, A., Gower, S. T., O’neill, K. P., et al. (2000). The role of fire in the boreal carbon budget. Glob. Change Biol. 6, 174–184.

Google Scholar

Hart, S. J., Henkelman, J., McLoughlin, P. D., Nielsen, S. E., Truchon−Savard, A., and Johnstone, J. F. (2019). Examining forest resilience to changing fire frequency in a fire-prone region of boreal forest. Glob. Change Biol. 25, 869–884. doi: 10.1111/gcb.14550

PubMed Abstract | CrossRef Full Text | Google Scholar

Hossain, M. F., Chen, W., and Zhang, Y. (2015). Bulk density of mineral and organic soils in the Canada’s arctic and sub-arctic. Inform. Process. Agric. 2, 183–190. doi: 10.1016/j.inpa.2015.09.001

CrossRef Full Text | Google Scholar

Houle, G. P., Kane, E. S., Kasischke, E. S., Gibson, C. M., and Turetsky, M. R. (2017). Recovery of carbon pools a decade after wildfire in black spruce forests of interior Alaska: effects of soil texture and landscape position. Can. J. Forest Res. 48, 1–10. doi: 10.1139/cjfr-2017-0236

CrossRef Full Text | Google Scholar

Ilisson, T., and Chen, H. Y. H. (2009). Response of six boreal tree species to stand replacing fire and clearcutting. Ecosystems 12, 820–829. doi: 10.1007/s10021-009-9259-z

CrossRef Full Text | Google Scholar

Johnson, K. D., Harden, J., McGuire, A. D., Bliss, N. B., Bockheim, J. G., Clark, M., et al. (2011). Soil carbon distribution in Alaska in relation to soil-forming factors. Geoderma 167, 71–84. doi: 10.1016/j.geoderma.2011.10.006

CrossRef Full Text | Google Scholar

Johnstone, J., and Chapin, F. (2006). Effects of soil burn severity on post-fire tree recruitment in boreal forest. Ecosystems 9, 14–31. doi: 10.1007/s10021-004-0042-x

CrossRef Full Text | Google Scholar

Johnstone, J. F., Hollingsworth, T. N., and Chapin, F. S. III (2008). A Key for Predicting Postfire Successional Trajectories in Black Spruce Stands of Interior Alaska. General Technical Report PNW-GTR-767. Portland, OR: Pacific Northwest Research Station, USDA Forest Service, 37.

Google Scholar

Johnstone, J. F., Hollingsworth, T. N., Chapin, F. S., and Mack, M. C. (2010). Changes in fire regime break the legacy lock on successional trajectories in Alaskan boreal forest. Glob. Change Biol. 16, 1281–1295. doi: 10.1111/j.1365-2486.2009.02051.x

CrossRef Full Text | Google Scholar

Johnstone, J. F., and Kasischke, E. S. (2005). Stand-level effects of soil burn severity on postfire regeneration in a recently burned black spruce forest. Can. J. Forest Res. Rev. Can. Rech. Forest. 35, 2151–2163. doi: 10.1139/x05-087

CrossRef Full Text | Google Scholar

Jonsson, M., and Wardle, D. A. (2010). Structural equation modelling reveals plant-community drivers of carbon storage in boreal forest ecosystems. Biol. Lett. 6, 116–119. doi: 10.1098/rsbl.2009.0613

PubMed Abstract | CrossRef Full Text | Google Scholar

Kane, E. S., and Vogel, J. G. (2009). Patterns of total ecosystem carbon storage with changes in soil temperature in boreal black spruce forests. Ecosystems 12, 322–335. doi: 10.1007/s10021-008-9225-1

CrossRef Full Text | Google Scholar

Keeley, J. E. (2009). Fire intensity, fire severity and burn severity: a brief review and suggested usage. Int. J. Wildland Fire 18, 116–126.

Google Scholar

Lenth, R., Singmann, H., Love, J., Buerkner, P., and Herve, M. (2019). emmeans: Estimated Marginal Means, aka Least-Squares Means. Available online at: https://CRAN.R-project.org/package=emmeans

Google Scholar

Li, F., Lawrence, D. M., and Bond-Lamberty, B. (2017). Impact of fire on global land surface air temperature and energy budget for the 20th century due to changes within ecosystems. Environ. Res. Lett. 12:044014. doi: 10.1088/1748-9326/aa6685

CrossRef Full Text | Google Scholar

Luyssaert, S., Schulze, E.-D., Börner, A., Knohl, A., Hessenmöller, D., Law, B. E., et al. (2008). Old-growth forests as global carbon sinks. Nature 455, 213–215. doi: 10.1038/nature07276

PubMed Abstract | CrossRef Full Text | Google Scholar

Oris, F., Asselin, H., Ali, A. A., Finsinger, W., and Bergeron, Y. (2013). Effect of increased fire activity on global warming in the boreal forest. Environ. Rev. 22, 206–219. doi: 10.1139/er-2013-0062

CrossRef Full Text | Google Scholar

Pan, Y., Birdsey, R. A., Fang, J., Houghton, R., Kauppi, P. E., Kurz, W. A., et al. (2011). A large and persistent carbon sink in the World’s forests. Science 333, 988–993.

Google Scholar

Ping, C., Michaelson, G. J., Jorgenson, M. T., Kimble, J. M., Epstein, H., Romanovsky, V. E., et al. (2008). High stocks of soil organic carbon in the North American Arctic region. Nat. Geosci. 1, 615–619. doi: 10.1038/ngeo284

CrossRef Full Text | Google Scholar

R Development Core Team (2018). R: A Language and Environment for Statistical Computing. Vienna: R Foundation for Statistical Computing.

Google Scholar

Rogers, B. M., Veraverbeke, S., Azzari, G., Czimczik, C. I., Holden, S. R., Mouteva, G. O., et al. (2014). Quantifying fire-wide carbon emissions in Interior Alaska using field measurements and Landsat imagery. J. Geophys. Res. Biogeosci. 119:2014JG002657.

Google Scholar

Seedre, M., Taylor, A. R., Brassard, B. W., Chen, H. Y. H., and Jõgiste, K. (2014). Recovery of ecosystem carbon stocks in young boreal forests: a comparison of harvesting and wildfire disturbance. Ecosystems 17, 851–863. doi: 10.1007/s10021-014-9763-7

CrossRef Full Text | Google Scholar

Subedi, N., and Sharma, M. (2013). Climate-diameter growth relationships of black spruce and jack pine trees in boreal Ontario, Canada. Glob. Change Biol. 19, 505–516. doi: 10.1111/gcb.12033

PubMed Abstract | CrossRef Full Text | Google Scholar

Tarnocai, C., Canadell, J. G., Schuur, E. A. G., Kuhry, P., Mazhitova, G., and Zimov, S. (2009). Soil organic carbon pools in the northern circumpolar permafrost region. Glob. Biogeochem. Cycles 23:GB2023.

Google Scholar

Turetsky, M. R., Donahue, W. F., and Benscoter, B. W. (2011a). Experimental drying intensifies burning and carbon losses in a northern peatland. Nat. Commun. 2, 1–5.

Google Scholar

Turetsky, M. R., Kane, E. S., Harden, J. W., Ottmar, R. D., Manies, K. L., Hoy, E., et al. (2011b). Recent acceleration of biomass burning and carbon losses in Alaskan forests and peatlands. Nat. Geosci. 4, 27–31. doi: 10.1038/ngeo1027

CrossRef Full Text | Google Scholar

US EPA (2015). Ecoregions of North America. Data and Tools. Avaialble online at: https://www.epa.gov/eco-research/ecoregions-north-america (accessed November 25, 2015).

Google Scholar

Veraverbeke, S., Rogers, B. M., Goulden, M. L., Jandt, R. R., Miller, C. E., Wiggins, E. B., et al. (2017). Lightning as a major driver of recent large fire years in North American boreal forests. Nat. Clim. Change 7, 529–534. doi: 10.1038/nclimate3329

CrossRef Full Text | Google Scholar

Walker, X. J., Baltzer, J. L., Bourgeau-Chavez, L. L., Day, N. J., De groot, W. J., Dieleman, C., et al. (2020). ABoVE: Synthesis of Burned and Unburned Forest Site Data, AK and Canada, 1983–2016. Oak Ridge, TN: ORNL Distributed Active Archive Center.

Google Scholar

Walker, X. J., Baltzer, J. L., Cumming, S. G., Day, N. J., Ebert, C., Goetz, S., et al. (2019). Increasing wildfires threaten historic carbon sink of boreal forest soils. Nature 572, 520–523. doi: 10.1038/s41586-019-1474-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Walker, X. J., Baltzer, J. L., Cumming, S. G., Day, N. J., Johnstone, J. F., Rogers, B. M., et al. (2018a). Soil organic layer combustion in boreal black spruce and jack pine stands of the Northwest Territories, Canada. Int. J. Wildland Fire 27, 125–134.

Google Scholar

Walker, X. J., Rogers, B. M., Baltzer, J. L., Cumming, S. G., Day, N. J., Goetz, S. J., et al. (2018b). Cross-scale controls on carbon emissions from boreal forest megafires. Glob. Change Biol. 24, 4251–4265. doi: 10.1111/gcb.14287

PubMed Abstract | CrossRef Full Text | Google Scholar

Whitman, E., Parisien, M.-A., Thompson, D., Flannigan, M., Whitman, E., Parisien, M.-A., et al. (2018). Topoedaphic and forest controls on post-fire vegetation assemblies are modified by fire history and burn severity in the northwestern Canadian boreal forest. Forests 9:151. doi: 10.3390/f9030151

CrossRef Full Text | Google Scholar

Zuur, A., Ieno, E. N., Walker, N., Saveliev, A. A., and Smith, G. M. (2009). Mixed Effects Models and Extensions in Ecology with R. Berlin: Springer Science & Business Media.

Google Scholar