Impacts of wildfires on boreal forest ecosystem carbon dynamics from 1986 to 2020

Wildfires significantly change boreal forest ecosystem carbon balance through both direct combustion and post-fire carbon dynamics. Affected vegetation influences soil thermal regime and carbon cycling by impacting the surface energy balance of boreal forests. This study uses a process-based biogeochemistry model to quantify carbon budget of North American boreal forests during 1986–2020 based on satellite-derived burn severity data. During the study period, burn severity generally increases. Fires remove ecosystem carbon of 2.4 Pg C and reduce net ecosystem production (NEP) from 32.6 to 0.8 Tg C yr−1, making the forest ecosystems lose 3.5 Pg C, shifting a carbon sink to a source. The canopy’s cooling effect leads to lower soil temperature and lower net primary production due to lower nitrogen mineralization and uptake. Post-fire NEP decreases from 1.6 to 0.8 Tg C yr−1. This reduction accounts for 50% of the simulated NEP when the effects of fire-affected canopy are not considered. Our study highlights the importance of wildfires and their induced-canopy changes in soil thermal and ecosystem carbon dynamics of boreal forests.


Introduction
Boreal forests store more than one-third of global terrestrial carbon (Kasischke and Stocks 2000).This large carbon pool is vulnerable to fire disturbance (Helbig et al 2016) and will potentially release large amounts of carbon into the atmosphere through direct combustion emission and post-fire decreased ecosystem production (Kurz and Apps 1999, Amiro et al 2006, Yin et al 2020, Zhao et al 2021).The warming climate under anthropogenic forcings has increased burned area and burn severity in boreal forests (Gillett et al 2004, Iglesias et al 2022, Zheng et al 2023), suggesting there are positive feedbacks between fires and global climate (Moubarak et al 2023).
Wildfires influence carbon dynamics in boreal forests by removing vegetation carbon via crown fires and through consumption of surface organic matters, including soil organic carbon (Turquety et al 2007, Turetsky et al 2011, De Groot et al 2013, Rogers et al 2015).With changing ecosystem structure, soil moisture, soil temperature, carbon dynamics can be significantly affected (Sullivan et al 2011, Li et al 2017).Soil organic matter (SOM) combustion can release a large quantity of carbon into the atmosphere, accounting for up to 90% of the total carbon emission for severe fires (Walker et al 2018).In addition to direct combustion emissions, post-fire soil respiration can also impact the carbon budget by changing soil temperature and moisture (Kulmala et al 2014).The postfire net primary production (NPP) can also decrease significantly due to the loss of the canopy (Zhao et al 2021).
After fire, the loss of canopy and SOM influences soil temperature, in turn, affecting carbon dynamics (Drüke et al 2021, Martín Belda et al 2022, Oppen et al 2022, Xu and Zhuang 2023).Boreal forests underlain with permafrost are especially vulnerable to this effect since their physical properties and carbon storage are sensitive to soil temperature changes (Hayes et al 2014, Johnston et al 2014).After fire, the loss of canopy can alter the albedo of land and impact the local radiative energy budget.This is because canopy can reflect solar radiation and reduce the short-wave radiation directly absorbed by the land (Gu et al 1999).Canopy also influences the near-surface viscosity and change the sensible heat flux between land and the atmosphere (Gu et al 1999).Observations indicate that sub-canopy temperature is 2.1 • C lower than free atmosphere in summer across Europe (Zellweger et al 2020, Haesen et al 2021).Canopy presence induces soil temperature changes and therefore the regional carbon budget, especially for post-fire areas where the canopy changes during the recovery process.For example, heterotrophic respiration (R H ) and nitrogen mineralization are directly influenced by soil temperature, impacting the soil carbon pool sizes.In addition, SOM loss depends on burn severity and carbon storage.Fires influence soil thermal dynamics during the forest recovery (Sullivan et al 2011, Zhao et al 2021).Taken together, fires will directly affect carbon release, canopy structure, altering surface energy balance and soil thermal regime, indirectly impacting carbon balance.
Previous studies have modeled fire or canopy's influences on carbon dynamics with limitations.Many of them mainly focus on direct combustion emission (Conard and Ivanova 1997, Amiro et al 2001, French et al 2002), and have not considered the impacts of burn severity because it is unavailable or difficult to obtain (Kasischke et al 2005, Balshi et al 2007).Some studies have tried to overcome these limitations by using process-based models (Zhao et al 2021) but did not consider the impacts of fire-induced canopy changes on soil thermal regime and carbon budget.For studies that have included canopy's influence, many of them have not used explicit burn severity information (Martín Belda et al 2022, Xu andZhuang 2023).In all, carbon modeling studies can be improved by using the data of satellite data of burn severity and fire-induced canopy changes.
This study revises a sophisticated process-based model, the terrestrial ecosystem model (TEM) (Zhuang et al 2002) to quantify the influence of fire and canopy change on regional carbon dynamics in North American boreal forests.The interactions between canopy, fire and soil thermal dynamics are modeled.We incorporate a soil surface energy balance model (Xu and Zhuang 2023) to a previous TEM model (Zhao et al 2021).The revised-TEM is then used to evaluate the spatial and temporal carbon dynamics based on satellite burn severity data during 1986-2020.The revised-TEM is expected to better estimate the soil temperature under the influence of fire-induced canopy changes, the carbon storage changes, combustion emissions, and net ecosystem carbon exchanges.

