Vegetation Mortality within Natural Wildfire Events in the Western Canadian Boreal Forest: What Burns and Why?

Abstract

Wildfires are a common disturbance event in the Canadian boreal forest. Within event boundaries, the level of vegetation mortality varies greatly. Understanding where surviving vegetation occurs within fire events and how this relates to pre-fire vegetation, topography, and fire weather can inform forest management decisions. We used pre-fire forest inventory data, digital elevation maps, and records of fire weather for 37 naturally-occurring wildfires (1961 to 1982; 30 to 5500 ha) covering a wide range of conditions in the western Canadian boreal forest to investigate these relationships using multinomial logistic models. Overall, vegetation mortality related to a combination of factors representing different spatial scales. Lower vegetation mortality occurred where there was lower fuel continuity and when fires occurred under non-drought conditions. Higher classification accuracy occurred for class extremes of no mortality (i.e., unburned areas within the burn event) and high mortality; partial vegetation mortality classes were harder to distinguish. This research contributes to the knowledge required for natural pattern emulation strategies, and developing responses to climate change.

1. Introduction

The boreal forest provides a vast array of goods and services, including timber products, wildlife habitat, water quality and quantity, recreation, carbon sequestering, hunting and fishing opportunities, minerals, and gas and oil [1]. Concerns about the sustainability of the boreal forest and its ecosystem services [2] has inspired a management approach where knowledge of historical disturbance patterns provides the primary foundation for planning activities [3]. The natural pattern (NP) concept proposes that landscape compositions and structures that are similar to historic conditions are more likely to maintain historical levels of biodiversity [4]. Further, the NP concept is the cornerstone of ecosystem-based management (EBM) [5] and has been embraced by provincial regulators in Canada [6] and forest certification agencies [7].

In the Canadian boreal forest (hereafter termed “boreal”), wildfires are responsible for the majority of the historical disturbance patterns [8,9]. For many years, the boreal was referred to as a stand-replacing ecosystem, meaning that most fires are highly severe, killing a large percentage of the trees [8]. However, more recent evidence suggests that historical boreal wildfires left a large amount of unburned and partially burned vegetation within the fire boundary [10,11]. Moreover, the amount and spatial patterns of surviving remnant vegetation within the area of a wildfire event (hereafter referred to as “remnants”) can vary greatly both between and within fires. For example, remnants in wildfires in the Boreal Plains Ecoregion tend to create stepping stones, while those of the Foothills Ecoregion tend to form spatially contiguous corridors [12]. The amount, but also the spatial arrangement of remnants is important for forest succession, habitat, refugia, and revegetation patterns [13,14].

In general, our understanding of how much vegetation survives in boreal wildfires far exceeds that of the spatial patterns of remnants and the associated underlying factors. Survival levels within boreal fires have been linked to fire weather, size and duration, pre-fire forest structure, vegetation types, topography, the locations and density of water bodies, and fuel types and arrangement [15,16,17,18,19]. Despite the breadth of research, the natural range of variation (NRV) of both the amount and spatial locations of remnants within boreal fires remains largely unexplained.

One possible explanation for this limited understanding is under-estimating the underlying complexity of the relationships between vegetation survival and causal factors. Across large landscapes and over several decades, individual fires respond to local differences in fire weather, topography, ecological zone, and fuel-type conditions [12,20,21]. This means that capturing the NRV of wildfire patterns requires extensive information across very large areas. To address this, a number of studies have used forest inventory data (e.g., [22]), the national wildfire database (e.g., [23]), and satellite imagery to capture NRV of historical fire patterns (e.g., [15,24]). While all of these data are readily available, none of these data capture smaller remnants and/or partial mortality. Recent studies have shown that this fine-scale complexity accounts for a significant part of wildfire remnants in the western boreal [12]. To capture the high complexity and variability of remnant patterns among fire events, large sample sizes are needed. While detailed analyses of burning conditions on individual or a small number of fires can provide valuable local insights (e.g., [17,25]), the findings of such studies are only a subset of NRV. Similarly, studies that include only a few explanatory variables (e.g., [26]) could lead to misleading conclusions about which factors are most important with regards to vegetation survival.

Further contributing to the challenges of understanding wildfire burning patterns is the inconsistent terminology. The meanings of “fire severity”, “burn severity”, “vegetation survival” and “remnants” not only differ but may also vary from one study to another. Many studies using remotely sensed data use the Normalized Difference Vegetation Index (NDVI) and Normalized Burn Ratio (NBR) to represent fire severity and burn severity, respectively. These indices capture vegetation mortality, but also reflect losses of dead wood and organic matter, ash deposition, and changes to soil structure and chemistry [18,27]. The Composite Burn Index (CBI) also captures burn severity, but does not capture vegetation mortality levels consistently. Jain and Graham [28] found that burn severity can vary greatly for canopy versus ground levels. Definitions of wildfire remnants have ranged from large, entirely unburned islands of trees within a burned polygon (e.g., [22]) to all types and levels of surviving vegetation within the general vicinity of a wildfire event (e.g., [23]). Andison [29] found that a three-fold difference in remnants can occur with even very subtle definition differences. As an example, Eberhart and Woodard [30] found that remnants comprised 5%–15% of the area of historical wildfire events, compared to 20%–60% by Andison and McCleary [12] for fires from a similar study area. These inconsistencies can lead to pattern artefacts [31], can reduce the ability to compare and compile NRV knowledge across studies, and create unnecessary confusion.

In this paper, we examined relationships between physical and environmental variables linked to fire behaviour at a range of spatial scales and levels of vegetation mortality within historic fire events to inform NRV management strategies in the western Canadian boreal. For this, we used an existing, spatially extensive, high resolution photo-based historical wildfire dataset from western boreal Canada that also includes pre-fire vegetation data. We then obtained information on a number of possible biotic and abiotic factors related to fire-caused vegetation mortality levels including topography metrics and local fire weather data. Using these variables at the fire event and within-event spatial scales, we built multinomial logistic (i.e., multinomial logit) models to predict vegetation mortality classes as a mechanism for examining the roles of these factors in determining the amounts and locations of post-fire remnants.

2. Materials and Methods

2.1. Study Area

The historic wildfire data used for this study was obtained from an existing database of 129 fires distributed across the Boreal Plains, Canadian Shield, Foothills, and Rocky Mountains Ecoregions of Canada (Figure 1).

