Carbon 5–60 Years After Fire: Planting Trees Does Not Compensate for Losses in Dead Wood Stores

Wildfire is a natural disturbance in many forested biomes, with the loss of carbon to the atmosphere and mortality of trees actively sequestering carbon of global concern as a contribution to climate change. Natural regeneration is often successful at reestablishing a forest in ecosystems adapted to fire, but there is increasing concern that the changing size, frequency and severity of wildfire is causing regeneration failures or inadequate densities of trees that sequester and store carbon following these disturbances. It remains unclear whether the action of planting trees accelerates carbon storage following fire compared to forests established through natural regeneration. The central interior of British Columbia recently experienced multiple years of record-breaking fire activity. Rehabilitation planting focused on reestablishing trees in the managed forest but was also prescribed in previously unmanaged forests to initiate carbon sequestration. Planting is often accompanied by other stand treatments such as salvage harvesting or snag removal and debris clearing to ensure planter safety. Here, we determine carbon recovery and stores in 21 wildfires across a chronosequence from the early 1960s to 2015. We measured above and belowground carbon pools to determine the effect of time since fire and planting treatments on carbon. Tree planting did not increase total ecosystem carbon over time, but rather decreased carbon through the loss of dead wood from site preparation. All carbon pools were affected by time since fire except the mineral soil pool, which was best predicted by soil clay content and coarse fragments positive effects. Live tree carbon increased over time, with more stored in planted stands over 60 years compared to stands that were not planted. Projecting growth to 100 years since fire suggests we may see increasing divergence in carbon stores in planted stands over a full fire-return interval, but these differences remain relatively small [mean (sd): 140.8 (19.6) Mg⋅ha–1 in planted compared to 136.9 (27.5) Mg⋅ha–1 in not-planted stands], with 1.4 Mg⋅ha–1 year–1 sequestered in not-planted compared to 1.5 Mg⋅ha–1 year–1 in planted stands. To meet carbon objectives, replanting trees on average sites in burned forests of BC’s central interior would require preserving the carbon legacy of fire, including dead wood.

Wildfire is a natural disturbance in many forested biomes, with the loss of carbon to the atmosphere and mortality of trees actively sequestering carbon of global concern as a contribution to climate change. Natural regeneration is often successful at reestablishing a forest in ecosystems adapted to fire, but there is increasing concern that the changing size, frequency and severity of wildfire is causing regeneration failures or inadequate densities of trees that sequester and store carbon following these disturbances. It remains unclear whether the action of planting trees accelerates carbon storage following fire compared to forests established through natural regeneration. The central interior of British Columbia recently experienced multiple years of record-breaking fire activity. Rehabilitation planting focused on reestablishing trees in the managed forest but was also prescribed in previously unmanaged forests to initiate carbon sequestration. Planting is often accompanied by other stand treatments such as salvage harvesting or snag removal and debris clearing to ensure planter safety. Here, we determine carbon recovery and stores in 21 wildfires across a chronosequence from the early 1960s to 2015. We measured above and belowground carbon pools to determine the effect of time since fire and planting treatments on carbon. Tree planting did not increase total ecosystem carbon over time, but rather decreased carbon through the loss of dead wood from site preparation. All carbon pools were affected by time since fire except the mineral soil pool, which was best predicted by soil clay content and coarse fragments positive effects. Live tree carbon increased over time, with more stored in planted stands over 60 years compared to stands that were not planted. Projecting growth to 100 years since fire suggests we may see increasing divergence in carbon stores in planted stands over a full fire-return interval, but these differences remain relatively small [mean (sd): 140.8 (19.6) Mg·ha −1 in planted compared to 136.9 (27.5) Mg·ha −1 in not-planted stands], with 1.4 Mg·ha −1 year −1 sequestered in not-planted compared to 1.5 Mg·ha −1 year −1 in planted stands. To meet carbon objectives, replanting trees on average sites in burned forests of BC's central interior would require preserving the carbon legacy of fire, including dead wood.