Burn severity data development
Historical fire burn areas were retrieved from 1986 to 2020 for 45 degrees latitudes northward in North America, based on the data from the Bureau of Land Management & Alaska Fire Service Alaska Interagency Coordination Center, the Canadian National Fire Database, and the United States geological survey combined wildfire dataset (Welty and Jeffries 2020).The Alaska, Canada, and contiguous U.S. perimeters were created by their respective agencies through a combination of methods over the period of interest, with older fires (pre-Landsat) typically requiring field data and aerial imagery, and more recent fires being satellite-derived.All NA fire perimeters were standardized through a series of sequential dissolves to fix issues regarding fire ID uniquity using GIS-based tools available in QGIS or ArcGIS Pro.Fire perimeters were then simplified to a 60-meter resolution to improve run time of subsequent processing steps.
A google earth engine (GEE) script was developed to calculate a pre-and post-normalized burn ratio (NBR) value (Miller and Thode 2007) from within each perimeter, as well as from within a 300 m buffer ring offset 1.5 km from the perimeter.Because Landsat Collection 2 was not yet available on GEE at time of analysis, data were sourced from Landsat Collection 1. Use of Landsat Collection 2 GEE product in a future data product release would provide image calibration and other improvements relative to Collection 1 (Crawford et al 2023).Landsat sensors 4-8 were used after applying a cloud, cloud shadow, snow, and water mask.NBR for each fire was calculated as the ratio between the near infrared (NIR) and short-wave infrared bands for the available Landsat mission: Pre-fire data use the median pixel values of the image collection, with images filtered to the approximate snow-free fire season of 15 June to 15 September, from the two years before the fire.Postfire data use the median pixel values from the image collection 15 June to 15 September the year following the fire.The median pixel values from within the perimeter and buffer ring averaged to create a pre-fire NBR and post-fire NBR for each perimeter and buffer ring.These NBR values were written out of GEE along with each unique fire ID.Methods in the GEE code for masking and pooling the Landsat 4-8 data were sourced from Holsinger et al (2021).
The values in the buffer ring are used to correct for differences between the imagery not related to the fire (e.g.phenology, plant health), as this land area is assumed to have not burned.One caveat is that in some areas that had nearby fires within the same year, there may be some overlap between the buffer ring and another fire perimeter, however these cases are uncommon.The difference between the pre-fire NBR and the postfire NBR is known as the delta NBR (dNBR) and can be used to represent burn severity for a fire: dNBR = NBR prefire − NBR postfire . (2) In addition to calculating dNBR for each individual fire over 1986-2020, data were aggregated at a 0.5 • × 0.5 • grid cell for each year in order to incorporate these data into TEM simulations.