The Boreal Plains Ecoregion is generally a flat or gently rolling landscape with thick glacially shaped soils and subtle upland areas interspersed with lower-lying wetlands. In the south, pure aspen (Populus tremuloides (Michx.) and other Populus spp.) stands are common, with patches of white spruce (Picea glauca (Moench) Voss) associated with moist soils. Wetlands are often dominated by shrubs. Further north in the Boreal Plains Ecoregion, upland vegetation is composed of mixed stands of aspen, balsam poplar (Populus balsamifera (L.)), and white spruce, with extensive lowland wetland areas with black spruce (Picea mariana (Mill.) B.S.P.), larch (Larix laricina (Du Roi) K. Koch), shrubs, and sedges. Throughout the region, stands of jack pine (Pinus banksiana (Lamb.)), lodgepole pine (Pinus contorta (Douglas) var. latifolia (Engel.) Critchfield), and hybrids of these two species can occur on sandy well-drained soils [32,33]. The Canadian Shield is characterized by glacially scoured granite bedrock, often barren and exposed at the surface or covered by thinner soils, resulting in rolling hills, coniferous or mixed forest stands where the soil is sufficiently thick, and numerous lakes and wet areas. Black spruce dominates the area, mixing with jack pine on uplands, and larch on lowland areas. Hardwoods and mixed stands are found on certain sites with favourable soil and aspect [32,33]. The Foothills Ecoregion is distinguished from the Boreal Plains Ecoregion by steeper slopes, higher altitudes, and more common lodgepole pine stands. Upland stands are typically deciduous or mixed-wood, with wetlands dominated by stunted black spruce, larch, or shrubby vegetation. In the Rocky Mountains, there is mountainous terrain with a larger range of elevation, steeper slopes, and orthographic precipitation. At low elevations, closed conifer stands are typical, while at higher elevations, more open conifer stands and herbaceous meadows, grasslands, or low shrubs typically occur. At high elevations, the growing season may be too short to support trees, resulting in alpine and subalpine environments. Given the steep slopes and large terrain features, vegetation cover can be strongly influenced by slope location and microsite (in particular, the availability of moisture and solar energy); for example, open coniferous stands with grasslands may be more commonly found on south facing slopes and closed conifer stands on north facing slopes [33].

2.2. Data

Thirty-seven fire events had available data on pre-fire vegetation, topography, and fire weather. Fire events in our sample were distributed across four Ecoregions: Boreal Plains (24); Canadian Shield (seven); Foothills (three); and the Rocky Mountains (three) (Figure 1).

2.2.1. Mortality Maps

The wildfire database used for this study had strict eligibility criteria, specifically:

• No evidence of pre-fire anthropogenic activity;
• Minimal or no fire suppression;
• No post-fire salvage logging; and
• High quality, high resolution aerial photos available within five years post-fire (see [29] for details).

The photo negatives required to cover the area of each fire were obtained and digitally scanned at 10 μm and imported into Softcopy [36]. The scale of the aerial photographs was either (a) better than 1:20,000 or (b) large format plates greater than 1:31,860 of sufficient quality to resolve individual trees (see [12,29] for details). In the Softcopy environment, photo-pairs were used to define the outer boundaries of each fire (i.e., the “shell”), and map all remnants using six mortality classes: class 0 = no mortality, 0%; class 1 = 1% to 25%; class 2 = 26% to 50%; class 3 = 51% to 75%; class 4 = 76% to 94%; and class 5 = ≥95%. Vegetation mortality was measured as:

• The percent dead tree crowns attributed to the fire event for forested areas;
• The percent cover of dead shrubs and bushes for non-forested areas with other woody plants; or
• The percent area scorched for non-forested areas with only grass or bryophytes.

For consistency, the same person interpreted all fire events. A minimum mapping unit of 10 m × 10 m was used, representing an area occupied by no more than three tree crowns. The resulting mortality maps were scanned, digitized, and spatially registered to 1:50,000 topographic maps. Thus the remnants tested in this study are the equivalent of “island remnants” sensu Andison [29].

2.2.2. Pre-Fire Vegetation

Pre-fire data were collected using high resolution photos acquired within five years of the fire event [29]. The photo negatives covering each fire were digitally scanned at 10 μm. The extent of vegetation interpretation included a 100 m buffer on the outside of the shell of each fire to capture the boundary area. Interpretation protocols followed the standardized vegetation inventory methods as defined by both Alberta and Saskatchewan [37,38]. The parameters used to define polygons included percent crown cover, species composition, height class, and age class for forested areas, vegetation type for non-forested areas, and soil moisture class. Since the specific criteria for photo-interpretation varied by province, variables chosen for this study were harmonized (

Table 1. Biotic and abiotic attributes.

Type	Variable	Description
Site	Ecoregion	Ecological land classification [32,33]
	Event area (m2)	Event area within boundary defined by vegetation mortality including unburned islands.
Vegetation	Age (years)	Year of the fire—estimated year of stand origin (photo-interpreted). Zero for non-forested polygons.
polygon	Fuel category	See Table 2.
	Height (m)	Average height of the dominant and codominant trees of the leading species. Zero for non-forest.
	Overstory crown closure (%)	Percent ground area covered by crowns of the dominant and codominant trees. Zero for non-forest.
	Soil moisture class	Interpreted classification of dry, mesic, or wet.
	Understory crown closure (%)	Cover of trees forming a distinct understory layer. Zero for non-forest or no understory layer.
	Understory height (m)	Average height of trees forming a distinct understory layer. Zero for non-forest or no understory layer.
	Ladder fuel index	Ratio of overstory and understory height multiplied by understory crown closure (following Alexander et al. [39]).
Topography: All variables calculated using the 30 m Digital Elevation Model (DEM).	Elevation (m)
	Slope (degrees)
	Curvature	Surface curvature. Positive is convex (ridges), zero is flat, and negative is concave (valleys) [40].
	Solar aspect index (TSRAI)	Circular aspect: 0 on northeast slopes, 0.5 for flat ground to 1 on southwest slopes [41].
	Compound topographic index (CTI)	Potential soil moisture, based on slope and catchment area. High values have potentially high soil moisture [42].
	SCOSA + SSINA	Slope * cosine (aspect) and slope * sine (aspect). Represents the “northness” and “eastness” of a pixel, respectively [43].
	Slope position	Relative elevation. Areas with values higher than 1 are higher than the surroundings and vise-versa [44].
Fire weather	Fine fuel moisture code (FFMC)	Moisture content of litter and fine materials indicating the ease of ignition and flammability of fine fuels.
	Duff moisture code (DMC)	Related to the average moisture content of organic soil layers to a moderate depth.
	Drought Code (DC)	Indicates seasonal drought and the expected smouldering in deep organic soil layers and large logs.
	Initial Spread Index (ISI)	Expected rate of fire spread given wind and FFMC. ISI threshold is the proportion of fire days with ISI > 8.7 [45].
	Build Up Index (BUI)	A rating of the total amount of fuel available for combustion by combining DMC and DC.
	Fire Weather Index (FWI)	Index of potential fire intensity given ISI and BUI. FWI threshold is the proportion of fire days FWI > 19 [45].

). Additionally, midpoints of classes were used to represent percent height class, age class, and crown cover instead of including these as class variables in all analyses.