INTRODUCTION
Forests are an important global store of carbon, with boreal ecosystems sequestering 0.5 Pg·C·year −1 and temperate forests 0.72 Pg C year −1 from 1990 to 2007 (Pan et al., 2011). This terrestrial forest carbon is stored in aboveground wood and vegetation, belowground in the roots of living and dead trees, and belowground in vegetative litter, soil organic layers, and mineral soil horizons. The amount of carbon stored in forests changes over time, as living vegetation establishes, grows, dies, and decays. Large-scale disturbance events, such as insect outbreaks and wildfire result in larger forest carbon emissions. Climate change may increase forest carbon emissions through increased size, severity or frequency of wildfires (Wang et al., 2015;Coogan et al., 2019), with an estimated 2.2 Pg·C·year −1 of carbon emitted from fire globally (1979-2016Van Der Werf et al., 2017). Capturing the impacts of fire on carbon dynamics requires quantifying pyrogenic carbon emissions (Campbell et al., 2007;Walker et al., 2018;Bowman et al., 2020), the carbon legacies from fire (Gough et al., 2007;Santín et al., 2015;Ferster et al., 2016;Furniss et al., 2019;Stenzel et al., 2019) as well as rates of carbon recovery following fire (Seedre et al., 2014;Eskelson et al., 2016;Stevens-Rumann and Morgan, 2019;Kelly et al., 2021). Fire legacies include varying amounts of live, dead, scorched or charred trees, consumed organic soils and exposed mineral soil (Walker et al., 2018). Trees may be relatively undamaged from a fire event, remaining as live trees in a stand, or may be combusted to varying degrees, left standing dead or falling to the ground. The live legacies continue to sequester carbon, while the standing and downed dead wood will serve as attenuated release carbon stores as they decompose (Auclair and Carter, 1993).
Although many forests are adapted to wildfire (Burton et al., 2008;Day et al., 2020), some forested regions are experiencing a decline in forest recovery after fire, either from a lack of trees establishing (Harvey et al., 2014;Stevens-Rumann and Morgan, 2019;Turner et al., 2019), or shifts in forest community composition (Baltzer et al., 2021;Mack et al., 2021). As a result, there is growing concern about the potential for biome shifts, where previously forested land becomes deforested, or slow rates of tree establishment following fire, impacting carbon sequestration potential. Regenerating trees are especially vulnerable to growing season drought (Harvey et al., 2016) and minimum temperatures (Guz et al., 2021). To mitigate this risk of slow or absent tree regeneration and kickstart carbon sequestration after fire, actions such as planting trees are increasingly used to mitigate the carbon lost to fire (Povak et al., 2020). These nature-based solutions to climate change are immediate actions to prevent worst case climate forecasts (Fargione et al., 2018;Drever et al., 2021), but also require ongoing estimation of the carbon consequences of these actions (i.e., Graves et al., 2020). Tree planting after fire may coincide with additional treatments that occur prior to planting. For instance, salvage harvesting, where dead trees are removed shortly after a disturbance for wood products, or site preparations, such as knocking down, piling, and then burning dead trees to mitigate overhead hazards for tree planter safety.
These post-fire treatments affect the biomass remaining after fire, impacting post-fire carbon emissions released directly through combustion (pile burning), or decomposition over time of carbon left on site (Hagemann et al., 2009;Bradford et al., 2012;Kishchuk et al., 2015). Carbon stored in the forest floor, or belowground in roots and in the mineral soil can also be affected by post-fire site preparations if these include soil disturbance, such as mounding or trenching to increase planting microsite availability.
Before the uncertain impact of climate change, BC's subboreal forests would have been expected to naturally regenerate after fire. A stand-replacing fire regime dominates the region, with fire return intervals of 100 years at lower elevations, and 150 years at higher elevations (DeLong, 1998). While a standreplacing regime suggests complete canopy mortality, partial canopy mortality and unburned patches within fire perimeters are common (DeLong and Kessler, 2000;Meyn and Feller, 2006) with a matrix of varying fire severities (Guindon et al., 2020). Lodgepole pine (Pinus contorta var. latifolia Dougl. Ex Loud.), interior spruce [Picea glauca x engelmannii (Moench) Voss)], subalpine fir [Abies lasiocarpa (Hook) Nutt.] and trembling aspen (Populus tremuloides Michx.) dominate these forests and generally regenerate well after fire, particularly the serotinous pine and suckering aspen.
Forests become carbon sources as a result of fire both through the direct combustion of vegetation, wood, and if the carbon emitted from decomposing dead material after the fire exceeds the rate of photosynthesis from newly established, or legacy live vegetation (Auclair and Carter, 1993). The rate at which a forest returns to a carbon sink depends on many factors including previous stand composition , and management history (Dieleman et al., 2020), fire frequency (Gough et al., 2007), the amount of carbon combusted during fire (Walker et al., 2018), rate of decomposition (Auclair and Carter, 1993), and rate of reforestation (Kashian et al., 2006). Logged stands in sub-boreal systems can transition from carbon sources to sinks between 8 and 10 years after harvest (Fredeen et al., 2007), or after 5 years in European hemiboreal (Rebane et al., 2020), whereas burned stands can be more variable in the time required to return to a carbon sink (Amiro et al., 2010). Much of the carbon stored in sub-boreal forests is found belowground (Fredeen et al., 2005;Kranabetter, 2009), requiring carbon accounting of both above and below ground carbon pools to understand the implications of post-fire rehabilitation by tree planting and associated site preparations.
Our study area is at the southern edge of the boreal biome of British Columbia, Canada, where tree planting after fire has emerged as a common action to mitigate the risk of regeneration failure and post-wildfire erosion. Record-breaking wildfires in 2017 and 2018 followed by another significant year in 2021 (BC Wildfire Service, 2021) reinforce the need to quantify the impact of rehabilitation planting on carbon and other aspects of forest recovery. Additional goals of planting programs include establishing healthy timber crops with selected seed and species, fast-forwarding succession, and accelerating carbon storage after fire (Forest Carbon Initiative, 2020). The impact of planting on forest carbon has not been sufficiently understood in these sub-boreal forests to determine which of these goals are being met. Post-fire forest management, such as salvage harvesting, can also influence tree establishment (Povak et al., 2020) and there is a growing need to understand how stands managed and unmanaged after fire compare to each other. It is the rate of postfire carbon recovery that we examine here, and whether it can be accelerated by planting trees. To quantify carbon recovery, we sampled the major carbon pools across a 60-year chronosequence of forests recovering after wildfire of differing severities and project carbon in these stands forward to 100 years post-fire. Specifically, we determine whether planting trees after fire results in greater total ecosystem carbon over time compared to stands that are left to regenerate naturally.

Study Area
The study area is located within the central interior of British Columbia (Figure 1), a sub-boreal landscape with varying topography that is heavily influenced by the proximity to the Coastal and Rocky Mountain chains. The climate is continental with seasonal extremes in temperature; mean annual temperature ranges from 0.3 to 5 • C and mean annual precipitation ranges from 335 to 900 mm, of which 25-50% is snow (Meidinger and Pojar, 1991). Upland forests are dominated by interior spruce, subalpine fir, lodgepole pine, and trembling aspen, with paper birch (Betula papyrifera Marsh.), Douglas-fir [Pseudotsuga menziesii (Mirb.) Franco], and cottonwood (Populus trichocarpa Torr. & Gray) also present.

Research Design
Between 1960 and 2015, over 584,000 ha of forest burned across 952 fires in the study area. We selected 21 fires in this area stratified across a chronosequence of time since fire (5-60 years; Figure 1). Fires were also selected based on accessibility, with most fires containing both planted and notplanted areas, with not-planted areas regenerating naturally after fire (Supplementary Table 1). We verified historical fire boundaries and tree planting history for the older fires using imagery through Google Earth timelapse and field verification (see Skakun et al., 2021; Supplementary Figure 1 for example).