Model description
TEM is a process-based biogeochemical model (Zhuang et al 2002) that quantifies carbon and nitrogen dynamics in terrestrial ecosystems.TEM has been used to simulate fire disturbance impacts on carbon dynamics in North America (Zhuang et al 2002, Zhao et al 2021).In addition to the direct combustion emission, TEM can also estimate carbon dynamics during post-fire recovery (Zhao et al 2021).The foliage recovery is assumed linear for the first 5 years and sigmoid for the following years.In TEM, the soil thermal module uses surface air temperature as upper boundary condition to simulate soil thermal dynamics, which might cause large uncertainties (Haesen et al 2021), especially in summer when solar radiation is very strong (Xu and Zhuang 2023).In this study, we incorporate a surface energy balance scheme into TEM by considering the influences of plant canopy on soil surface temperature.The model is then used to evaluate the impacts of fire disturbance on regional carbon dynamics.More details of this scheme can be found in the supporting information (figure S1).
The calculated dNBR, representing burn severity for fires, is used to estimate the proportion of vegetation and soil carbon consumption in TEM.In the model, fires remove vegetation and soil carbon based on the burn severity (Zhuang et al 2002, Zhao et al 2021).Fires with greater burn severity remove more carbon and nitrogen from the ecosystem.Following fire, soil organic carbon composition can vary more drastically than vegetation carbon due to reduction in litterfall while the vegetation recovers.This reduced vegetation carbon and soil nitrogen storage results in lower ecosystem production during the recovery phase.In addition, fires impact soil thermal dynamics (i.e.soil temperature), together with the canopy's influence on soils, further changing carbon and nitrogen dynamics, such as R H and vegetation nitrogen uptake.

Model validation
The revised-TEM is calibrated using observation data from black spruce ecosystems in Alaska (Zhuang et al 2002, Zhao et al 2021).The canopy energy balance scheme is calibrated using publicly available soil temperature observation data, collected and published by AmeriFlux (Xu and Zhuang 2023).This scheme better estimates the soil thermal regime.The parameters used in this study are from previous studies (Zhao et al 2021, Xu andZhuang 2023).Canopy reflects short-wave radiation and has a cooling effect on the soil surface, especially in summer when the shortwave radiation is the strongest in a year.This influence is close to zero when it is winter (Xu and Zhuang 2023).Here, we compare simulated summer soil surface temperature with observation data from AmeriFlux (Black 2016a, 2016b, 2023, Ueyama et al 2023).Simulated soil surface temperature in July is generally higher than observations when canopy's effect is not considered (figure 1).When canopy's effect is considered, the simulated soil surface temperature is much closer to the observations.We compare the simulated carbon dynamics with field measurements in July 2015 from Köster et al (2017) at two sites burned in 1990 and 2012, respectively (table 1).The model can well estimate the vegetation carbon, soil carbon, soil organic nitrogen, and surface soil temperature.

Regional input data and simulation protocols
Data required to drive TEM includes monthly mean surface air temperature, cloudiness, precipitation, vapor pressure, and surface wind speed, which are from ERA5 from 1986 to 2020 (Hersbach et al 2020).Spatially-explicit soil texture (percentage of silt, clay and sand), elevation, plant function types and annual CO 2 concentrations of the atmosphere are also used (Melillo et al 1993, Zhuang et al 2002).The model is spun-up for 120 years before 1986 with cyclic climate data from 1986 to 2000 to achieve an equilibrium state.Then transient simulations from 1986 to 2020 are conducted for each grid cell at a spatial resolution of 0.5 • × 0.5 • for the North American boreal forests.
Two regional simulations are conducted, with and without fire disturbance considered, respectively.Fire polygons are dissected into each unit with unique fire history and then intersected with each grid cell when considering fire impacts.The output values for each grid cell are area-weighted mean of each fire polygon and no-burn area within the cell (Zhao et al 2021).An extra simulation with fire considered but without a canopy energy balance scheme is also conducted for quantifying the impact of the canopy.

