Impact and recovery of forest cover following wildfire in the Northern Rocky Mountains of the United States

Background Anthropogenic climate change is expected to catalyze forest conversion to grass and shrublands due to more extreme fire behavior and hotter and drier post-fire conditions. However, field surveys in the Northern Rocky Mountains of the United States show robust conifer regeneration on burned sites. This study utilizes a machine learning (GBM) approach to monitor canopy cover systematically on a census of burned areas in two large wilderness areas from 1985 to 2021, to contextualize these recent field surveys and create a monitoring baseline for future change. Results A predictive model was developed from coincident LiDAR and Landsat observations and used to create time series of canopy cover on 352 burned sites (individual wildfires subset by number of times burned), which were then summarized using fire impact and recovery metrics. Fire impact, defined as canopy cover loss relative to pre-fire condition, was highly correlated with burn severity (Spearman’s R = 0.70). Recovery was characterized by the following


Introduction
Fire is a natural force of succession and evolution on western North American landscapes.Prior to European colonization, both individual species and ecosystems were maintained by episodic fire and evolved under the selective pressure of fire (Pausas and Keeley 2009;Stevens et al. 2020).However, the nature of fire is changing under anthropogenic climate change and historic fire suppression.Wildland fuels are more abundant and continuous today than they were prior to European settlement (Hagmann et al. 2021).At the same time, warmer and drier conditions are associated with more area burned (Westerling 2016;Abatzoglou et al. 2021) and higher-severity fire (Parks and Abatzoglou 2020).
In light of shifting climatic conditions, there is a growing concern that some forests will not recover after fire and instead transition to shrub or grass systems.As more area burns at high severity (Parks and Abatzoglou 2020), stand replacing fire becomes more common on the landscape, which is the first step towards forest conversion (Parks et al. 2019).Failure to recruit a new cohort of trees in these areas ultimately leads to conversion (Davis et al. 2020), and small trees are more sensitive to climate extremes than their mature counterparts (Trouillier et al. 2019).Short-term drought (Young et al. 2019), chronically dry conditions (Urza and Sibold 2017;Boucher et al. 2020;Stewart et al. 2021), and temperature extremes (Kemp et al. 2019) all limit seedling recruitment in early seral forests.
In western North America, evidence is emerging of widespread conversion from forest to non-forest, currently and under expected future climate conditions.One third of forests in the Klamath region of northern California and southern Oregon could transition to shrublands by the end of the century (Serra-Diaz et al. 2018).In the boreal forests of Northern Alberta, Canada, wildfire is projected to facilitate the conversion of half of conifer forests to hardwoods or grass under some future climate scenarios (Stralberg et al. 2018).In the southwestern United States, increased fire activity is already associated with less conifer regeneration and greater shrub presence (Keyser et al. 2020).
Forests in the Northern Rocky Mountains of the conterminous United States may be more resistant to conversion than other systems in western North America.Under a 2 °C warming scenario, Montana forests are expected to experience less stand replacing fire and more post-fire recruitment than forests of the central and southern Rocky Mountains (Davis et al. 2020).Additionally, recent field surveys of burn scars in the Northern Rockies indicate abundant conifer regeneration and survival when close to a seed source (Clark-Wolf et al. 2022).Similarly, Jaffe et al. (2023) documented abundant conifer regeneration across all fire histories in mixed conifer forests of the Selway-Bitterroot wilderness.Further west, Povak et al. (2020) also found significant regeneration following wildfire in mixed conifer forests of the Okanogan Highlands and Eastern Cascades, which they partially attribute to mild postfire weather conditions that may have reduced water stress on seedlings.Collectively, these findings indicate healthy forest recovery after wildfire in the region, at least in mixed-conifer forests.However, by mid-century, large swaths of the Northern Rockies are expected to experience summer heat extremes too warm to support widespread conifer seedling regeneration (Kemp et al. 2019).Northern Rockies forests currently sit at a crossroads: current conifer recruitment appears optimistic, but widespread conversion and community reorganization is possible before the end of the century.
Remote sensing is one of a few tools available to monitor forest recovery systematically across large land areas.In particular, Landsat has been used extensively to study land surface change.Since the entire Landsat collection became free and publicly available in 2008, multiple metrics of analysis and methods for deriving ecological information from time series have been developed (Zhu 2017).Indices such as normalized difference vegetation index (NDVI) (Kriegler et al. 1969) and normalized burn ratio (NBR) (García and Caselles 1991) are spectrally derived metrics that are commonly proxies for vegetation through time.These metrics are simple, easy to calculate, and correlate with a wide range of vegetative characteristics in multiple biomes (Huang et al. 2021).However, the proliferation of cheap and accessible computing power and the development of new machine learning techniques have reduced reliance on spectral proxies and allowed for more direct modeling of forest attributes such as cover (Matasci et al. 2018a), height (Hudak et al. 2002;Pascual et al. 2010), and biomass (Pflugmacher et al. 2014;Sun et al. 2022) from satellite imagery.
Additionally, there has been extensive development in methodology for deriving disturbance information from a time series of metrics (Zhu 2017).LandTrendr (Kennedy et al. 2010) uses a segmentation strategy to divide time series into components meant to reflect periods of landscape stasis and change and has been used extensively for diverse monitoring tasks since its development (Bright et al. 2019;de Jong et al. 2021;Runge et al. 2022).It has been used to detect disturbance when the locations are unknown (Cohen et al. 2018) and has been designed for sensitivity to both abrupt and subtle change (Kennedy et al. 2010).Similarly, the Vegetation Change Tracker (VTR) is designed to track abrupt change utilizing thresholding (Huang et al. 2010).However, both approaches are univariate and unable to accommodate multiple spectral indices in a single analysis.Recent innovations in remote sensing utilize machine learning to aggregate multiple spectrally derived metrics (Moran et al. 2020;D'Este et al. 2021;Sun et al. 2022).Machine learning approaches can also integrate information from other sensors such as LiDAR.
Aerial LiDAR is well-suited for quantifying forest structure because it has capacity for fine spatial resolution and is the best remote sensing tool available for measuring in the vertical dimension (Hyde et al. 2006).The comparative advantage of LiDAR in machine learning is the massive training datasets it avails along with systematic and consistent measurements of forest structural attributes.However, LiDAR alone is limited as a monitoring tool because there is non-continuous spatial coverage and no temporal consistency compared to passive, satellite-based remote sensing systems such as Landsat and Sentinel.
LiDAR-Landsat covariance is the process of identifying relationships between coincident Landsat imagery and LiDAR observations to overcome the limitations of each respective sensor.This method can be used to extrapolate LiDAR-derived metrics beyond their original geographic footprints (Hudak et al. 2002;Wilkes et al. 2015;Matasci et al. 2018a) or to update older LiDAR data without collecting another acquisition (Matasci et al. 2018b).
This study uses Lidar-Landsat covariance to produce a systematic assessment of fire impact and forest recovery on a census of large, burned areas in the Northern Rocky Mountains.The purpose of the remote sensing approach is to contextualize field-based findings indicating widespread conifer regeneration with a growing body of literature describing the potential for future forest conversion (Stevens-Rumann et al. 2018;Stralberg et al. 2018;Serra-Diaz et al. 2018;Keyser et al. 2020).The units of analysis for this study are wildfire events subset by the number of times they experienced reburning.Previous pixel-scale research has focused on segmenting time series into ecologically interpretable units (Kennedy et al. 2010) or assigning categorical values to pixels corresponding to a modeled shapes (Moisen et al. 2016).Our study uses simpler, more direct metrics to summarizing fire impact and recovery at the scale of individual wildfires with known locations and dates.The value of this approach is to reduce non-linear time series into intuitive, reproducible metrics that can be compared across space and time at the scale of individual burned areas.Evaluating fire impact and recovery at event scale also increases relevance to management as post-fire remediation and restoration efforts are usually funded and implemented on a fire-by-fire basis.
Our impact and recovery metrics were derived from time series of modeled annual canopy cover.A gradient boosted regression model (GBM) was developed to predict canopy cover from five spectral indices and three topographic variables.The model was trained on forest cover measurements derived from ten LiDAR datasets encompassing 160,000 ha.Annual cover estimates were then produced for the period 1985-2021 at 30 m spatial resolution across the Selway-Bitterroot and Bob Marshall wilderness complexes, to create 352 individual 35-year time series depicting the progression of canopy cover on a census of burned sites, where each time series corresponds to an area with a unique fire history.Fire impact on canopy cover and recovery rate were then calculated from each time series.
Fire impact was defined as the difference between mean pre-fire canopy cover and the absolute minimum post-fire canopy cover, which could occur at any point after a fire event to account for secondary mortality.Sites were considered to be recovering if they were gaining canopy cover post-fire, as evidenced by a positive slope of a linear regression fit to annual canopy cover predictions occurring after the observed absolute minimum.On recovering sites, a recovery rate was also calculated by using the regression slope to project how many years it would take to recover to each site's pre-fire canopy cover.Sites were considered non-recovering if they never began gaining canopy cover post-event or if they had projected recovery times greater than 500 years.Because the prognosis for conifer regeneration is often correlated with fire severity (Davis et al. 2020) and fire history (number of times burned), impact and recovery were also evaluated as a function of area burned at high severity and number of times burned during the observation period.
The objectives of this study are to: 1. Monitor canopy cover on a census of 352 burned wilderness sites from 1985-2021 using machine learning to integrate forest observations from multiple sensors 2. Quantify fire impact and recovery of individual burned areas using direct, reproducible metrics 3. Determine if impact and recovery correlate with fire severity and number of times burned 4. Assess whether length of record influences the probability of observing recovery

Study area
The study area includes two adjacent wilderness areas (Fig. 1) representing a vast swath of relatively intact ecosystems with similar vegetation and fire regime characteristics and covering 1.13 million ha (Teske et al. 2012).
It is topographically complex, with elevations ranging from 530 to 3200 m.Precipitation ranges from 63 cm in the south to over 178 cm (Teske et al. 2012).Plant communities are diverse, with Pinus ponderosa (Ponderosa pine) occurring at low elevations, transitioning to Pseudotsuga menziesii (Douglas fir) and Larix occidentalis (western larch) on cooler sites (Parks et al. 2015).
Remnant patches of the pacific maritime community are widespread, including Tsuga heterophylla (western hemlock) and Thuja plicata (western red cedar) (Rollins et al. 2002;Larson et al. 2020).Abies lasiocarpa (subalpine fir) and Pinus contorta (lodgepole pine) are common at higher elevations (Rollins et al. 2002;Teske et al. 2012;Larson et al. 2020).These forests contain diverse fire regimes, characterized by both mixed severity and infrequent stand replacing fire (Teske et al. 2012), with observed fire rotations ranging from 35 to 350 years (Rollins et al. 2002).Total area burned area is primarily driven by large wildfires, though the vast majority of fires are small.Teske et al. (2012) estimate that 71-75% of fires in the BMC and SB are less than 405 ha.High-severity fire is common on the landscape, with roughly 35-40% of the land area burning in high severity events per century (Morgan et al. 2017;Teske et al 2012).Reburned patches are common but they tend to be spatially small, as existing burn scars tend to limit fire spread (Teske et al. 2012).

Fire history and units of analysis
A census of 352 individual sites (total area of 4523 km 2 ) were included in the analysis (Fig. 2).Of the 352 sites studied, 175 burned once, representing 89% of the studied land area.One hundred forty sites burned twice (10.5% by area), and 37 sites burned three or four times (0.38% by area).Portions of fires that reburned in subsequent incidents were treated as unique sites.For example, the Footstool Fire burned in the central Selway in 1988.A portion then reburned in the 2016 Moose Fire.The once-and twice-burned portions of the Footstool Fire are separate units of analysis for this study.

Response variable: LiDAR-derived canopy cover
Gradient boosted machine learning (Hastie et al. 2001) was used to predict LiDAR-derived canopy cover from Landsat and topographic inputs.Canopy cover was defined conventionally as the proportion of LiDAR returns where Z is greater than 2 m in height (Peterson et al. 2015;Moran et al. 2020).The model was trained on spatially coincident Landsat and LiDAR data from 10 sites (160,000 ha) across Montana and Idaho (Fig. 1), representing a wide range of forest structures and biomes present in the study area.Sites ranged from 4681 to 36,299 ha in area.Each of the ten sites was collected as part of different LiDAR acquisition initiatives over the course of a decade, so sensor and collection date varied.However, all LiDAR was collected during the summer months between 2009 and 2020.Each LiDAR acquisition had at minimum 10 cm vertical accuracy, more than 8 returns per m 2 , and nominal pulse spacing < 0.35 m.Canopy cover was calculated at 30 m using a Landsat raster as a template to avoid resampling the LiDAR once processed.Machine learning was performed using the R package H2O (2021).Raster processing and data cleaning and visualization were performed with packages raster (Hijmans 2022), terra (Hijmans 2023), dplyr (Wickham et al. 2022), and ggplot2 (Wickham et al. 2023).

Predictor variables: LANDSAT spectral indices and topographic metrics
Landsat and topographic metrics served as the predictor variables in model training and were also used to predict beyond the spatial and temporal extent of existing LiDAR acquisitions.Landsat Tier 1 Collection 1 images were acquired for the flight year of each LiDAR acquisition in Google Earth Engine (Gorelick et al. 2017).Pixels with a medium or high confidence of clouds, shadow, or snow and ice were removed.Median and maximum composite images were calculated for the May-October growing season of each year.Each image has 5 bands: normalized difference vegetation index (NDVI), normalized burn ratio (NBR), and Tasseled Cap brightness, greenness, and wetness, for a total of 10 bands per year.The LANDFIRE Existing Vegetation Type (EVT) (US Geological Survey and Department of Agriculture 2016) dataset was used to query out bare ground, urban development, and irrigated cropland.Pixels were removed for every year of analysis if they fell into one of these categories in 2016.
The GBM model was also trained on elevation, slope, aspect, and topographic heat load index.These factors were included to capture a range of geophysical site characteristics that control, in part, the distribution of canopy cover on the landscape.Elevation, slope, and aspect were taken from LANDFIRE (US Department of Interior and US Department of Agriculture 2022).Aspect was then used to calculate topographic heat load (McCune and Keon 2002) (Table 1).
The 35-year study period spanned three generations of Landsat sensors, incorporating imagery from Landsat 5, 7, and 8 at various points.Landsat 8 was used from 2013 to 2021, Landsat 5 was used from 1985 to 2012, and Landsat 5-7 composites were used in 2012 to offset the impact of the scan line corrector failure on Landsat 7 (Markham et al. 2004).Linear regression was utilized to transform the Landsat 8 Operational Land Imager (OLI) sensor raw spectral data to match Landsat 5 and 7 Enhanced Thermal Mapper (ETM+) prior to machine learning (Roy et al. 2016).The GBM model was also given a binary variable identifying the source sensor of each pixel to internally correct for subtle differences between the ETM+ and OLI sensors.

GBM modeling
GBM is a supervised machine learning method that iterates through a series of regression or classification trees and improves as it builds more trees.It is an ensemble method in the sense that final predictions are drawn from a composite of multiple regression trees (Hastie et al. 2001).GBM was chosen for this study because a previous study indicated that a single, general GBM model can perform well on varied landscapes (Moran et al. 2020).
Two grid searches were used for model tuning.First, a preliminary grid tested several tree depths while holding all other hyperparameters constant.This grid search was performed first because initial experimentation indicated that tree depth had the single biggest impact on model performance of any user-defined control.The number of trees in the ensemble was not limited to encourage the machine to build a shallow forest of many trees as a method to prevent overfitting (Elith et al. 2008) and was allowed to sequentially add trees until constrained by a stopping function.Then, once the optimal tree depth was defined, a second grid search was performed across learning rate (0.05 or 0.1), sample rate (testing values between 0.1 and 1 on 0.1 intervals), and column sample rate (values between 0.3 and 1 on 0.1 intervals) to find the best performing overall model.The two-step approach was chosen because it is a systematic, computationally efficient hyperparameter search.An 80/20 training/testing split was used during model development, independent of the 5% validation holdback.Ultimately, after all quality control and validation holdbacks, the model was trained on 4,005,816 samples, with at least 83,913 samples from each training landscape.

Sampling and validation
GBM models improve along a gradient relative to a loss function.In this study, model gradients were calculated relative to root mean square error.Unbalanced datasets are difficult to evaluate using this metric (Branco et al. 2019) so each site was up-sampled or down-sampled to the mean number of samples across all ten sites.Then, the sites were combined and shuffled.Five percent of the data was randomly held back for validation.Samples spaced at distances less than the estimated variogram ranges at each site were removed to address potential bias in the model performance evaluation due to spatial autocorrelation, leaving a count of 37,777 samples in the final validation dataset.

Predicting on burned landscapes and generating time series
The validated model was used to predict canopy cover of each burned pixel in the study area from 1985 to 2021 and then aggregated to each unit of analysis.Fire perimeters were acquired from MTBS.Portions of fires that fell outside the bounds of the study area were removed.
Intersecting polygons from multiple fires produced individual burned units that were analyzed separately.These areas were identified and separated using a spatial union between all fire polygons in the study area and categorized as once-burned, twice-burned, and > 3 × burned.
A time series of median canopy cover per year was calculated for each site.Time was normalized to the time of the incident, with year zero being the fire year.Units that burned more than once were normalized to their first fire.Fires were removed from analysis if the most recent burn occurred after 2016 to give adequate time to track recovery.Fires that occurred in 1985 were also removed so that only fires with a known pre-fire canopy cover were included.

Defining impact and recovery
Pre-fire mean canopy cover was calculated for each time series-providing the reference value for estimation of impact.The absolute minimum canopy cover (blue vertical dashed line, Fig. 3) was calculated from a three-year moving window to minimize the effects of atmospheric variations in satellite imagery year to year.
Impact was defined as the proportion of pre-fire canopy cover lost post-fire.It includes both primary loss from combustion (year 1 post-fire) and secondary mortality (> year 1 post-fire).The mathematical definition of impact was calculated as follows: Two metrics were used to assess recovery: (1) whether a burned area began gaining canopy cover (recovery) during the study period, as evidenced by a positive slope of a linear regression fit to annual canopy cover predictions occurring after the observed absolute minimum, and (2) how long it would take to reach prefire cover given observed gain rates (projected recovery time).Rates of gain were calculated by fitting a linear regression model between canopy cover and years since fire, beginning after the observed absolute minimum cover and including a minimum of 3 years post-minima.Projected recovery time was then defined as years to reach pre-fire canopy cover at the observed rate of change, plus the number of years of secondary mortality before the absolute minimum was reached.Sites were considered non-recovering if they never began gaining canopy cover over the course of the study or if their projected time to recovery was greater than 500 years.The methodology for projecting years to recovery is shown in Eq. 2. Recovery rate was derived by  1988, 2001, and 2007.This time series depicts the once-burned unit.The pre-fire canopy cover was 63%.Fifty-nine percent of pre-fire canopy cover was lost during the fire (impact).Absolute minimum canopy cover was 25.8%, which was observed 4 years post-event.After this point, the site began recovering (gaining canopy cover).It is expected to take 76 years to recover to its pre-fire canopy cover fitting a linear model to the portion of the time series that occurs after the absolute minimum canopy cover.
Nine burned areas increased in cover immediately following fire.Collectively, these nine sites made up 0.8% of the studied land area.Fire impacts on these latter sites were restricted to expanses of low cover, south-facing slopes that greened in the year following fire.These sites are likely non-forested based on their locations and response to fire and were not included for further analysis.

Correlation with severity and fire history
The proportion of area burned at high severity for each unit was calculated to assess whether severity correlates with impact or recovery.Severity was defined using the US Geological Survey Monitoring Trends in Burn Severity (MTBS) dataset (USDA Forest Service and US Geological Survey 2017).Severity was defined as the proportion of a unit's area that burned at high severity at any point in the study period and was calculated cumulatively on sites that burned multiple times.Thirty percent of the study area burned at high severity at least once during the study period.Fire history was considered by comparing burned area impact and recovery among once-burned, twice-burned, and ≥ 3x-burned landscapes.

Statistical testing
Because ignitions occurred across the 36-year time series, each site had a variable number of post-fire observations available for assessing recovery.Two statistical tests were performed to determine if length of observation post-event controls recovery.First, time series were truncated to 5, 10, and 15 years post-event to standardize time series length regardless of year of ignition.Then, a binomial GLM model estimated the probability of a site recovering based on year of ignition and length of observation post-event (5, 10, 15 years, plus entire available record).Second, a series of analysis of variance (ANOVA) tests were used to compare projected times to recover between ignition years when normalized to 5, 10, and 15years post-event.

Model performance
The best performing GBM model produced an overall root mean square error (RMSE) of 9.70% and an R 2 of 0.905.Performance varied across the ten training sites, with R 2 values ranging from 0.75 to 0.99 and RMSE values ranging from 1.33-13.1.The site with the R 2 value of 0.99 is a low canopy cover, short-statured ponderosa pine (2) Years to Recovery = Impact Recovery Rate +Years between fire event and absolute minimum CC and juniper woodland in northeastern Montana (Phillips County).The most important variables for predicting canopy cover were median wetness, median brightness, and median NBR.

Fire impact
Fires burned across a full range of pre-fire canopy cover but were concentrated in areas of 50-70% CC (Fig. 4).Burned areas lost a median 55% of pre-fire cover with loss occurring both immediately and for several years post incident.On many sites, these delayed losses were greater than the initial impact.The median proportion of impact that occurred after the first-year post-fire was 43%, meaning that secondary loss accounted for roughly half of total cover loss for half of the burned areas.The median time to begin recovering was 9 years, with 25% of sites starting in 4 years, 75% in 16 years, and 90% in 20 years.The longest period observed before recovery started was 31 years, on two sites in the SB.Both sites burned in separate 1986 incidents and then again in the 2015 Meeker Fire on a high elevation ridgeline north of the Selway River.Although high-severity fire made up a relatively small proportion of burned area (median: 14%; IQR 4-35%), it was highly correlated with fire impact (Spearman's R s = 0.70) (Fig. 4) and much less correlated with pre-fire mean canopy cover and subsequent years to begin recovering.Twice-and thrice-burned sites were associated with greater impact than their once-burned counterparts.However, the variability within fire histories is greater than the difference in impact observed in once-, twice-, and thrice-burned sites (Table 2).

Fire recovery
Sixty-two burned areas (14% of landscape by area) do not show a recovery signal.281 appear to be recovering, representing 85% of the studied area.Initial canopy cover, fire severity, and fire impact were similar between recovering and non-recovering sites (Table 1).Time was the primary delineating variable, with more recently burned areas less likely to show recovery.Of the sixty-two sites that are failing to recover, thirty-seven occurred since 2003 and 28 occurred since 2010.

Probability of recovery by year
The probability of a site recovering was calculated for each ignition year to evaluate if the length of observation period impacts the probability of recovery.The probability of a burned area gaining canopy cover (e.g., recovering) approaches one when given 30 years of recovery and declines to 0.5 in the last year of observation (2016) (Fig. 5, gray ribbon).However, when fires are normalized to record length by evaluating only the first 5, 10, and 15 years post-event (Fig. 5, colored ribbons), it becomes evident that the probability of recovery is similar for modern fires and fires occurring early in the record.The mid-part of the record (1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003) experienced longer periods of time between fire year and initiation of recovery regardless of length of observation, highlighting an element of temporal variability in recovery.However, once these fires began recovering, the total time to regain pre-fire cover was not different than other sites that burned earlier and later in the period (ANOVA, P < 0.05).
Table 2 Characteristics of recovering and non-recovering sites.Impact is defined as the proportion of pre-fire canopy cover lost during and after a fire.Pre-fire canopy cover, area burned at high severity, and impact were all similar between recovering and nonrecovering sites.Non-recovering sites generally burned more recently

Projected recovery time by year
Median recovery time was 40 years and was stable through the period of record (Fig. 6).No difference was observed in median recovery times between ignition years across the period of record (ANOVA, P < 0.05).
Even when normalized to 5, 10, and 15 years of observation post-fire, there was no difference in time to recover between fire years (ANOVA, P < 0.05), indicating that rates of change are stable regardless of how many years of observation are available.In sum, the year a fire occurred is not related to median recovery time.Ignition year, however, is related to both the probability that a fire will begin recovery and the length of time between the fire year and initiation of recovery (Fig. 4).

Patterns in recovery time
Seven burned areas have projected recovery times longer than 200 years (Fig. 6).Five of those fires experienced a single burn during the study period and two burned twice.The single longest recovery time observed was 365 years, on a 14-km 2 site in the South Selway that burned once in the 2012 Green Mountain Fire.The areas with projected recovery times greater than 200 years appear to be burning later in the record compared to their faster recovering counterparts.One of these slow recovering sites burned in 1988, while the remaining six experienced their first ignition after 2006.Fire impact on these outlier sites was also greater than the study area at large, with a median impact of 70% compared to 55% in the broader study area.

Additional elements of complexity
Repeated fire during the 35-year period of record impacted a modest proportion of the landscape, with 102 sites (13.5% of burned area) burning twice, and twentyeight sites (0.5% of area) burning three or four times.Time series were normalized to the first fire, so repeated fire had a cumulative effect on impact.Once-burned sites lost a median 51% pre-fire canopy, twice-burned lost 59% and thrice-burned lost 69%.Repeated fire was associated with non-recovery, though the cumulative effect of additional fire was small.Seventeen percent of single burned sites failed to begin recovery, compared to 18% of twiceburned and 21% of thrice-burned sites.On recovering sites, the median time to reach pre-fire canopy cover was 37 years on once-burned sites, 43 on twice-burned, and 47 on thrice-burned sites.
Median fire impact was greater on the BMC (64%) than the SBW (50%).However, the SBW experienced more variation in fire impact, with an interquartile range of 38% compared with 16% on the BMC.This is consistent with the more varied forest communities, fire regimes, and terrain on the SBW.The BMC also experienced more widespread recovery, with 89% of the burned area showing recovery compared to 79% of the SBW.However, once sites began recovering, recovery took longer on the BMC at 49 years, compared to 37 years on the SBW.

Discussion
Canopy cover change following wildfire is a complex phenomenon.We found that cover continues to decline for a period of time after fire, perhaps due to secondary mortality and falling snags.At the same time, surviving small trees release into canopy caps and intermediate and dominant stems expand their crowns in response to better access to light, water, and nutrients.In most cases, a new cohort eventually emerges, grows above 2 m in height, and contributes to cover increases.Given that burned forests are highly heterogeneous, all of the aforementioned processes are happening together at variable intensities and time scales across large burned areas.It is therefore difficult to ascribe the specific ecological causes of cover change consistently from one fire to another.For this reason, we caution against interpreting increases in cover after wildfire as the result of regenerating trees exclusively and we acknowledge that our definitions of impact and recovery are simplistic relative to the inherent ecological complexity of these systems.
The 2-m height threshold used to define cover is necessary to avoid over-estimating due to debris, shrubs, rocks, and other terrain features.While the LiDAR based response variable necessitated a 2-m threshold, the underlying Landsat indices of NDVI, NBR, and Tasseled Cap brightness, greenness, and wetness likely detect regenerating conifers before they reach 2 m in height.Thus, we do not know at what stage in early tree growth this model becomes sensitive to regeneration, particularly in the context of other simultaneous processes such as mature canopies expanding into gaps.Additionally, we are not able to distinguish the vegetation types contributing to cover, and it is possible that some of the early recovery detected is due to emerging shrubs such as snowberry, ceanothus, and serviceberry.Our findings are contextualized by post-burn field surveys indicating abundant tree regeneration of primarily western larch, ponderosa pine, and Douglas fir in the Bob Marshall Complex (Berkey 2020) grand fir and Douglas fir in the Selway (Jaffe et al. 2023) and Douglas fir and lodgepole pine in the non-wilderness forest between the two study areas (Wolf 2021).However, the inability to distinguish between short trees and tall shrubs is a limitation of our approach.
Burned sites were used as the unit of analysis in this study because they are the common unit for assessing, planning, and funding forest restoration.However, this approach does not capture intra-fire variability in recovery, which is of particular importance for evaluating fauna habitat (Sitters et al. 2016;Parkins et al. 2018).A pixel-based approach would enable a different analysis more focused on the biophysical factors related to impact and recovery such as elevation, temperature, water availability, and distance to a seed source.Such an approach would be a worthy follow-on investigation.

Impact
Impact was defined as the proportion of pre-fire canopy cover lost during and after an incident (Fig. 3).Our estimates found a median 55% impact across sites, with an interquartile range of 23% to 87%.Area burned at high severity was positively associated with fire impact.This is consistent with previous work showing that high-severity fire increases tree mortality (Hood et al. 2018), and may be of some management concern, as area burned under high severity is expected to increase under future climate scenarios (Parks and Abatzoglou 2020).
Many sites continued to lose canopy cover years after a fire.Roughly half of sites experienced more impact from secondary effects than immediate combustion, with a median 43% of impact occurring after year 2. Secondary mortality has not previously been quantified in the region.However, this finding is consistent with surveys of ponderosa pine in the American southwest that similarly found secondary mortality to be almost equal to initial impact (Stoddard et al. 2018).While immediate mortality is typically caused by thermal injury (Hood et al. 2018), the mechanisms of secondary mortality are less understood and may be related to pre-fire competition, drought stress, and post-fire conditions (van Mantgem et al. 2018).
Repeated fire was uncommon, with 14% of the studied area burning more than once.While there was a positive association between number of reburns and cumulative impact to cover, the additive contribution of each subsequent fire was modest, with an 8% difference in median impact between once-and twice-burned sites and 10% difference between twice-and thrice-burned sites.This impact also translated to present but small differences in recovery probability.Nine of the ten burned areas ≥ 100 ha that are not recovering or recovering slowly from a fire early in the period of record are multi-fire landscapes.The one exception is a part of the Helen Creek Fire in the northwest Bob Marshall which burned in 2000 in a rocky, high-elevation, subalpine forest that is expected to recover slowly.
Because the study period covers only 35 years, all observed fire-on-fire interactions occurred on relatively short intervals.Previous work has defined short interval fire in similar systems as 15 years in Glacier National Park (Hoecker and Turner 2022), 12 years in the central Cascades (Busby et al. 2020), and 30 years in Yellowstone National Park (Turner et al. 2019).Short interval fire is rare on this landscape (Rogeau and Armstrong 2017).Previous work on the SBW, BMC, and neighboring landscapes indicates that fire is generally self-limiting (Parks et al. 2014), and reburned patches tend to be small (Teske et al. 2012).
Parts of the study area likely experienced frequent, lowseverity fire prior to European colonization (Brown et al. 1994) and could reasonably be expected to be resilient to cover loss under repeated fire.Our findings indicate that reburned sites have only slightly less recovery potential than their once-burned counterparts, indicating some resilience to repeated fire.However, short interval fire can also be associated with departure from historic plant communities in the Greater Yellowstone area (Hoecker and Turner 2022) particularly in fire intolerant woody species (Enright et al. 2015) and species dependent on serotinous reproduction (Buma et al. 2013;Turner et al. 2019).Additionally, climate-fire models with strong feedback between fire and subsequent fuel constraints still indicate an increase in area burned between 2021 and 2050 (Abatzoglou et al. 2021), suggesting that the self-moderating effect of fire on the landscape might be dampened and reburning might soon affect more area.

Probability of recovery
Most sites, when given enough time, begin recovering canopy cover.Field surveys of the area have observed robust regeneration (Jaffe et al. 2023), particularly when near a seed source (Clark-Wolf et al. 2022).Our findings are consistent with these studies.When sites are observed for 5, 10, and 15 years post event, the probability of recovery is similar for most years of ignition.The exception is the period from roughly 1994-2003, during which the probability of recovery is dampened in the first 15 years post-fire compared to earlier and later years in the record.
The region experienced multiple exceptionally dry periods during the late 1990s into the 2000s, including an extended drought from Fall of 1999 to Spring of 2005 (NOAA National Centers for Environmental Information 2023).Multiple studies indicate that hot and dry conditions following wildfire tend to suppress conifer regeneration, though the actual timespans studied range from just the year of germination (Davis et al. 2019), 3 years (Harvey et al. 2016;Stevens-Rumann et al. 2018), or 5 years (Urza and Sibold 2017;Stewart et al. 2021).This drought could explain the observed decreased probability of recovery in the first 15 years after fire.However, despite their initial lag, these sites eventually began recovering at rates indistinguishable from those that burned in more normal moisture conditions.Although this finding indicates that forests in the region may be resilient to episodic drought despite initially delayed recovery, it suggests a sensitivity of forest recovery to climate.It also indicates that longer monitoring may be required to fully evaluate conifer regeneration, particularly if sites have only been observed in periods of drought.
Sixty-two sites (14% of studied land area) did not show a recovery signal during the period of record.These sites represent a relatively small portion of the landscape but may be of critical management importance as early harbingers of forest conversion.Most of these non-recovering sites burned recently, particularly from 2011 onward.Our analysis shows a high likelihood many of these recently burned sites will recover with time.However, the Northern Rocky Mountains have trended hotter and drier since the beginning of the study period (Hansen et al. 2006;Hornbach et al. 2016), with a notable increase in maximum summer temperatures (Joyce et al. 2018).While the widespread recovery from fires in 1980s through 2000s is reassuring, the trend has potential to break down as postfire conditions diverge from historic normals.

Time to recover
Median time to recover canopy cover was 45 years.Once sites begin the recovery process, there was no difference in projected time to recover when observing sites for 5, 10, and 15 years post incident, indicating that most sites are stable in their rate of recovery.Median recovery time was also stable through time (Fig. 6), but outliers became more common in recent years.Six fires occurring after 2006 had projected recovery times greater than 200 years, compared to one fire pre-2006.
Parts of the study area have fire return intervals up to 350 years (Teske et al. 2012).This study observed the landscape for 35 years, and almost certainly captured more fire in systems that burn frequently and have evolved resistance to impact (Stevens-Rumann et al. 2018).It is possible that the long recovery outliers in this study represent stand replacing fire in a climax forest community that will naturally take a long time to recover but still exist within the bounds of the historic normal.
Additionally, this study utilizes pre-fire state as the baseline to which recovery is calculated.While wilderness areas have generally been allowed to burn unhindered since the late 1960s, they still experienced fire suppression and subsequent fuel buildup in the decades prior (Berkey 2020).Thus, pre-fire state on these six outlier sites could potentially be an overstocked relic of historic fire suppression, artificially inflating the time to recovery.However, their clustering in recent years raises suspicion that a biophysical or climatic factor may be slowing recovery.Further work is needed to monitor these six sites going forward and contextualize them within their ecological settings relative to factors such as distance to seed source, patch severity, and water availability.

Conclusion
This research builds on a legacy of remote sensing for studying forest impact and recovery from disturbance (Kennedy et al. 2010;Moisen et al. 2016).Time to recover from disturbance, expressed in years, has not been examined in the Northern Rocky Mountains of the US, although similar spectral-based approaches have been used to monitor boreal forests north of the study area (Schroeder et al. 2011;Bolton et al. 2015;Pickell et al. 2016).The work's novelty is in the integration of multiple sensors to inform a growing debate around the prognosis of forest recovery.It harnesses a training set of four million pixels to produce consistent cover estimates with known accuracy, and it can be updated easily to monitor into the future.Designated wilderness is well-suited for this investigation because its unique management context provides a natural laboratory for observing fire and vegetation dynamics relatively unfettered by human influences such as harvesting, planting, and reseeding (Kreider et al. 2022).
Forest recovery after wildfire is broadly evident in the Selway-Bitterroot and Bob Marshall Wildernesses, consistent with field surveys showing widespread conifer regeneration across a variety of fire histories.This study identifies some sites that are either failing to recover or have projected recovery times greater than 200 years, most of which have occurred since the mid-2000s.Continued monitoring is required to determine if these sites will track similarly to burn scars of the 1980s and 1990s or if future environmental conditions will fundamentally alter their trajectories.However, in the broader context of the nearly 1.2 million ha of burned forest analyzed for 35 years, we conclude that forests in the study area are currently recovering back to a forested state.Given the high likelihood of a warmer, drier future, the primary strength of this study is in defining baseline conditions and providing a framework for monitoring to detect changes in patterns of recovery as early as possible.

Fig. 1
Fig. 1 Map of the study area with fire history polygons, which includes the Selway-Bitterroot (SBW, top left) and Bob Marshall Complex (BMC, top right).Darker colors correspond to areas that reburned during the study period.The gradient-boosted machine learning model was trained on ten LiDAR sites shown in dark grey on the bottom map.The study area falls in Idaho and Montana, United States Fig. 3Example time series of the Conger Creek Fire.Conger Creek yielded four separate units of analysis: a once-burned unit that burned in 2007, a twice-burned unit with fires in 2007 and 2012, a twice-burned unit with fires in 1988 and 2007, and thrice-burned area with fires in1988, 2001,  and 2007.This time series depicts the once-burned unit.The pre-fire canopy cover was 63%.Fifty-nine percent of pre-fire canopy cover was lost during the fire (impact).Absolute minimum canopy cover was 25.8%, which was observed 4 years post-event.After this point, the site began recovering (gaining canopy cover).It is expected to take 76 years to recover to its pre-fire canopy cover

Fig. 4
Fig.4Spearman's correlation between the proportion of a site burned at high severity (Y axis) and pre-fire mean canopy cover (left), fire impact (center), and years to begin recovering canopy cover (right).The box plots depict the distribution of each variable, with the notch in the box corresponding to the median.Impact is proportion of pre-fire canopy cover lost and includes both immediate and delayed loss

Fig. 5
Fig. 5 Probability of sites recovering by ignition year, determined by logistic regression.The beginning of recovery is defined as the point post-event at which a site begins gaining more canopy cover than it loses.Lines represent probabilities, and ribbons represent standard error.The grey band depicts the probability of recovery when observed for the entire available record, the length of which varies based on ignition year.The colored bands represent the probability of recovering when observed for 5, 10, or 15 years post-event.All probabilities are calculated by a binomial glm model

Fig. 6
Fig. 6 Time to recover by year of ignition in conventional box-and-whisker notation.Non-recovering sites are depicted as a simple count below the X axi.s

Table 1
Covariates used in model training and prediction.A median and maximum value was calculated for each Landsat covariate, for a total of 10 Landsat and 3 LANDFIRE observations per pixel.Only pixels with observations available for all 13 covariates were used