Treatments and Sampling
Within each of the 21 fires selected for this study, we identified potential sampling areas based on access and availability of stands that were either planted or not-planted (naturally regenerated). In these potential sampling areas, we generated random points with package sf in R (Pebesma, 2018), resulting in 39 planted and 36 not-planted stands sampled across the 21 fires. To reduce the variability driven by site productivity, we only considered potential sampling areas if satellite or topographic features were consistent with mesic conditions (i.e., Harper et al., 2005). We did not stratify sampling by previous dominant species due to FIGURE 1 | The study area spans the central interior of British Columbia, with study fire perimeters colored by fire year, and black points indicating plot locations within fires. the wide range of fire years, and many fires predating digital forest inventories.
Fires burning through the central interior of BC contain mosaics of managed and unmanaged forest ( Supplementary  Figure 1), which we captured through our random sampling approach. This resulted in nine stands with management history (harvesting) that pre-dated fire, with the remaining 65 stands without a previous management history (i.e., fire burned mature forest that was then planted or left to naturally regenerate). Two of the previously managed stands were harvested prior to the fire but then not re-planted after the fire, so are considered part of the not-planted treatment, and seven stands that were harvested prior to the fire were replanted after fire (planted treatment; Supplementary Table 1).
Post-fire site treatments in planted areas could include removal of live or dead standing trees (salvage harvesting) or a mechanical or manual tree knockdown, followed by debris management, and possible soil preparation for planting. More recently, post-fire treatments include underplanting, where legacy trees are left on site (Supplementary Table 1). With the wide range in pre-planting site treatments applied in different combinations, we could not determine the impact of any one particular site preparation treatment on post-fire carbon recovery. Rather, we look at the broader treatment of planting vs. not-planting, with planting encapsulating a wide range of associated site preparation activities. Available harvest records for our stands suggest that six stands were salvage harvested after fire and up to 15 had tree knockdown and/or debris pile burning. Four or five stands could be considered underplanted (standing live and dead trees left on site). There were no site preparation records for the remaining 22 planted stands, but field sampling and evidence of forest clearing from satellite imagery if available, suggested many of these stands had similar activitiesknockdown or salvage -after fire.
Within stands, we selected plot locations by walking to a random point and confirming the intended treatment (planted or not-planted) and site productivity (ca. mesic) at the random point. We based our sampling on the Canadian National Forest Inventory (CFIC, 2008) using a nested circular plot design with three plots (11.28, 5.64, and 3.99 m radii) established around plot center. Two 50-m line transects ran through the plot center, with the first following a random compass bearing, and the second perpendicular to the first. At the beginning and end of each transect we established a satellite subplot (3.99 m radius) to increase sampling of tree regeneration when density was low. At each plot we assessed fire severity and measured carbon pools (live trees, dead trees, downed wood, litter, forest floor and mineral soil).

Trees
Live and dead trees were inventoried within the nested circular plots. Large trees (diameter at breast height; dbh ≥ 12.5 cm) were measured within the 11.28 m radius plot, and small trees (dbh < 12.5 cm and >1.3 m tall) were measured in the 5.64 m radius plot. We recorded species, dbh, height, and decay class

Woody Debris
The line intersect technique (Van Wagner, 1968) was used to sample coarse woody debris (CWD) and fine woody debris (FWD). CWD was considered dead wood that is not selfsupporting with a tilt angle < 45 • from the ground plane and with a diameter ≥ 7.5 cm, including otherwise supported or entangled fallen trees not yet resting on the forest floor. FWD was considered dead wood, either small boles or branches, with a diameter < 7.5 cm. Transect distances were corrected for slope. For CWD, diameter, species, decay class (1-5), and depth of char were recorded for each piece intersecting the entire transect. For FWD, pieces were tallied by diameter class (1.1-2.5, 2.6-5.0, and 5.1-7.5 cm) at the first and last 5 m of each transect.

Litter, Forest Floor, and Mineral Soil
Litter, and forest floor were sampled at 11 m from plot center on each transect (totaling four sub-samples per plot). Within a 0.2 m × 0.2 m square frame, litter and forest floor samples were collected for laboratory analysis. Litter was defined as the relatively fresh organics that were easily identifiable plant material such as leaves, wood, or twigs resting on the surface of the forest floor, while forest floor was defined as the organic residues (leaves, branches, bark, and stems) in various stages of decomposition present at the top of mineral soil (Klinka et al., 1997). The four sub-samples for litter or forest floor were bulked together into one sample per plot for laboratory analysis. Litter and forest floor were combined into a single carbon pool (forest floor).
Mineral soil samples, including coarse fragments and roots, were collected by excavation to 0.2 m in depth in one of the four sub-sample locations (0.2 m × 0.2 m) for bulk density determination. The excavated volumes were measured with silica crystals.

Trees
We estimated carbon in standing live and dead trees using species specific allometric equations from Ung et al. (2008) to calculate biomass and then applying a carbon factor of 0.5 (Lamlom and Savidge, 2003). The biomass of standing dead trees was further reduced by applying a species and decay class specific density reduction factor (Harmon et al., 2011) and a decay class specific structural reduction factor (Domke et al., 2011). We also estimated carbon in seedlings (height < 1.3 m) by extending the Ung et al. (2008) biomass allometric equations to trees that had not yet reached breast height. We used the mode of the height class (0.15 cm for 0-30 cm height class, and 80 cm for 31-130 cm height class) and estimated a hypothetical dbh value for each seedling size class by calibrating the resulting biomass estimates against the biomass of similar European seedlings and saplings species (Annighöfer et al., 2016). Based on this, we used a dbh of 0.1 cm for the 0-30 cm height class, and 1 cm for 31-130 cm height class and then applied Ung et al. (2008) speciesspecific allometric equations to estimate biomass, followed by a factor of 0.5 to estimate carbon with a 5% carbon reduction for dead seedlings.

Roots
We calculated roots associated with live trees using the following formula from Li et al. (2003).
Where C ri is the estimated carbon (Mg·ha −1 ) in the roots from live tree i. and C ti is the estimated above-ground carbon of tree i. We did not include root carbon in the dead wood pool, because we could not confidently assess dead tree biomass where site preparation treatments had removed stumps.