Using these pre-fire vegetation and soil attributes, a fuel category was assigned to each vegetation polygon within the buffered area of influence of each fire based on the Canadian Fire Behavior Prediction (FBP) fuels types [46]. For forested areas, fuel categories were assigned based on overstory species and stand age; several new categories were assigned to represent local conditions (i.e., FX for balsam-fir (Abies balsamea (L.) Mill.) and LX for larch (

Table 2. Classification of fuel category based on FBP fuel type [46].

Type	Variable	Description
Forest	C2	Boreal spruce: Spruce (Picea spp.) as the leading species.
	C3	Mature pine: Pine (Pinus spp.) as the leading species, age > 40 years.
	C4	Immature Pine (Pinus spp.) as the leading species, age ≤ 40 years.
	D1	Leafless aspen: Trembling aspen (Populus tremuloides (Michx.)) and other types of aspen as the leading species, fire event started before June 1st (i.e., commonly before leaves flushed).
	D2	Leaf-on aspen: Trembling aspen (Populus tremuloides (Michx.)) and other types of aspen (Populus spp.) as the leading species, fire event started after June 1st (i.e., commonly after leaves flushed).
	FX	Fir: Balsam fir (Abies balsamea (L.) Mill.) as the leading species.
	LX	Larch: Larch (Larix laricina (Du Roi) K. Kosh) as the leading species.
	M1	Boreal mixedwoods, leafless: Conifer and deciduous species mixed in the overstory, fire event started before June 1st.
	M2	Boreal mixedwoods, leaf-on: Mixed conifer and deciduous species in the overstory, fire event started after June 1st.
Non-forest		Seven classes: Open shrub, closed shrub, open muskeg, treed muskeg, brush and alder, bryophytes, and grassland.

). Additionally, codes for non-forest vegetated areas were added.

2.2.3. Topography

Digital elevation data were downloaded to cover the buffered area of influence for each fire event from the Canadian Digital Elevation Dataset (CDED) which has a spatial resolution of 30 m (available from http://geobase.ca, accessed 27 October 2014). A suite of topographic indices were calculated that relate to elevation, slope, aspect, and landscape position (Table 1) using the Toolbox for Surface Gradient and Geomorphometric Modeling extension software version 2.0-0 for ArcGIS version 10.3.1 [40].

2.2.4. Fire Weather

Historic fire weather data from Amiro et al. [47] were used for this research. These data represent daily fire weather for Canadian wildfires >200 ha from 1959 to 1999, created by spatially interpolating weather station measures for 21 days following the recorded fire ignition date. Using these data, six daily fire weather indices were calculated based on the moisture availability of fuels and fire behaviour (Table 1). For this study, only three of the six available fire weather indices were selected, representing different aspects of fire behaviour [45,48]: ISI, to represent short-term changes in the availability of fine fuels; DC, to represent longer-term seasonal weather conditions and the availability of deep organic soil and large logs; and FWI as an overall index of potential fire energy output. For ISI and FWI, Podur and Wotton [45] calculated threshold values that distinguish fire spread events, where fire intensity is high, tree crowns are consumed, and the rate of spread is rapid, and non-spread events, where fire intensity is low, most burning occurs on the ground surface, and the rate of fire spread is low. In this study additional summary statistics were calculated for ISI and FWI as the proportion of days where each index exceeded the respective threshold for spread conditions (ISI = 8.7 and FWI = 19; following Podur and Wotton [45]). For each of these indices, daily values were summarized (i.e., median, maximum, minimum, and range) for either the duration of the fire minus one day, or for the full 21-days from the beginning of the fire event if the event was longer than 21 days.

2.2.5. Other Data

Other data collected and used in the analyses included Ecoregion, as defined by a blending of the Alberta and Canadian Ecoregions (see [12] for details). The event area was calculated (including unburned islands) and assigned to each pixel within each respective fire event.

2.3. Multinomial Logistic Models for Vegetation Mortality Class

For each fire event, the vegetation mortality class and all described biotic and abiotic attributes were tabulated by intersecting all layers with the center point of 30 m pixels, which was the lowest common spatial resolution of the vegetation polygons, vegetation mortality, and topography data. Although originally included in the data, pixels from the 100 m buffer around each fire boundary were not used for model training and testing because large portions of these boundary regions were subject to fire event-ending weather conditions that were different than the summary fire weather. From the remaining 433,475 pixels, 60% were randomly selected for model fitting to provide a computationally manageable balance between the data used for training and testing, with the majority used to train the model. Multinomial logistic (also known as multinomial logit) models were then fit to predict vegetation mortality class from subsets of biotic and abiotic attributes. These models provide cumulative probabilities of each class (i.e., class 1 versus others, classes 1 and 2 versus others, etc.) which can then be used to obtain the probability of each class via subtraction [49]. Because the focus of this study was to investigate which biotic and abiotic attributes affect vegetation mortality levels within fire events, equal weights were assigned to each mortality class in model fitting. Since there were a large number of attributes that might affect vegetation mortality and these attributes are likely related, particularly with the groups defined in Table 1 (i.e., vegetation, topography, and fire weather), the process used to select among predictor variables to obtain a final model was:

• Vectors of means (for continuous variables in Table 1) and proportions (for ecoregion plus other class variables in Table 1) were calculated by vegetation mortality level using all pixels to investigate univariate relationships with vegetation mortality and to indicate importance of each possible predictor variable. The results from this were used to provide an overview that described the general relationships.
• A correlation matrix was used to identify pairs of variables with moderately high multicollinearity (r > 0.7) as an aide to pre-selection among highly related variables. Using this as a guide along with interpretability of relationships as reported in other studies, a subset of the fire weather variables was retained in further analyses, namely: DC (median) to represent the drought conditions for the fire event; FWI (min) to represent the minimum potential energy output during the burning period; and ISI (proportion of days above threshold) to represent the proportion of days during the burning period that had high intensity fire-spread conditions [45]. Similarly, from the topography variables, slope location and surface curvature were strongly correlated; slope location was selected since it related to biomass consumption in fire events in previous studies [50]. From the vegetation variables, the understory index was strongly correlated with understory cover, so the simpler and more interpretable understory cover variable was retained.
• Following Hosmer et al. [51], multinomial models were initially fitted using each predictor variable (i.e., univariate models) and ranked based on the Akaike Information Criterion (AIC) [52] along with the percent of correctly classified pixels using the model fitting data.
• For each variable group (i.e., site, vegetation, topography, and weather), a stepwise selection process was followed by including all variables of the variable group and then removing variables one at a time. Variables previously dropped were then considered for entry back into the model at later steps. Again, AIC and the percent of correctly classified pixels were used to evaluate variable importance (i.e., whether to retain or drop a variable). However, supporting literature regarding relationships between mortality and biotic and abiotic variables was also used in deciding whether to retain or drop a variable.
• Since class variables can affect the coefficient associated with each continuous variable, interaction terms were then added and evaluated as in Step 4. Some interactions resulted in singular matrices and were removed. For example, some fuel categories only occurred within given ecoregions (e.g., muskeg did not occur in the Rocky Mountains or Foothills), and others only occurred within certain terrain (e.g., muskeg only occurred on flat ground) resulting in no differences in some attributes (i.e., all grasslands had an overstory height of 0).
• Once a subset of predictor variables was selected from each group of variables, these were merged together to obtain a model using all variable groups. Interactions between class variables and continuous variables across variable groups were then evaluated with regards to model improvements.
• Since predictor variables eliminated in previous steps might become important in a later step (for example, in combination with ecoregion or other variables), these were added again to the overall combined model and evaluated.

The vegetation mortality class with the highest predicted probability using the selected multinomial logistic model was assigned for the classification. The model that was finally selected was further assessed using confusion matrices for the reserved test data (40% of pixels within the fire event). The kappa coefficient was calculated [53] to indicate classification performance relative to a random assignment of pixels (−1 to +1; +1 indicates perfect agreement; <0 indicate performance is no better than random). For the final mode, the weighted kappa coefficient [54] was calculated to assess the degrees of agreement between cells in the confusion matrix.

3. Results

3.1. General Relationships

The highest vegetation mortality class accounted for 68% of pixels overall, although that ranged from 64% in the Boreal Plains to 82% in the Foothills (

Table 3. Numbers of fire events and pixels by Ecoregion, and percentages by vegetation mortality class.

	Number of		Vegetation Mortality
Ecoregion	Fire Events	Pixels	0%	1%–25%	26%–50%	51%–75%	76%–94%	≥95%
Boreal plains	24	246,951	3%	5%	7%	8%	13%	64%
Canadian Shield	7	166,425	1%	5%	7%	7%	8%	72%
Foothills	3	8378	2%	7%	7%	2%	0%	82%
Rocky Mountains	3	11,721	3%	5%	6%	2%	5%	78%
All	37	433,475	2%	5%	7%	7%	10%	68%

). Fires in the Boreal Plains and Canadian Shield had the highest levels of partially burned areas (33% and 27% respectively, compared to 16% and 18%) (Table 3). No clear trend was evident for fire event area.

For the variables from the vegetation polygons, the only variable that showed a consistent trend was overstory crown closure, which increased with higher mortality classes ranging from 42.5% in mortality class 0 to 49.6% in class 5 (See

Table A1. Means (standard deviations) of vegetation variables (see Table 1 for definitions) by vegetation mortality class.

Variable	Vegetation Mortality
	0%	1%–25%	2%–50%	51%–75%	76%–94%	≥95%
Age (years)	39.9 (20.4)	42.6 (24.4)	40 (21.9)	41.4 (20.4)	41.7 (19.8)	38.8 (21)
Height (m)	11.5 (5.4)	9.9 (4.4)	9.5 (4.3)	10.3 (4.6)	11.0 (4.7)	10.4 (3.7)
Ladder fuel index (no units)	8.2 (13.8)	6.3 (10.6)	6.6 (10.1)	6.4 (9.6)	6.1 (9.3)	8.5 (11.6)
Overstory crown closure (%)	42.5 (23.6)	43.6 (23.9)	43.2 (21.9)	46.4 (22.1)	49.3 (22.2)	49.6 (23.8)
Understory crown closure (%)	14.3 (20.4)	11.9 (18.8)	13.1 (18.7)	12.4 (17.1)	12.2 (17.3)	16.8 (21.1)
Understory height (m)	3.3 (3.7)	2.5 (3.2)	2.6 (3.0)	2.9 (3.2)	2.8 (3.1)	2.9 (3.0)

). For age, height, understory height, and soil moisture, no obvious trend was evident. For ladder fuel index and understory crown closure, partially burned pixels were common. Partially burned pixels had lower ladder fuel index (6.1–6.6) than pixels with no (8.2) or high (8.5) mortality. Partially burned pixels had also lower understory crown closure (11.9%–13.1%) than either high (16.8%) or low mortality (14.3%) pixels (Table A1).

Sixty-eight percent of all forested pixels were in the highest mortality class, although the numbers varied by fuel-type. Immature pine (C4) had the highest proportion of pixels with high vegetation mortality at 79%. Boreal spruce (C2), mature pine (C3), and boreal mixedwoods (M1 and M2) pixels had between 67%–71% of their pixels in the high mortality class. The lowest levels of high vegetation mortality occurred in fir (FX) and larch (LX) (49% and 51%, respectively) although the number of pixels in FX was very small (

Table A2. Relative frequency (as a %) and percentage of pixels by vegetation mortality class for each fuel category (see Table 2 for definitions).

Fuel Category	Vegetation Mortality
	0%	1%–25%	2%–50%	51%–75%	76%–94%	≥95%	Frequency
Forest
C2—Boreal spruce	2%	6%	6%	6%	9%	71%	27.62%
C3—Mature pine	2%	5%	6%	8%	11%	67%	17.73%
C4—Immature pine	1%	4%	4%	5%	7%	79%	8.1%
D1—Leafless aspen	5%	11%	4%	4%	17%	59%	0.83%
D2—Leaf-on aspen	1%	7%	15%	4%	6%	67%	0.41%
FX—Fir	16%	16%	6%	13%	0%	49%	0.03%
LX—Larch	1%	4%	14%	16%	14%	51%	6.95%
M1—Boreal mixedwoods, leafless	1%	3%	8%	7%	11%	70%	2.81%
M2—Boreal mixedwoods leaf-on	3%	4%	6%	8%	12%	67%	26.23%
All forest	2%	5%	7%	8%	11%	68%	90.71%
Non-forest
Open shrub	1%	2%	0%	1%	0%	96%	0.06%
Closed shrub	6%	15%	14%	10%	13%	42%	2.27%
Open muskeg	10%	10%	3%	9%	5%	62%	1.21%
Treed muskeg	2%	6%	7%	4%	10%	72%	5.49%
Brush and alder	0%	15%	0%	25%	0%	60%	0.06%
Bryophytes	0%	0%	0%	0%	100%	0%	0%
Grassland	1%	4%	21%	28%	3%	43%	0.2%
All non-forest	4%	9%	8%	7%	10%	63%	9.29%
All forest and non-forest	2%	5%	7%	7%	10%	68%	100%

).

Non-forested fuel types within the highest mortality class accounted for 63% of the pixels, which is 5% lower than the overall percentage for the forested class (Table A2). Although the FBP fuel-types do not differentiate non-forested fuel types, the percentage of high mortality pixels in closed shrub is only 42%, compared to 62% for open muskeg, and 72% for treed muskeg (Table A2).

Areas with no vegetation mortality were observed at higher elevations (>700 m) than areas with partial and high mortality (<615 m) (

Table A3. Means (and standard deviations) for topography variables (see Table 1 for definitions) by vegetation mortality class. Slope position index, surface curvature, and SCOSA + SSINA were scaled by a factor of 1000; TSRAI was scaled by a factor of 10.

Variable	Vegetation Mortality
	0%	1%–25%	2%–50%	51%–75%	76%–94%	≥95%
Elevation	711.8 (234.5)	561.4 (249.3)	606.1 (262)	614.1 (252.1)	611.5 (260.4)	598.5 (268.3)
CTI	9.0 (2.1)	9.5 (2.3)	9.1 (2)	9.1 (2)	9.0 (2)	8.9 (2)
Slope (degrees)	2.9 (2.6)	2.1 (3.2)	2.5 (3.3)	2.5 (3.2)	2.8 (3.7)	2.9 (3.6)
Slope position	−7.2 (145.7)	−7.0 (161.8)	−3.0 (188)	1.4 (174.8)	2.7 (183.8)	11.0 (196.7)
Curvature	−4.1 (84.2)	−3.9 (93.8)	−1.7 (108.5)	0.8 (101.1)	1.5 (106.5)	6.2 (113.4)
TSRAI	5.1 (3.6)	5.5 (3.2)	4.9 (3.3)	5.1 (3.4)	5.1 (3.3)	5.3 (3.4)
SCOSA	−7.1 (46.4)	−7.4 (49.2)	−4.7 (49)	−2.1 (39.4)	−5.8 (57.4)	−7.1 (56.8)
SSINA	5.7 (48.8)	−1.3 (46.7)	3.0 (55.8)	−6.1 (62.2)	−0.8 (62.7)	−1.5 (61.8)

). Areas with steeper slopes were associated with both unburned areas and areas with very high vegetation mortality (2.8–2.9 degrees), while slopes ranged between 2.1 and 2.5 degrees in areas mortality classes 1, 2, and 3 (Table A3). More prominent slope positions (i.e., + slope position values) and convex slopes (i.e., + curvature values) had higher vegetation mortality (classes 3, 4, and 5) than less prominent positions and concave slopes (e.g., gullies and valleys) (Table A3). For SCOSA and SSINA, representing “eastness” and “northness”, respectively (see Table 1. for definitions), southeast slopes (negative SCOSA and positive SSINA) were associated with unburned areas, and south slopes associated with higher vegetation mortalities (negative SCOSA, SSINA averaging near zero). There was no obvious pattern between vegetation mortality and CTI or TSRAI (Table A3).

In terms of fire weather variables, high median, minimum, and maximum values of ISI, DC, and FWI were associated with low to moderate levels of vegetation mortality (

Table 4. Means (and standard deviations) of daily fire weather variables (see Table 1 for definitions) by vegetation mortality class. For all variables, higher values indicate more favourable burning conditions; for the range variables, positive numbers indicate changes toward more favourable burning conditions and negative numbers indicate changes to less favourable burning conditions.

Variable	Vegetation Mortality
	0%	1%–25%	26%–50%	51%–75%	76%–94%	≥95%
ISI Median	4.8 (2.6)	6.2 (4.2)	5.4 (3.4)	5.2 (3)	4.5 (2.2)	4.3 (2.1)
ISI Min	1.6 (2.3)	2.7 (4.7)	1.6 (4)	1.3 (3.3)	0.8 (1.9)	0.7 (1.7)
ISI Max	15.5 (3.2)	13 (4.2)	12.7 (3.5)	13.5 (3.8)	13.4 (4.4)	12.2 (4.1)
ISI Range	−13.1 (5.7)	−7.7 (7.8)	−9.6 (6.7)	−10.7 (7)	−10.7 (8)	−8.7 (8.5)
ISI Threshold	0.2 (0.2)	0.3 (0.3)	0.2 (0.2)	0.2 (0.2)	0.1 (0.1)	0.1 (0.1)
DC Median	207.2 (61.9)	254.2 (118.9)	290.1 (110.6)	282.9 (107.2)	275.1 (105.1)	276.2 (104.8)
DC Min	166.1 (56.4)	217.3 (109.4)	247 (105.7)	238.5 (103.9)	225.6 (103.4)	229.6 (102.2)
DC Max	15.5 (3.2)	13.0 (4.2)	12.7 (3.5)	13.5 (3.8)	13.4 (4.4)	12.2 (4.1)
DC Range	62.9 (51.5)	52.6 (65.4)	74.9 (58)	80.4 (53)	79.6 (62.6)	70.4 (76.2)
FWI Median	14.3 (6.9)	17.5 (8.4)	16.5 (6.2)	16.3 (6.1)	14.2 (6.1)	13.8 (6.3)
FWI Min	5.1 (7.4)	7.8 (10.1)	4.7 (8.1)	4.1 (7.2)	2.8 (4.8)	2.7 (5.1)
FWI Max	35.8 (5.3)	31.5 (8.6)	32.5 (6.4)	34.1 (6.4)	33.1 (8.2)	31.2 (8)
FWI Range	−29.2 (11.9)	−18.1 (18.2)	−24.3 (16.7)	−26.1 (17.1)	−25.2 (18.9)	−20.8 (21.2)
FWI Threshold	0.3 (0.3)	0.4 (0.4)	0.4 (0.3)	0.4 (0.3)	0.3 (0.2)	0.3 (0.2)

). Wider ranges of ISI, DC, and FWI were associated with both unburned areas and increasing levels of vegetation mortality. For both the ISI and FWI threshold values, a greater number of days above the threshold for spread conditions was associated with low and moderate levels of vegetation mortality (Table 4).

3.2. Model Fitting and Classification

For the site variables, preliminary analyses indicated a nonlinear trend between vegetation mortality classes and the fire event area. As a result, the logarithm of fire event area was used in all models as in previous studies (e.g., [15]). The model using fire event area, Ecoregion, and the interaction term between them resulted in the best model fit (i.e., lowest AIC relative to the null model) and classification performance (

Table B1. Fit statistics for models using Ecoregion and fire event area (area).

Model	Variables a	AIC	Δ AIC b	% Correct	Kappa
1 (Null)		932,037	0	17%	0.00
2	Ecoregion	929,929	2107	30%	0.03
3	Area	932,030	7	32%	0.01
4	Ecoregion, area	929,230	2807	38%	0.04
5	Ecoregion * area	920,267	11,769	37%	0.07

^a Interaction terms indicated by *. For example x * z means that the model includes variables x and z, as well as the interaction between them; ^b Δ AIC is relative to the null model with a larger value indicating a better model.

, Model 5). The interaction term indicates that the effects of fire event area vary with Ecoregion. Based on these results, both Ecoregion and fire event area were considered for inclusion in multivariable models using variables across variable groups.

Of the models using a single vegetation variable, fuel category had the best fit (

Table B2. Fit statistics for models using vegetation polygon variables.

Model	Variables a	AIC	Δ AIC b	% Correct	Kappa
1 (Null)		932,037	0	17%	0.00
2	Fuel category	925,179	6858	20%	−0.01
3	Overstory crown closure	928,397	3640	36%	0.02
4	Soil moisture class	931,447	590	44%	0.02
5	Understory crown closure	931,840	197	34%	0.02
6	Understory height	931,912	125	36%	−0.02
7	Age	931,918	119	26%	0.01
8	Overstory height	931,948	88	32%	0.00
9	Fuel category, overstory crown closure, understory crown closure, understory height, age, overstory height	917,106	14,930	43%	0.07
10	Fuel category, age, overstory crown closure * overstory height, understory crown closure * understory height	915,557	16,479	42%	0.07
11	Fuel category, overstory crown closure, understory crown closure, understory height, age	917,529	14,507	43%	0.07

^a Interaction terms indicated by *. For example x * z means that the model includes variables x and z, as well as the interaction between them; ^b Δ AIC is relative to the null model with a larger value indicating a better model.

, Model 2); however, the classification performance was poor. Despite the poor classification, fuel category was included for consideration in multivariable models, given this impacts fire behaviour [46]. Variables related to fuel continuity (overstory and understory crown closure) had poorer fits, but better classification results (Table B2, Models 3 and 5). As a result, the fuel continuity variables were also selected for possible inclusion in the multivariable models. Of the two variables representing soil moisture, CTI was selected over soil moisture classes since it is a continuous variable at a smaller spatial scale (Table B2, Model 4 vs.

Table B3. Fit statistics for models using topographic variables.

Model	Variables a	AIC	Δ AIC b	% Correct	Kappa
1 (Null)		932,037	0	17%	0.00
2	Elevation	930,250	1787	36%	0.02
3	CTI	931,265	771	41%	0.02
4	SCOSA, SSINA	931,726	310	29%	0.02
5	Slope	931,741	296	23%	0.00
6	Slope position	931,760	277	33%	0.01
7	TSRAI	932,038	−1	39%	0.00
8	Elevation, CTI, SCOSA, SSINA, Slope, Slope position	928,148	3889	33%	0.02
9	Elevation, CTI, SCOSA, SSINA, Slope	928,282	3755	34%	0.02
10	Elevation, CTI, SCOSA, SSINA	928,454	3583	34%	0.02
11	Elevation, CTI	929,482	2555	34%	0.02
12	Elevation * CTI	929,186	2851	32%	0.03

^a Interaction terms indicated by *. For example x * z means that the model includes variables x and z, as well as the interaction between them; ^b Δ AIC is relative to the null model with a larger value indicating a better model.

, Model 3). Understory height and age followed with regards to model fit. Both were considered for the multivariate models, given the impacts of ladder fuels on fire behaviour [55] and of stand age on fire impacts [19,56]. Finally, overstory height was selected for multivariable models, given the importance of crown structure for fire behaviour [55]. Using all of these selected variables resulted in a better fit than using individual variables (Table B2, Model 9). The classification performance of the model with all selected variables (Table B2, Model 9) was not improved by adding interaction terms between the canopy variables for the two respective layers (Table B2, Model 10). It was not possible to add further interaction terms for the fuel category variable, due to rank deficiency and linear combinations. Removing the overstory height variable slightly improved the model classification and resulted in fewer terms in the model (Table B2, Model 11).

The topographic variables with the best fit and classification accuracy were elevation and CTI (Table B3, Models 2 and 3). Using the two highest ranked individual variables together had nearly as good of a fit and classification performance as using all topographic variables (Table B3, Model 11). Adding an interaction term between the variables further improved fit and classification performance (Table B3, Model 12). Adding the previously removed variables to the model did not improve fit or classification accuracy.

Within the fire weather variable group, DC had the best model fit, while ISI had the best classification accuracy (

Table B4. Fit statistics for models using fire weather variables.

Model	Variables a	AIC	Δ AIC b	% Correct	Kappa
1 (Null)		932,037	0	17%	0.00
2	DC	923,964	8072	43%	0.04
3	FWI	924,210	7827	51%	0.05
4	ISI	926,688	5349	53%	0.06
5	DC, FWI, ISI	920,298	11,739	34%	0.02
6	DC * FWI * ISI	914,502	17,535	47%	0.03

^a Interaction terms indicated by *. For example x * z means that the model includes variables x and z, as well as the interaction between them; ^b Δ AIC is relative to the null model with a larger value indicating a better model.

). DC was the best variable for classifying unburned areas (Table B4, Model 2; kappa = 0.05, 0, 0, 0, 0, and 0.07 for classes 0 through 5, respectively), FWI was the best variable for classifying areas with intermediate levels of vegetation mortality (Table B4, Model 3; kappa = 0.03, 0, 0.03, 0.05, −0.04, and 0.10 for classes 0 through 5, respectively), and ISI was the best variable for classifying areas with high vegetation mortality (Table B4, Model 4; kappa = 0.01, 0, −0.02, 0, 0.02, and 0.15 for classes 0 through 5, respectively). All three variables were selected for the final multivariable model for this group of variables (Table B4, Model 5). Adding interactions between all selected variables improved model fit (Table B4, Model 6).

Combining selected variables from each variable group improved both fit and classification performance particularly when interaction terms between variables and Ecoregion were included (

Table 5. Fit statistics for models using combinations of variables across variable groups.

Model	Variables a	AIC	Δ AIC b	% Correct	Kappa
1 (null)		932,037	0	17%	0.00
2	Ecoregion * area, fuel category, overstory crown closure, understory crown closure, understory height, age, elevation * CTI, DC * FWI * ISI	888,482	43,555	40%	0.07
3	Ecoregion * (overstory crown closure, understory crown closure, understory height, age, elevation * CTI), DC * FWI * ISI	870,911	61,126	39%	0.10

^a Interaction terms indicated by *. For example x * z means that the model includes variables x and z, as well as the interaction between them; ^b Δ AIC is relative to the null model with a larger value indicating a better model.

, Model 3). No further improvements were obtained by including variables eliminated in previous steps. In particular, using soil moisture (classes) instead of CTI did not improve the model (AIC 1352 higher and equal overall kappa), nor did using ISI alone without the other fire weather variables (AIC 6649 higher and overall kappa 0.01 lower). Adding overstory height improved the fit (AIC 1744 lower), but not the classification performance (equal kappa). Adding additional topography variables made relatively small changes to the fit (changes in AIC < 1750) and did not improve the classification performance (equal kappa).

Using the final model (Table 5, Model 3) resulted in a classification accuracy of 39%. The percent of pixels correctly classified was best for the unburned and highest vegetation mortality categories (the kappa coefficients were 0.18, 0.06, 0.03, 0.03, 0.00 and 0.19 for class 0 through class 5, respectively). For intermediate levels of vegetation mortality, the classification accuracy was only slightly better than random, and class 4 had the worst classification accuracy (

Table 6. Confusion matrix for the final model using a combination of site, vegetation, topography, and fire weather variables.

Predicted Vegetation Mortality Class	Actual Vegetation Mortality Class
	0%	1%–25%	26%–50%	51%–75%	76%–94%	≥95%
0%	66%	31%	14%	16%	16%	7%
1%–25%	19%	19%	15%	13%	15%	9%
26%–50%	2%	7%	12%	13%	9%	8%
51%–75%	1%	5%	8%	8%	5%	5%
76%–94%	7%	19%	27%	24%	22%	23%
≥95%	4%	19%	25%	25%	33%	48%

). Many of the misclassifications were in adjacent categories. For example, the weighted kappa, which takes into consideration misclassifications in neighbouring classes by assigning a weighted penalty proportional to the squared distance to the correct class, was considerably better at 0.26. Finally, by counting classifications within ±1 category as correct, which is justifiable from a management standpoint, 64% of pixels were classified in the correct classes with kappa coefficients of 0.26, 0.29, 0.21, 0.37, 0.31, and 0.32 for mortality classes 0 through 5, respectively.

Application of the selected model to visually compare predicted to actual vegetation mortality classes for a randomly selected fire from each Ecoregion provided further insights (Figure 2). Overall, the actual and predicted vegetation mortality class maps are similar. However, the fire event in the Boreal Plains Ecoregion (Figure 2a) had long and narrow patches with high vegetation mortality (high length to breadth ratio, see [46]) oriented from southwest to northeast, suggesting that it may have been driven by high-velocity winds. Further, assuming fire propagation followed this wind direction, remnants formed on the lee-side of non-flammable areas, while high vegetation mortality occurred in windward locations. In this case, the direction of fire propagation may have also been important, given that it appears that live vegetation remnants formed while the fire was spreading in a downhill direction. However, without fire propagation maps and detailed information on wind direction, it is impossible to know how wind velocity, topography, and fire propagation affected vegetation mortality. The fire event in the Canadian Shield had large areas of very high vegetation mortality (Figure 2b). Within the classification, several low lying wet areas, indicated by predicted low-levels of mortality in typical meandering patterns at the bottom of drainages, were misclassified as having lower vegetation mortality than was observed (Figure 2f). Within the Foothills, the classification had very similar patterns of vegetation mortality; in particular, remnants appeared as linear shapes in both the observed and predicted maps. These remnants could serve as wildlife corridors and were present in both the observed and predicted maps (Figure 2c,g). Finally, for sampled fire events in the Rocky Mountains, the model appeared to make misclassifications within fuel types that may have been related to differences in aspect, which was not accounted for in the final model (Figure 2d,h).

4. Discussion

4.1. Factors Contributing to Vegetation Survival

The results suggest that a combination of factors influenced fire-caused vegetation mortality at two distinct scales. At the event scale, overall vegetation mortality was moderately influenced by Ecoregion reflecting the broad ecological patterns first identified by [12], and found by others [15,16]. Fires within both the Foothills and Rocky Mountains had lower proportions of remnants in partial mortality classes and higher proportions of high mortality than other Ecoregions. In terms of fire weather, the influence of DC likely reflects a logical relationship between fewer remnants and higher levels of seasonal drought conditions associated with lower moisture content of large woody debris and deep organic soils. However, the direct relationship between FWI minimum values and the amount of low to moderate vegetation mortality was unexpected. One possible explanation for this is the nature of diurnal burning patterns. For example, during weather conditions favourable to even high-energy output wildfires, patches with low to moderate vegetation mortality may be created at night under less favourable burning conditions (i.e., lower temperatures and higher relative humidity). Another possible explanation is that fuel-types that normally have a low risk of flammability (e.g., wetlands) will burn under highly favourable fire weather conditions, but even then, only partially.

Within events, the most influential factor affecting vegetation mortality was fuel type, which is not surprising given its prominence in fire prediction models (e.g., [46]). Combining fuel-type with forest structure variables such as crown closure further improved the ability to predict vegetation mortality levels, which is also consistent with what we know about the strong relationship between crown closure and crown fires [57]. Elevation and soil moisture were the next most important within-event scale variables, although one could argue that elevation also acts at the event scale. The relationship between vegetation mortality and elevation also varied greatly between Ecoregions: in the Canadian Shield, higher vegetation mortality was associated with higher elevations, while in the Foothills, higher elevations were associated with lower mortality. In the Canadian Shield, high elevation sites tend to be height-of-land ridges, peaks and knolls, which tend to be hot and dry. In the Foothills, “high-elevation-sites” is a relative term, in that all fire events lie within the eastern slopes of the Rocky Mountains, which tend to be colder and wetter.

4.2. Complexity of Relationships

Many of the relationships with remnants found in this study have been noted by others, and are intuitively consistent. For example, Ecoregion was previously associated with the proportion of vegetation mortality within fire events [12,26], certain fuel types and percent cover variables have been linked with vegetation mortality and burn severity [28,46,56], topography has been related to burn severity [15,50], and fire weather variables exhibited weak relationships to burn severity [21]. While this study confirms many of these previously found relationships, it reveals a new level of complexity in terms of the processes involved in the formation of remnants. The relative influences of, and interactions between, both large-scale and within-fire attributes were examined in this study. Rarely is there a single physical or environmental factor influencing the amount or location of remnants acting alone. Acknowledgment of this complexity helps guide both research and management in several ways. First, studying within-fire vegetation survival improves our knowledge of factors affecting fire-related mortality and thereby improves our ability to manage forests using the NP concept. To facilitate this, access to high resolution spatial data is an advantage in that it allows a more precise evaluation of these relationships. The use of high quality spatial data sets should become the rule, not the exception. Second, given the large variability of remnant patterns noted here both within and between fire events, extensive and detailed data sets that cover the range of conditions are needed across large areas in order to capture the NRV of remnant patterns at local and landscape scales. Finally, both factors that change slowly in time and those varying within and between years affect vegetation mortality within wildfires. Ecological type and topography variables change only over a very long time-scale. In contrast, fire weather attributes are highly variable in shorter time periods. Although improving our knowledge of within-fire vegetation mortality is critical for altering factors that impact vegetation mortality, some factors will always be beyond the control of forest managers (e.g., fire weather). However, a decision support tool that helps predict the impacts of fire weather and other factors on vegetation mortality would be very valuable given expected impacts of climate changes on future wildfires.

From a management perspective, decision support tools can and should differentiate between variation that can be accounted for (i.e., deterministic) and that which cannot (i.e., random). The deterministic part of the natural variability attributable to regional-scale factors such as ecological zone, and local-scale factors such as vegetation type and topography can be modelled, predicted and emulated by careful planning. The random portion of remnant NRV attributable to fire weather can be applied more broadly across a region or landscape as intrinsic variability. Integrating both types of variation into forest planning is consistent with an NRV approach. The division of remnant NRV may also have value as a tool for understanding and potentially mitigating the impacts of climate change. The relationship between remnants and fire weather may help us better predict how wildfire remnant patterns might shift in the future, and how forest management activities might offset those changes as desired.

4.3. Classification Accuracy and Model Complexity

While the overall classification accuracy found in this study is comparable to those from other studies, the comparison is misleading for several reasons. First, this study included fire events covering an extensive land area with diverse ecological, topographical, and climatic conditions. Other broad regional studies tended to report lower accuracy with more complex models [21,31,58] than studies using individual or a very limited numbers of fire events [25,50]. Second, the increased variation represented in this more extensive area coupled with finer-scale variables enabled a more thorough investigation of factors and highlights the importance of considering local and regional conditions. Lastly, we included six mortality classes, compared to two or three classes for most other comparable studies. If we ignore all of the partial mortality classes, our classification accuracy was much higher than any of those reported in other studies.

Our models struggled most with identifying partial mortality. Partial mortality is a well-documented weakness of imagery-based interpretations [24,59]. We had hoped to overcome this challenge by using more accurate and precise mortality data based on high resolution photographic negatives, the results suggest that there are other factors involved in predicting vegetation mortality levels. One possibility is that mortality continues over several years post-fire [60] which means that all post-fire “snapshots” (including those used in our study) represent just one of several remnant vegetation patterns over time. Another possibility is that moderate mortality is associated with areas that experienced less severe fire conditions, such as overnight fires and/or the back or flanks of a fire where the relative influence of fire weather, topography and fuel-type may change, or even reverse [25]. Finally, the factors affecting vegetation mortality may differ among levels. Fire weather conditions may be critical for the formation of moderate levels of vegetation mortality, while fuel-type and topography differences may over-ride fire weather and play a more dominant role for determining areas that are either unburned or burned with high mortality. The higher proportion of fires that burned with high levels of island remnants found in the Foothills relative to the Boreal Plains reported by Andison and McCleary [12] support this concept. In any case, given the prevalence of partial mortality noted recently in boreal wildfires, this would seem an obvious target for future research.

We also learned some important lessons from a methodological perspective. One of the issues that arose in this and other related studies is how to examine the many variables that may be related to fire-caused vegetation mortality. Random Forests [61] is a commonly utilized machine learning approach mainly used in conjunction with remote sensing-based response variables [16,21,25,50], while very few studies have used logistic regression and vegetation mortality-based response variables [57]. In contrast, the multinomial logit model approach used in this research required a deliberate and thoughtful process for the selection of dependent variables. While documenting this process can be tedious, it highlights the complexity of relationships for vegetation mortality and importance of interpreting the relationships. The process we used included separating variables into logical groups which also improved our ability to interpret relationships.

5. Conclusions

Understanding the processes involved in the formation of surviving remnant vegetation within boreal forest wildfires has been more challenging than we have collectively assumed. Many prior fire remnant vegetation studies relied on available spatial data (of varying quality), a small number of fire events, and only a few predictor variables, all of which could impact the results and interpretation, ultimately resulting in an incorrect representation of the true nature of historical, natural patterns. In this study, we used high resolution and high quality spatial data for a large number of fire events over an extensive area to improve post-fire remnant vegetation mortality prediction models over results previously reported for regional studies. Moreover, we demonstrated that vegetation survival occurs as a result of a complex interaction between several biotic and abiotic factors at different spatial and temporal scales. More specifically, we found that the most influential general variables such as Ecoregion and drought code (DC) interacted with fuel-type, canopy closure, and soil moisture in very specific ways. While the influences of some variables (such as soil moisture, fuel-type, and drought code) were universally applicable, their relative influence in our predictive model varied by fire event and ecological region.

The findings from this study suggest that adapting fine-scale forest management patterns that are truly ecosystem-based involves a more holistic interpretation of fire patterns knowledge. However well-intentioned, traditional regulatory models, that mandate the amount and location of remnants based on deterministic rules are unlikely to capture the full historic range. Historic bounds are defined very generally by smaller-scale, in situ attributes such as fuel-type, soil moisture and canopy closure, many of which can be manipulated by forest managers. However, our results suggest that a substantial portion of the natural variation may not be directly related to in-situ variables. Fire weather, acting in concert with in-situ factors, cannot generally be manipulated by forest managers. Managing towards achieving natural variation objectives will require a more results-based approach that is challenging to integrate into practice from a regulatory perspective, but at the same time, this will offer the flexibility needed for management of a broader range of ecosystem values.