Regional burn area and severity
Most grid cells (0.5 • × 0.5 • ) are burned for less than 15% while the maximum burn proportion in a grid cell can exceed 80%, which mainly occurred in Alaska, Central Canada and Quebec (figure 2).The average dNBR value is 200-400 for the most burned area (figure S2

2008
).The overall area-weighted value is 368.0, which is slightly higher than the estimates in Zhao et al (2021).There is a slightly increasing trend for the annual area-weighted dNBR value (figure S2(b)), indicating the burn severity is increasing in this region.

Fire impacts on soil thermal regime
Compared with the no-fire simulation, simulated mean soil temperature increases (figure 3

Fire impacts on carbon dynamics
When there is no fire, the whole ecosystem acts as a net carbon sink in this region with highest net ecosystem production (NEP) in east Canada and Alaska (figure 4(b)).The regional 35 year mean NEP is 32.6 Tg C yr −1 (table 2).When there is fire, the NEP significantly decreases with many regions changing from a carbon sink to a source (figure 4 3), enhancing R H . On the other hand, fires also significantly remove soil carbon (figure S4(d)), decreasing R H . Overall, fires reduce R H by 31.5 g C m −2 yr −1 on average.When canopy's influence is not considered, TEM estimates the regional mean NEP is 1.6 Tg C yr −1 (table 2).It is slightly lower than the estimate by Zhao et al (2021), which might be due to the relatively higher burn severity in our data.When canopy's influence is considered, the mean NEP under fire in this region is 0.8 Tg C yr −1 (table 2) while the combustion emissions are 67.7 Tg C yr −1 (table 2).That means although the boreal forests show a carbon sink during the post-fire recovery, the whole region still acts as a carbon source and will release 2340.5 Tg C into the atmosphere during the study period (66.9 Tg C yr −1 ).Compared with the no-fire simulation, the whole ecosystem carbon storage is reduced by 3481.5 Tg C or 99.5 Tg C yr −1 (figure S4).

Impacts of fire-induced canopy changes on soil thermal regime and carbon dynamics
The fire-induced changes in canopy can significantly influence ecosystem carbon dynamics (Martín Belda et al 2022, Oppen et al 2022, Xu and Zhuang 2023).
Here we compare the simulations between with and without the canopy energy balance scheme in the model in addition to considering fire impacts.
Spatial differences of the simulated monthly mean soil temperature in July (when solar radiation is the strongest) for the period of 1986-2020 due to the influences of the canopy are significant (figure 2).Soil surface temperature at 5 cm soil depth (figure 3(a)) is consistently lower with a magnitude of about 1 • C for most grid cells (Zellweger et al 2020, Haesen et al   3(b)) exhibits a similar pattern to the soil surface with smaller magnitude of changes since the deep soil is less impacted by the air and land surface processes.The regional mean difference is about 0.9 • C for soil surface temperature (figure S3(a)) and 0.6 • C for deep soil temperature, indicating a persistent cooling effect from 1986 to 2020.Both NPP and R H decreased due to the canopy's cooling effect (figure 5).The magnitude of the changes shows a similar pattern to the soil temperature changes (figure 3), which means stronger cooling led to higher carbon changes.The regional differences of NPP and R H (figure S5) remain negative for each year during 1986-2020.In all, the whole ecosystem has a decreased NEP of 0.8 Tg C yr −1 (table 2) since the decreased regional mean NPP is higher than R H , resulting in more carbon being released to the atmosphere due to the canopy's influence.The estimated vegetation carbon and soil carbon are lower (figure S6), making the whole ecosystem lose 670.3 Tg C/year with 370.8 Tg C yr −1 from vegetation and 299.5 Tg C yr −1 from soils during 1986-2020.

Canopy's influence on carbon dynamics
Fire-induced canopy changes significantly affect carbon dynamics both directly and indirectly.First, the fire-induced canopy changes directly affect GPP, thus NPP.Second, the canopy changes affect soil surface energy balance and soil temperatures, influencing soil respiration (R H ) and nitrogen mineralization and plant nitrogen uptake, thus NPP.While our previous analysis indicated that R H decreases due to canopy's cooling effect on soil temperatures (Jiang et al 2016), the effects of fire-induced canopy changes have not been analyzed.Here we show that nitrogen mineralization and uptake decreased due to the canopy's cooling effect.The reduced nitrogen uptake inhibits GPP and NPP (figure 5(a)), resulting in a net carbon source.In addition, the spatial patterns of NEP are consistent with soil temperature changes (figure 3).
Our simulated NPP and R H are both lower considering fire-induced canopy changes comparing with the simulations without considering canopy changes (table 2).The estimated decrease of NEP with 0.8 Tg C yr −1 accounts for about 50% of the originally simulated NEP without considering fire-induced canopy effects (table 2).The lower ecosystem carbon storage due to canopy's influence (figure S6) causes a lower combustion emissions (table 2).Overall, canopy's cooling effects significantly affect GPP and NPP and regional carbon dynamics (Martín Belda et al 2022, Oppen et al 2022, Xu and Zhuang 2023), especially when fire disturbance exists.

Comparison with previous estimates
In comparison with previous studies (table S1), our estimated combustion emissions in Alaska agree well with Wiedinmyer and Neff ( 2007), but higher than Zhao et al (2021) and Kasischke and Hoy (2012), lower than Chen et al (2017), Goetz et al (2012).These differences are potentially due to most of these studies focusing on Alaska while our study covers the boreal forests in whole northern North America.For the Canadian region, our estimates are generally higher (Goetz et al 2012, Zhao et al 2021) but agree with Chen et al (2017).For the whole North America, our estimates of combustion emissions generally agree with some newly-published studies (Phillips et al 2022, Potter et al 2023) but is a little higher than some previous studies (van der Werf et al 2010, Zhao et al 2021).This might be explained by relatively higher burn severity in our data compared with Zhao et al (2021).

Uncertainties and limitations
Using dNBR to represent burn severity to estimate fire impacts on carbon might have resulted in uncertainties of carbon dynamics since other factors such as moisture, elevation and time of burn are not considered in this study (Tan et al 2007, Kasischke and Hoy 2012, Zhao et al 2021).For example, it tends to be wetter and more difficult to burn in lower elevation forests.Additionally, the relationship between dNBR and combustion proportion is established based on black spruce dominated forests while North American boreal forests also include white Y Xu et al spruce or pines and other forest types.This might bias burn severity estimations since black spruce tends to show higher dNBR than white spruce with the same burn severity (Rogers et al 2015).For simplicity, this study did not distinguish between these two types but more detailed vegetation types of specific burn severity is needed to better simulate the regional forests response to fire (Amiro et al 2001, Stinson et al 2011).Our model was originally calibrated with field data from forests in Alaska.It might also induce uncertainties when the parameters are applied for whole boreal forests in North America (Zhao et al 2021).
In developing our burn severity data, we used a burn severity product as a (Loboda et al 2018).This dataset does not cover all of North America boreal forests, just the western half of the region and its temporal coverage ends in 2015.Thus, our carbon analysis is based on more spatially and temporally complete burn severity data that are expanded from Loboda et al (2018).
While we have considered vegetation and moss regrowth and soil carbon recovery and their impacts on soil thermal and carbon dynamics (Zhuang et al 2002), we have not explicitly modeled the shift of vegetation type after fire (Gewehr et al 2014, Stuenzi et al 2022).Severe fires might shift dominant vegetation type from slow-growing black spruce to fast growing deciduous trees, resulting in a difference of carbon storage (Mack et al 2021).In addition, the post-fire recovering process of moss (Zhuang et al 2002) and soil surface organic matter recovering and their impacts on soil thermal regime were modeled with a generic formulation of the recoveries.However, due to lacking information on regional moss dynamics, such as moss thickness, the regional simulations have not considered these effects in a spatially explicit manner.Finally, the burn severity might influence the presence of permafrost (Lucash et al 2023), which has not been evaluated in this study while our focus is on examining how fires-induced canopy changes affect soil temperature changes in upper soil layers, in turn, influencing carbon dynamics.These complex vegetation and soil dynamics and their interactions after fire shall be incorporated into future ecosystem carbon modeling for this region.

Conclusions
Using satellite-derived burn severity data, our quantification shows that boreal forest wildfires significantly reduce ecosystem carbon storage through combustion emissions at 67.7 Tg C yr −1 with total 2.4 Pg C during 1986-2020 and reduce post-fire ecosystem production from 32.6 Tg C yr −1 to 0.8 Tg C yr −1 .As a result, wildfires cause the North American boreal forests to lose 3.5 Pg C during the study period and change the boreal forests from a carbon sink to a source.With lower soil temperature simulated by considering the canopy influence on soil thermal dynamics, the revised model estimates a lower NPP due to inhibited nitrogen uptake and mineralization and lower heterotrophic respiration due to lower soil temperature.Despite their counteractive effects, fireinduced canopy changes lead to a reduced post-fire NEP from 1.6 Tg C yr −1 without considering canopy's influence to 0.8 Tg C yr −1 with considering canopy's influence during 1986-2020.This reduction accounts for 50% of the originally simulated post-fire NEP that does not consider canopy influences.This study highlights the importance of considering forest burn severity and fire-induced canopy changes in regional carbon dynamics in boreal forests of North America.
(a)), which agrees with previous studies (Key and Benson 2006, Allen and Sorbel

Figure 3 .
Figure 3. Differences of the simulated mean soil temperature with and without canopy's influence in July (with canopy simulation minus without canopy simulation) and fire (with fire simulation minus without fire simulation) at 5 cm depth (a), (b) and 20 cm depth (c), (d) from 1986 to 2020.
(b)) (Zhao et al 2021) since burning heats the soil surface and removes the ground moss layer which serves as a heat insulator.Deep soil temperature change follows the spatial pattern of surface soil temperature but with a smaller magnitude (figure 3(d)).Overall fires increase mean soil surface temperature about 0.1 • C for burn areas in 2020 (figure S3(b)).

Figure 4 .
Figure 4. Spatial pattern of carbon dynamics during 1986-2020: (a) cumulative NEP considering fires; (b) cumulative NEP without considering fires; (c) differences of NEP with and without considering fires (a minus b); (d) time-mean ecosystem carbon storage (vegetation carbon plus soil carbon) considering fires; (e) ecosystem carbon storage without considering fires; (f) differences of carbon storage with and without considering fires (d minus e).

Table 1 .
Comparison between observed and modeled variables.
Figure2.Fire burn proportion and severity (delta normalized burn ratio, dNBR value) in North American boreal forests from 1986 to 2020: (a) burn proportion (units: %), (b) burn severity (mean dNBR).The black lines mark the boundary of boreal forests.
), acting as a net carbon sink.With the fire disturbance, there is an opposite trend, resulting in a net carbon source.With and without fire, NPP, R H and NEP are synchronous with each other (figureS4).Fires reduce NPP by 154.8 g C m −2 yr −1 on average, which agrees with previous studies(Hicke et al   2003, Sparks et al 2018).Fires increase soil temperature (figure

Table 2 .
Regional mean NPP, RH, NEP and combustion emissions from three simulations.