Woody Debris
We estimated CWD carbon (Mg·ha −1 ) following methods similar to Russell et al. (2015) using the following formula.
Where V CWD is the volume (m 3 ·ha −1 ) of CWD for species s and size class m, derived using the estimator outlined in Woodall and Monleon (2008); WD sm is the species, s, and size class, m, specific wood density (g·cm −3 ; Harmon et al., 2008); DCRF sm is the species, s, and size class, m, specific decay reduction factor (%; Harmon et al., 2008) SRF is a structural reduction factor for decay class k (%; Fraver et al., 2013) and Cconc is species s and decay class k carbon content (Harmon et al., 2013).
Fine woody debris species were unidentifiable so an average of all present species in the study area was used to calculate the WD m and DCRF m . and we estimated FWD carbon (Mg·ha −1 ) using the following formula: Where V FWD is the estimated volume (m 3 ·ha −1 ) of FWD (Woodall and Monleon, 2008); WD m is the size class, m, specific wood density (g·cm −3 ), DCRF m is the size class, m, specific decay reduction factor and Cconc is the carbon concentration of 50% (Harmon et al., 2013).

Litter, Forest Floor, and Mineral Soil
Litter, forest floor and mineral soil samples were kept cold until they could be air-dried and then oven dried at 70 • C. We removed large woody pieces (pieces that had transitioned from CWD decay class 5 to forest floor) and pieces of charcoal from the forest floor and weighed them separately. We ground the litter and forest floor with a Wiley mill for weighing, sieved mineral soils and separated these into the fine fraction (<2 mm), coarse fragments (>2 mm), roots, and charcoal (>2 mm) for weighing.
We determined total carbon and nitrogen content (g·g) on finely ground mineral (C ff ), forest floor, and litter samples by combustion with a Flash 2000 elemental analyzer (Carter and Gregorich, 2008). We used the hydrometer method to determine particle size by sedimentation on mineral soils samples and used H 2 O 2 pre-treatment for samples with >2.5% carbon (Carter and Gregorich, 2008). We calculated bulk density (BD ff in kg·m 3 ) and soil organic carbon (SOC ff ) of mineral soil for the fine fraction with an assumed specific gravity of the coarse fragments of 2.65 g·cm −3 . The fine fraction mineral soil carbon (Mg·ha −1 ) was calculated using the formula: Where 2 represents the depth of sampling (0.2 m) multiplied by a conversion factor of 10 and CFrags is the coarse fragment content (m 3 ·m 3 ).

Fire Severity
We assessed fire severity in the field by capturing fire effects on trees and soil where possible: depth of burn, effects on roots, scorch height, % ground scorch, and branches burned. We then classed all plots into a five categories of fire severity: low, low-moderate, moderate, moderate-high, and high, based on a combination of these factors. For many of the plots, however, time since fire and forest management (dead tree removal and site treatments) following fire affected the evidence of fire and our ability to accurately classify severity in the field. As a result, we did not incorporate the effect of fire severity in our analysis of the effect of planting on carbon pools for plots from fires older than 1985. For stands burned between 1985 and 2015, we included fire severity using the differenced Normalized Burn Ratio (dNBR) from Guindon et al. (2020), which represents satellite-derived canopy mortality from fire (Supplementary Table 2). We extracted a dNBR value for the plots that burned after 1985 (N = 53) based on intersection of the plot location and dNBR raster. We then validated remote sensing derived fire severity with the field based fire severity assessment (trees and soil) by binning continuous dNBR values (Key and Benson, 2006) and comparing burn severity class to those assigned from the suite of variables assessed in the field.

Environmental Covariates
We examined 18 annual climate variables from ClimateNA (v7.0; Wang et al., 2016). With high correlation between climate variables, we selected growing degree days (GDD) and mean annual precipitation (MAP) to capture the influence of climate on biomass accumulation across the study area (Supplementary Figure 2). We estimated slope (% slope) in the raster package and solar radiation (heat load index; HLI) from the spatialEco package (Evans, 2021) to capture effects of topographic site position. Heat load index accounts for higher heat on southwest slopes and slope steepness (McCune and Keon, 2002). We also used the carbon to nitrogen ratio (C:N), coarse fragment content (CFC; kg m −2 ) and clay content (%) of mineral soil samples collected in the field to estimate site productivity. These environmental covariates were tested for correlation to avoid multicollinearity in predictors. We processed these spatial data in R with raster (Hijmans, 2020) and sf (Pebesma, 2018) packages.

Simulating Future Carbon
To determine the rate of carbon storage at the time scale of the fire return interval, we initiated a SORTIE-ND run and projected the biomass accumulation in live carbon pools up to 100 years since fire for each sampled stand. SORTIE-ND (hereafter SORTIE) is an individual-based, spatially explicit stand dynamics model parameterized for BC's sub-boreal forests. 1,2 In these forests, parameters for the processes of tree establishment, growth and mortality have been estimated for lodgepole pine, interior spruce, subalpine fir and trembling aspen (LePage et al., 2000;Coates et al., 2001;Canham et al., 2004;Greene et al., 2004;Astrup et al., 2008). We updated the sub-boreal SORTIE parameter file (Astrup et al., 2008) with estimates of adult and sapling growth across edaphic gradients (Lilles and Astrup, 2012;Coates et al., 2013), using a parameterization for mesic stands (Supplementary Table 3). We estimated the minor components in this study of cottonwood using trembling aspen parameters, and Douglas-fir and western larch using interior spruce parameters. We initiated runs for each of our 75 stands using stems·ha −1 in 2 cm diameter class increments and two seedling height classes (<30, 30-130 cm) randomly initiated to generate a 4 ha simulated stand. Based on the time since fire at field sampling date, we ran each stand to 100 years since fire and analyzed the trees from a 1 ha area in the center of each simulated stand. We estimated live carbon by applying biomass functions and carbon factor (50% of live biomass) described in Biomass and Carbon Calculations -Trees for individual trees in the simulated runs. We used SORTIE-ND version 7.5.04 (Canham et al., 2005) to run the model and rsortie (Beukema and Clason, 2022) for SORTIE initiation and output processing in R.

Effect of Planting on Carbon
To determine whether tree planting affected carbon storage over time since fire, we fit linear mixed effects model with a normal probability distribution (Bolker, 2008) to five carbon response pools: (1) Total, (2) live aboveground, (3) dead aboveground (woody debris and standing dead trees), (4) forest floor (litter and forest floor), and (5) mineral soil. We included a random intercept (varying among fires) to capture the hierarchical structure in our data collection. All environmental covariates were handled as fixed effects, centred to the mean and scaled prior to analysis to facilitate model fitting and results interpretation. As our study was designed to understand the effect of tree planting, we compared two alternate models for each carbon pool: one with and one without the planting treatment. If including the planting treatment improved the fit of the model so that the AIC dropped by more than 2 units, we considered the planting treatment significant. The remaining covariates (time since fire and all environmental covariates) were included in all models (with or without planting) to account for these processes and conditions that influence forest carbon. Although coarse fragment content was used to calculate mineral soil carbon, it was still included as a covariate for that pool for consistency. We describe the relative 1 http://www.sortie-nd.org/ 2 https://bvcentre.ca/sortie-nd influence of each covariate on each carbon pool by their slope estimates in the best model (with or without planting). We tested the impact of fire severity (dNBR) across carbon pools for each plot from fires after 1985 (N = 53) with the same model design (one model with planting treatment included and one without). We used base R, and fit likelihood models using lme4 (Bates et al., 2015) for carbon pool linear mixed model fitting. We calculated goodness-of-fit (R 2 ) for the overall model, a model with only fixed effects, and a model with only random intercepts using the R package rsq (Zhang, 2021).

Change in Live Carbon Over Time
To determine the change over time since fire in live carbon stores with and without tree planting, we fit linear and parametric nonlinear (Chapman-Richards) functions (Bolker, 2008) to the live aboveground carbon (i.e., excluding roots). We fit these models to simulated stands in SORTIE to capture the long-term (100year fire return interval) carbon implications of stands that were planted versus not-planted. We used the likelihood (Murphy, 2015) package to fit models projected carbon to 100 years since fire as this packages enables flexible function fitting. We verified model convergence and selected the best functional shape based on lowest AIC score.

Regeneration
We tested whether both the mean (t-test) and variance [Levene's test; car package (Fox and Weisberg, 2019)] of seedling density differed between young (<30 years) planted and not-planted stands and used a Wilcoxon rank sum test to determine whether the dominant species (pine, spruce, fir, and aspen) seedlings densities differed between treatments. We present figures of regeneration on a log scale for better visual representation, but regeneration data were not analyzed on a log scale. All data were processed and analyzed in R Core Team (2020), with code available at https://github.com/aclason/Frontiers_ FireRehab.

RESULTS
Total ecosystem carbon was lower in planted stands compared to not-planted stands over 60 years since fire (Table 1 and Figure 2). Of the stands that were 50-60 years old at the time of our study, planted stands contained an average of 142 Mg·ha −1 compared to 147 Mg·ha −1 in not-planted stands. This difference in total ecosystem carbon was largely driven by a significantly greater abundance of dead wood carbon in not-planted stands compared to planted (Table 1 and Figure 2). Including the planting treatment improved the model for all carbon pools ( Table 1). Planting increased live carbon and mineral soil carbon, but decreased dead, forest floor, and total carbon (Table 1 and Figure 2). Time since fire was important in predicting carbon across pools, except for mineral soil (Table 1 and Figure 2). Total, live and forest floor carbon increased, and dead carbon decreased over time since fire (Table 1 and Figure 2).
Some carbon pools were better fit by the treatment and covariates included in our models compared to other carbon pools. The mineral soil and live carbon pools had the best model fit (overall mineral soil R 2 = 0.70, overall live R 2 = 0.75). Fire identity explained little of the variance across most pools as a random effect, with the highest impact on the dead carbon pool (R 2 = 0.11). Overall model fit was lower for total, dead, and forest floor carbon pools (overall R 2 = 0.60, 0.46, and 0.34, respectively), suggesting greater variability, or missing covariates that might capture more of the observed variability for those pools. Including fire severity as a covariate for the years available (N = 53) did improve model fit for total (R 2 = 0.67), dead (R 2 = 0.53) and forest floor (overall R 2 = 0.49) carbon pools (Supplementary Table 4).

Live Aboveground Carbon
Live carbon accumulated in stands over time and comprised half of the total carbon pool of planted stands by 60 years after fire (Figure 2). Including the planting treatment improved the model, indicating that planting increased live carbon, although by a small amount-an estimated 4.1 Mg·ha −1 by 60 years for live aboveground biomass (Tables 1, 2). Young not-planted stands did maintain some live carbon immediately after fire, but it was a relatively small contribution to the total ecosystem carbon (Figure 2). When we simulated these stands forward in time, both linear and non-linear models had lower AICs when planting was included ( Table 3), suggesting that the small but detectable difference between planted and not-planted stands at 60 years would continue over a 100-year fire return interval. We found that the Chapman-Richards function had a better overall fit than a linear function (Figure 3 and Table 3), suggesting that biomass accumulation over time since fire has a slightly sigmoidal shape (Figure 3). While climate variables (growing degree days and mean annual precipitation) were not important predictors of carbon stores across this study, estimating the relationship between carbon pools, particularly live tree carbon, and climate (Table 1) is important for projecting future carbon stores with climate change. Heat load index was an important predictor for live carbon, with the negative slope suggesting less live carbon stored on warm, SW-facing aspects ( Table 1).
Over time, both planted and not-planted stands were dominated by lodgepole pine and interior spruce in the largest diameter classes, but deciduous species (trembling aspen and paper birch) were present throughout the chronosequence in many plots in both planted and not-planted stands (Supplementary Figure 3).

Dead Aboveground Carbon
There was a decrease in dead wood carbon resulting from planting treatments (Tables 1, 2, and Figure 2). When we included all plots across the 60-year chronosequence, there was also a significant decline in dead wood over time (Figure 2 and Table 1). For the 53 plots with severity estimates from dNBR (Guindon et al., 2020), increasing fire severity increased dead wood carbon (Supplementary Table 4), particularly in standing snags (Figure 4B), and time since fire became less important as a predictor of dead wood carbon.
There were significantly more standing snags and downed coarse woody debris in stands that were not planted compared to those that were planted after fire, and these differences remained over time ( Figure 4A and Table 2). The CWD pool was substantially larger in not-planted stands across the 60-year chronosequence, but in particular for the 20-30-year old stands.  Table 1). Fire identity explained most of the variance in the forest floor carbon pool when fire severity was included in the model (R 2 FIGURE 2 | (A) Carbon (Mg·ha −1 ) over time since fire with and without tree planting across five carbon pools: Total, live, dead, forest floor, and mineral soil. (B) Cumulative contribution of four carbon pools (live, dead, forest floor, and mineral soil) to total carbon with and without tree planting, and over time since fire. Data are averaged across stands within planting treatments in 10-year time since fire bins (<11, 11-20, 21-30, 31-40, 51-60; no stands in 41-50 years bin).

Forest Floor Carbon
fire identity = 0.40; Supplementary Table 4 and Figure 3). Other environmental covariates did not have significant effects on forest floor carbon.

Mineral Soil Carbon
The mineral soil pool contained more than half of the total carbon for planted stands less than 15 years post-fire ( Figure 2B).
Mineral soil carbon was higher in planted stands, increased with clay content and the C:N ratio, but decreased with coarse fragment content (Table 1). Although carbon in the mineral soil was higher in planted stands, these stands also had lower carbon in the forest floor by about the same amount (difference in intercepts was 3.9 for mineral and 3.1 for forest floor). Mineral soil carbon did not change over time (Table 1). 41-50 n/a n/a n/a n/a n/a 51-60 59.9 (15.7) 13.3 (3.5) 1.7 (4) 1.9 (3.8) 1.9 (3.1) 100* 140.8 (19.6) n/a n/a n/a n/a 41-50 n/a n/a n/a n/a n/a 51-60 55.8 (21.2) 12.4 (4.7) 7.2 (6) 11.7 (7.3) 8 (6.1) 100* 136.9 (27.5) n/a n/a n/a n/a *Simulated using SORTIE and does not include roots.

Regeneration in Young Stands
Tree regeneration occurred across all sampled postfire planted and not-planted stands (seedlings density > 0; Figure 5). On average young postfire stands that are planted and not-planted have similar tree regeneration densities and variability (Figure 5). There was no difference in the mean seedling density of young stands (<30 years since fire) between not-planted and planted stands (t = −1.58, p = 0.13), nor was there a difference in the variance between these treatments (F = 2.4, p = 0.1). There was one young stand (6 years post fire) with no conifer regeneration, but deciduous regeneration was present. Regeneration throughout the chronosequence was dominated by lodgepole pine, but with a diversity of species in both planted and not-planted stands. The dominant species did not differ in total seedling density between planted and not-planted stands (lodgepole pine W = 873.5, p = 0.06; interior spruce W = 757.5, p = 0.55, subalpine fir W = 798.5, p = 0.21, and trembling aspen W = 668, p = 0.45; Figure 6). Western larch (Larix lariana [Nutt.]) and Douglas-fir (Pseudotsuga menziesii [Mirb.] Franco) regeneration were only present in the planted stands (Figure 6). After 40 years, stands that were not planted had The bold values indicate the best models. The non-linear Chapman-Richards function best fit the data.
ongoing lodgepole pine, interior spruce, and subalpine fir in the regeneration layer, while the planted stands had mostly interior spruce and subalpine fir in this layer (Figure 6).

DISCUSSION
Wildfire leaves a heterogeneous legacy of forest structure and forest floor depth across landscapes which reflects variation in fire severity and affects post-fire recovery (Turner et al., 1997;Burton et al., 2008). In managed forest landscapes, post-fire recovery can also be affected by pre-fire management histories of stands and the legacy of post-fire salvage or rehabilitation activities (Marcolin et al., 2019;Povak et al., 2020). We found that across planted and not planted stands that varied in fire severity and time since fire, planting trees did not increase total ecosystem carbon after fire. Site preparation activities associated with planting removed dead wood and there were insufficient gains in the rate of live tree carbon sequestration from planted trees to overcome the loss of the dead wood carbon. While globally only 8% of carbon is found in dead wood (Pan et al., 2011), the majority of biomass remaining after fire is contained in this carbon pool (Gough et al., 2007;Donato et al., 2013;Stenzel et al., 2019;Lutz et al., 2020). We found the largest difference in carbon pools between planted and not-planted stands was the higher dead wood carbon, with significantly higher carbon in dead standing snags and downed wood stored in stands that were not planted. To reduce slash and eliminate overhead hazards for tree planters, post-fire management practices in BC have included site preparations such as knocking trees down, pile and burning dead wood, or salvage harvesting. Soil disturbances such as mounding, or trenching can also be used to prepare the site for planting. Whereas underplanting the standing snags would maintain the important dead carbon pool in planted stands, we did not find these practices common in our study area (only four stands were underplanted). Continued use of knockdown, pile and burning and soil preparation treatments in the central interior following more recent fires in 2018, suggests an ongoing need to explore options that leave dead wood on site as legacy carbon storage and provide safe working conditions for tree planting (i.e., North et al., 2019). Maintaining the dead wood carbon pool may conflict with other forest objectives, such as managing fuels for future wildfire risk Lutz et al., 2020). Understanding the tradeoff between carbon and fuels will be important in determining the desired postfire management practices across different portions of managed forest landscapes.
Snags constitute a large proportion of the dead biomass following fire (Campbell et al., 2007;Donato et al., 2016;Stenzel et al., 2019). Dead standing trees store carbon after fire (Stenzel et al., 2019), often decaying at a slower rate than downed wood due to lower moisture and decreased soil-borne decomposer activity (Angers et al., 2011;Kaytor, 2016). Planted stands had fewer snags, which is not surprising because site treatments or salvage harvesting would have removed or knocked down standing trees. We found that CWD contributed more carbon than snags to the dead carbon pool both in planted and not-planted stands over time which may be the result of prefire tree mortality. Our study occurs across a landscape heavily affected by a mountain pine beetle outbreak starting in the mid-1990s (Aukema et al., 2006). Fires that burned after the outbreak would have encountered a higher abundance of standing dead trees (Lewis and Thompson, 2011;Rhoades et al., 2020) and a higher abundance of CWD post epidemic (50-90% of MPB killed trees fall 12-16 years after death; Harvey, 1986;Mitchell and Preisler, 1998;Audley et al., 2021). We found higher carbon stored in snags versus downed wood only in not-planted stands that burned at high severity, which supports the finding that increased burn severity can decrease snag fall rates after fire (Angers et al., 2011). As we observed changes over time using a chronosequence and did not monitor stands over time, it is not possible to determine trends in snag fall, CWD recruitment, or CWD decay rates.
Combustion of the woody boles of large trees is low during fire, with more complete combustion of the forest floor across fire severity (Campbell et al., 2007). Consequently, the gradual development of the forest floor is a key process for post-fire carbon recovery. We found that the forest floor carbon pool was lowest shortly after fire, and increased over time, with slightly more carbon in not-planted stands compared to planted. Using our linear model from not-planted stands to predict forest floor carbon content at 190 years would yield 41-44 Mg·ha −1 , which is a reasonable estimate compared to longer-term studies [32-47 Mg·ha −1 at 190 years in Kranabetter (2009) and 35-78 Mg·ha −1 for >140 years in Fredeen et al. (2005) and Kranabetter and Macadam (2007)]. Although a linear response simplifies the recovery process, which is tied to changes in inputs from herbaceous cover to moss cover to tree cover during succession (Deluca and Boisvenue, 2012;Andrieux et al., 2018), others have also successfully fit linear relationships to forest floor carbon increases after fire (Nalder and Wein, 1999 for aspen stands; Andrieux et al., 2018 for black spruce stands). Conceptual models suggest that this carbon accumulation likely reaches an equilibrium and asymptotes with time (Seedre et al., 2011), although an impressive 5,000+ year chronosequence in the Swedish boreal forest suggests that in the absence of disturbance, belowground carbon can continue to accumulate (Wardle et al., 2012). In the model with fire severity as a covariate, in stands 5-35 since fire, the random effect of fire identity explained more variation than any fixed effect, which supports existing models that use the characteristics of individual fires and fire weather indicators such as soil moisture and drought indices to predict forest floor consumption (De Groot et al., 2009). Fire season and intensity, particularly the patterns of surface and crown fire, are also potential drivers of forest floor consumption that we were not able to include in our models. The satellite-derived fire severity index dNBR did not explain forest floor carbon, perhaps because fire effects on tree crowns (which can be remotely sensed by satellite imagery), are not highly correlated with consumption of the forest floor in all stand types (Alonzo et al., 2017).
Post-fire recovery of mineral soil carbon in boreal forests is depicted as a small gradual increase (Seedre et al., 2011), or unchanging (Deluca and Boisvenue, 2012). When soil carbon remains constant, as we observed, regenerating stands create carbon inputs through litter and fine root exudates and turnover as they grow, but these inputs are balanced by losses through decomposition. Small changes in mineral soil carbon would have been difficult to detect across our study area without a greater sampling effort because of stand differences in clay content, coarse fragment content and the C:N ratio that had strong effects on soil carbon (despite our efforts to only sample within a narrow range of site productivity). Site variability in the boreal is a common issue, although carbon increases with time can still be observed (Pregitzer and Euskirchen, 2004), sometimes even in a relatively short period- Seedre et al. (2014) found an increase in mineral soil carbon within 27 years after fire in boreal mixedwoods. Over a 300-year chronosequence in black spruce stands, Andrieux et al. (2018) detected indirect effects of time on pH (which increased metal oxides thereby increasing mineral soil carbon), demonstrating the complexity of the processes connecting belowground carbon with aboveground disturbance. Below and aboveground interactions can also effect carbon sequestration via productivity, but we did not observe an effect of mineral soil properties on the live carbon pool (even though soil C:N ratio can be an important predictor of tree growth; Omari et al., 2021). Although the slightly higher biomass of the planted stands in this study would have produced greater amounts of belowground carbon from root litter and exudates, net primary productivity and soil organic matter accumulation do not have a simple linear relationship (Jackson et al., 2017). If live biomass was a strong driver of mineral soil carbon we would have also expected to see an increase in carbon over time since fire. Instead, the most parsimonious explanation for the higher mineral soil carbon and lower forest floor carbon in planted stands is that site preparation treatments incorporated the forest floor (and possibly other dead carbon) into the mineral soil.
We found that live carbon was higher in planted stands compared to not-planted stands by 60 years since fire. This effect was small [decrease in AIC by 5 with the inclusion of planting treatment, and a difference in intercept (mean) of 7.9 Mg·ha −1 ], but detectable. More consistent and earlier establishment of tree regeneration, coupled with site preparation, and the use of seedlings produced through tree improvement programs, may explain the higher carbon stored in live trees in planted stands. Live tree carbon can also be higher in plantations after harvesting compared to wildfire, not because of increased growth rates, but from more live trees retained during harvest for wildlife habitat, compared to wildfires with complete tree mortality (Seedre et al., 2014). This was not the case in our study, where partial mortality from fire led to more live trees remaining in the not-planted stands. The planted stands had lower live carbon in larger diameter classes shortly after fire likely because of post-fire salvage and site preparation removing the live trees that survived wildfire. We estimate that a linear projection of growth and carbon accumulation in live stems to 100 years since fire is 0.1 Mg·ha −1 year −1 faster in planted stands (1.4 Mg·ha −1 year −1 ) compared to naturally regenerating stands (1.3 Mg·ha −1 year −1 ). The nonlinear model fit suggests an earlier peak in biomass accumulation in planted compared to not-planted stands. Similarly, in southern boreal ecosystems, biomass accumulation rate may peak later in naturally regenerating stands (77 years; Pare and Bergeron, 1995) compared to young managed stands (23-27 years; Repo et al., 2021).
Sub-boreal forests in BC are adapted to frequent stand replacing or mixed severity fire events (DeLong, 1998). These ecosystems support serotinous lodgepole pine that can retain aerial and soil seed banks even after widespread mountain pine beetle outbreaks (Teste et al., 2011) and regenerate successfully after fire even after the outbreak (Talucci et al., 2019). We also find support that these forests regenerate effectively after fire, even in the younger fires that would have occurred after a mountain pine beetle outbreak. We did not find evidence of natural regeneration failures following wildfire, with naturally regenerating stands containing similar densities of trees as planted stands. Recruitment of lodgepole pine can come as a rapid pulse (Charron and Greene, 2002) and/or as ongoing regeneration for upward of 10 years after fire (Stevens-Rumann et al., 2018). Regeneration densities over our 60-year chronosequence suggest that not-planted stands had a longer period of regeneration, but these stands all eventually became stocked. While we did not measure the spatial pattern of regeneration, the longer period of regeneration and lack of tree planting (even spacing) would likely result in different size and within stand spatial structure (Ziegler et al., 2017). This type of temporal and spatial variability in seedling establishment following insect outbreaks and fire, creates landscape-scale heterogeneity that can decrease the future risk of large-scale outbreaks (Seidl et al., 2016).
A common concern of planting trees compared to allowing natural regeneration to establish a forest, is lower tree species diversity in planted stands (Dampier et al., 2007). We found similar tree and seedling species composition in both planted and not-planted stands. In both cases, lodgepole pine dominated the regeneration layer shortly after fire, and a mix of trembling aspen, interior spruce and subalpine fir occurred at lower densities in both treatments. Interior spruce and subalpine fir are more dependent on the fire coinciding with a mast year (Pounden et al., 2014) and distance to seed source than pine, due to the lack of a seed bank (Greene and Johnson, 2000), but still established well after fire in our stands, and continued to establish in the understory of older planted and not-planted stands. While generally composition was similar, planted stands included Douglas-fir and western larch seedlings, as a result of assisted migration. Species distribution models suggest these species will be well suited to future climates with the sub-boreal becoming more suitable for temperate species (MacKenzie and Mahony, 2021). Climate may change more rapidly than species are able to migrate to occupy new suitable climatic habitat (Aitken et al., 2008), so planting trees based on future suitability may assist in creating more resilient forests to climate change (MacKenzie and Mahony, 2021).
Although our results suggest that risk of post-fire establishment failure was low for the past 60 years in central interior BC, risks could change in the future with climate change. Growing season conditions, such as increased drought, as well as more frequent, short-interval, overlapping fire events, can affect tree establishment and forest recovery following fire (Mansuy et al., 2012;Harvey et al., 2016;Hansen et al., 2018;Stevens-Rumann and Morgan, 2019). Summer droughts will impact both natural and planted tree seedling survival (Ouzts et al., 2015), so a greater understanding of the ecological context for the risk of tree establishment failure following wildfire is important in deploying reforestation programs effectively (Stevens-Rumann et al., 2018).
Previous forest cover is an important determinant of post-fire forest composition and carbon recovery in natural regenerating stands . We did not consider pre-fire tree composition in our analysis because of data limitations (variable availability and quality of forest inventory pre-fire) and because pre-fire composition will have less influence on post-fire recovery in stands that are planted. Our study was also challenged by the highly variable and potentially confounded effects of site preparation methods and time. Techniques used in the forest industry to prepare a site for planting have changed over time, with recently increased efforts in fire rehabilitation to leave dead wood on site. Anecdotally and with field observations, site preparations may have been more intense in older planted stands prior to planting compared to younger planted stands. We were not able to reconstruct silviculture histories for all stands from available records, and we also did not sample to capture the impacts of site preparation method. Changes in site preparation techniques may then result in different carbon recovery trajectories in planted stands in the future.

CONCLUSION
Forest carbon is not uniformly or completely combusted during a fire (Campbell et al., 2007;Donato et al., 2013;Stenzel et al., 2019), with as much as 83% of the biomass remaining in temperate forests following fire (Donato et al., 2013). Fire legacies include dead standing and down trees resulting in landscapes that do not start at zero carbon following fire. These legacies are critical to account for in post-fire carbon dynamics. The rehabilitation of recently burned forests through tree planting may be an effective nature-based solution to climate change (Graves et al., 2020), particularly when high disturbance severity or frequency would lead to long delays in live carbon recovery (Fu et al., 2017). However, we found that in BC's central interior, natural regeneration reforested stands effectively after fire, and the removal of dead wood from planting activities resulted in a loss in total ecosystem carbon that was not recovered through faster growing live carbon stores.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://github.com/ aclason/Frontiers_FireRehab.

AUTHOR CONTRIBUTIONS
AC and EL conceived the study. AC designed and ran the study with contributions from EL, contributed to data processing, conducted the analysis with input from EL and IF, and analyzed the simulation model. IF collected the data, conducted lab work, and contributed to data processing. AC, EL, and IF contributed to the manuscript and approved the submitted version.

FUNDING
We gratefully acknowledge the financial support of the Government of Canada through the Low Carbon Economy Leadership Fund and Province of British Columbia through the Ministry of Forests, Lands, Natural Resource Operations and Rural Development, with additional funding contributed by Canada Summer Jobs, and ECO Canada.

ACKNOWLEDGMENTS
Thanks to the project committee which includes John DeGagne and Anne-Marie Roberts at FLNRORD, Phil Burton at UNBC, and Brian Watson and Sara Lazuruk from the Forest Carbon Initiative. Thanks to Lori Daniels at the University of British Columbia for sharing Wildfire Field Plot Protocols, which modify the NFI protocols to support field sampling. Thanks to Sam Coggins for a literature and post-fire forest management practices review. Thanks to John DeGagne, Glen Buhr, Ljiljana Knezevic, Bruce LaHaie, Aaron Benterud, Richard Forget, Jacek Bankowski, Richard Kabzems, and Dawn Stronstad for feedback on postfire tree planting practices in BC. Thanks to Maliya Cassels, Ian Russell, and Jocelyn Biro for assistance in the field. Thanks to Phil Burton, Brian Watson, Anne-Marie Roberts, and John DeGagne for providing context and feedback on study objectives and design. Thanks to Wayne Salewski for insight into the Kenney Dam fire, and Elisabeth and Joe Doerig for accommodation and perspective on living with recent wildfires. Thanks to Dave Coates (ret. FLNRORD) and Sybille Haeussler for insight into the Swiss Fire and post-fire rehabilitation in central BC. We also thank reviewers for comments and suggestions that greatly improved this manuscript.