Spatially-Explicit Prediction of Wildfire Burn Probability Using Remotely-Sensed and Ancillary Data

Abstract

Wildfire is a critical process shaping the structure and composition of forest landscapes of western Canada. Spatially-explicit forest disturbance history and forest structure estimated using remotely-sensed data enables the characterization of burn probability, defined as the susceptibility of landscapes to fire hazard over time. In this research, we leveraged the Landsat archive to determine the capacity of land cover, forest structure, and forest disturbance information, together with ancillary data, to estimate burn probability. We analyzed the interactions between a number of contributing factors and identified landscapes with high probability to burn across forested ecosystems of Saskatchewan, Canada. Overall, we found that forests composed of coniferous species, with shorter trees (<3 m), low canopy height variability, an open stand structure (<10% canopy cover), and low timber volumes (<50 m³/ha), had higher burn probabilities. A 2015 burn probability map indicated that forests that did burn in 2015 (determined using independently mapped wildfires) had a median predicted burn probability of 81%, while the median burn probability for unburned forest area was 19%. This paper demonstrates the potential to generate detailed and spatially-explicit burn probability maps from time series remote sensing data to inform wildland fire risk modeling, management, and mitigation.

RÉSUMÉ

Les feux de forêt façonnent la structure et la composition des paysages forestiers de l’Ouest canadien. L’historique des perturbations forestières sur le plan spatial et la structure forestière estimée à l’aide de données de télédétection permettent de caractériser la probabilité de combustion, définie comme la susceptibilité des paysages aux risques d’incendie au fil du temps. Dans cette recherche, nous avons tiré parti des archives de Landsat et de données auxiliaires pour déterminer la capacité de l’information sur la couverture terrestre, la structure forestière et les perturbations forestières d’estimer la probabilité de combustion d’une forêt. Nous avons analysé les interactions entre un certain nombre de facteurs contributifs et identifié des paysages ayant une forte probabilité de brûler dans les écosystèmes forestiers de la Saskatchewan, au Canada. Dans l’ensemble, les forêts composées de conifères, avec des arbres plus courts (3 m), une faible variabilité de leur hauteur, une structure de peuplement ouverte (10% de couverture) et de faibles volumes de bois (50 m³/ha), présentent des probabilités de combustion plus élevées. La carte de probabilité de combustion établie pour 2015 indique que les forêts qui ont brûlé effectivement en 2015 (selon une cartographie des feux de forêt déterminée indépendamment) ont une probabilité de combustion médiane de 81%, tandis que la probabilité pour la zone forestière non brûlée est de 19%. Cet article démontre le potentiel de générer des cartes détaillées et explicites des probabilités de combustion à partir de séries temporelles d’images satellitaires afin de faciliter la modélisation, la gestion et l’atténuation des risques d’incendie dans les forêts sauvages.

Introduction

Wildfire is considered the most dominant disturbance agent in western Canada (Flannigan et al. 2005; Wang et al. 2015; Wotton et al. 2017), regulating landscape structure, species diversity, forest development stage, as well as carbon and hydrological cycles (Kasischke and Stocks 2000; Wotton et al. 2010; de Groot et al. 2013). Recent estimates suggest that over a 30-year period, the mean annual area burned in Canada’s forested ecosystems for the period 1985–2015 was greater than 1.6 Mha, but varied markedly with large fire years often dominating the long-term trend (σ = 1.1 Mha; Coops et al. 2018). Under a changing climate, a number of wildfire characteristics such as area burned, size, occurrence, and seasonality are expected to change as weather becomes more fire conducive (Boulanger et al. 2014; Xi et al. 2019; Whitman et al. 2019), although the increase in area burned is expected to be gradual (Price et al. 2013). Already, there is some evidence of recent shifts in area burned consistent with these warmer climate patterns (Kasischke et al. 2010). As a result, scientists and policy makers are increasingly concerned with the dynamics and mechanics of forest fires, fire-induced ecosystem changes, as well as implications to forest management, especially with regard to carbon cycle and risk mitigation policies (Taylor et al. 2013; San-Miguel et al. 2017).

The factors that drive the probability, behavior, and effects of wildland fires are spatially and temporally variable (Gralewicz et al. 2012; Xi et al. 2019). The occurrence of fire can be attributed to weather and climate (Flannigan and Harrington 1988), the availability of fuel (Romme 1982; Finney 2001), and ignition sources, such as lightning (Malamud et al. 2005), as well as anthropogenic influences including accidental and intentional ignitions (Rollins et al. 2001). While the factors are well known, the stochastic nature of fire occurrence associated with local conditions including surface moisture and wind, are difficult to predict. In addition to wind and daily fire weather, landscape fuel conditions play a significant role in fire spread once ignited (Cheney et al. 1998; Fernandes 2001). These fuel conditions include the vegetation types as well as the spatial heterogeneity of these landscape characteristics, such as vegetation species and forest age (Ryu et al. 2007).

To understand the susceptibility of a landscape to fire hazard over a given period of time, burn probability models have been developed, encompassing a range of drivers, such as climate, topography, vegetation, and weather (Parisien et al. 2005; Finney et al. 2011). These models typically employ a simulation approach to quantitatively estimate the probability of burning in a spatially-explicit manner (Wang et al. 2016). For example, the Burn-P3 model (Parisien et al. 2005) makes use of an existing fire database and simulates the ignition and spread of fires over a large number of years. As a result, a burn probability surface can be derived from the frequency with which each grid cell burns over years of simulation. Burn probability models are sensitive to quality of input data, especially spatially-explicit fuel information, which has been reported as the most critical variable governing spatial patterns of burn probability (Parisien et al. 2013). Given that detailed, up-to-date, and accurate fuel data might not be available over large areas (Parisien, Dawe, et al. 2020), classifications are typically carried out to convert generic vegetation maps into fuel maps (Parisien et al. 2013). Statistical or machine learning models have also been used, to some extent, to model burn probability or burned areas, often at regional to continental scales (Moritz et al. 2012; Parisien et al. 2014; Davis et al. 2017). As a consequence, these models typically use open access vegetation products at coarse spatial resolutions to characterize vegetation, such as those derived from MODIS (Moderate Resolution Imaging Spectroradiometer) or AVHRR (Advanced Very-High-Resolution Radiometer). Additionally, these vegetation products, such as spatially generalized forest cover (Davis et al. 2017) or net primary productivity (Moritz et al. 2012) may not characterize fuel conditions with sufficient locational or categorical accuracy or to the extent possible with newly available datasets.

A variety of technologies and data sources have been used to map fire occurrence and spread across Canadian boreal forests, including GPS (i.e. fire boundaries recorded from aircraft), aerial photography, and satellite imagery (Stinson et al. 2011). As Canadian forests are composed of both managed and unmanaged land crossing provincial and territorial borders, multiple jurisdictions are involved in stewardship, leading to issues in both temporal and spatial compilations (Stocks et al. 2002; Boulanger et al. 2012; Bolton et al. 2019) with an emphasis on larger and more severe fires, potentially introducing bias in national level trends and statistics (Taylor et al. 2013; Anderson et al. 2015). Remote sensing approaches provide spatially-explicit information on land cover and surface characteristics in a repetitive fashion, which makes it suitable and cost-effective for mapping burned forests, especially over large areas (Chu et al. 2016; Soverel et al. 2010). Free and open access to the Landsat archive, in particular, has enabled the characterization of an entire area of interest with the same information source, thereby allowing for the transparent, systematic, and repeatable production of spatially exhaustive information products (Roy et al. 2010; Griffiths et al. 2013; White et al. 2014). Recent research has demonstrated that by mining this extensive archive of long term spectral observations of the Canadian boreal, disturbance objects, such as fire can be accurately identified, described, and quantified (Hermosilla et al. 2016; White et al. 2017; Coops et al. 2018), at a level of spatial detail commensurate with information needs supporting post-fire monitoring, reporting, modeling, and possible remediation.

This capacity to systematically and exhaustively map burned areas provides a marked improvement in our understanding of the spatial and temporal occurrence of wildfires and the related characteristics (Coops et al. 2018). In this study, we developed a burn probability model using a statistical approach for forested landscapes in Saskatchewan, Canada, rather than using simulation techniques as often applied in existing studies. The input of our model included an extensive, open access, database of Landsat mapped fires, and a suite of fire occurrence and spread drivers derived from remotely-sensed data, including land cover, forest structure, and forest disturbance information. Our specific objectives are to (i) assess the importance of these Landsat time series derived drivers and other environmental factors (biotic and abiotic) for characterizing burn probability; (ii) understand the nature of the relationship between wildfires and these factors, and (iii) showcase the wall-to-wall map of predicted burn probability and assess its agreement with actual burned area mapping.

Study area

The forested area of Saskatchewan, Canada, occupies over 37.8 Mha representing approximately 58% of the provincial land base (Figure 1). Saskatchewan covers four distinct ecozones: Taiga Shield, Boreal Shield, Boreal Plains, and Prairies ecozones (Ecological Stratification Working Group, 1995), with the first three above noted ecozones being forest dominated. Forests in Saskatchewan are dominated by conifer species, including black spruce (Picea mariana) and Jack pine (Pinus banksiana) in the north and quaking aspen (Populus tremuloides), eastern larch (Larix laricina), and white spruce (Picea glauca) in the south (Steele et al. 1997; Archibold et al. 2000). While southern Saskatchewan has a focus on agriculture, the more northern reaches of the province are sparsely populated and primarily characterized by unmanaged forests subject to wildfires (Stocks et al. 2002; Wulder et al. 2004).

Figure 1. (A) Study area in Saskatchewan, Canada, shown with false color (Red: shortwave infrared band, Green: infrared band, Blue: red band) Landsat best-available-pixel (BAP) composite relating 2015 conditions overlaid with forested ecozones. Inset map displays its location in Canada. (B) Spatial distribution of burned forests in 1985 to 1994 and 1995 to 2014 (the latter period was used to build the burn probability model) detected from Landsat time-series data in relation to unburned forests over the 1985 to 2014 period. (C) Topography of the study area. (For interpretation of the references to color in this figure legend, the reader is referred to the online version of this article.)

Methods

The methodology presented herein enables the characterization of burn probability through the integration of Landsat time series data and ancillary data (Figure 2). From these datasets we derived predictor variables characterizing both biotic (i.e., land cover, forest structure, spectral trajectory metrics, and presence of historical wildfires) and abiotic (topographic features and geographic coordinates) conditions of forests. Disturbance information derived from Landsat time series was used to inform and annually stratify the sample selection for the burn probability model generation.

Figure 2. Workflow for integrating Landsat and ancillary data for modeling burn probability.

Data

Spectral trajectory metrics, Forest disturbance information, and land cover

In order to examine the relationship between environmental conditions and forest fire, annual surface reflectance image composites from 1985 to 2015 were used. These image composites were derived in conjunction with forest disturbance information following the Composite2Change (C2C) approach (Hermosilla et al. 2016). Given that the open Landsat data archive contains imagery acquired by different sensors with varying acquisition dates and atmospheric conditions, a best-available-pixel (BAP) compositing technique was used to select pixels of the highest quality over our study area. For the time period of interest, we used surface reflectance data acquired by Landsat-5 Thematic Mapper (TM), Landsat-7 Enhanced Thematic Mapper Plus (ETM+), and Landsat-8 Operational Land Imager (OLI) imagery (Masek et al. 2006; Vermote et al. 2016). To best characterize forest environments, a target date of August 1st was used, ensuring that only pixels acquired within a ±30-day window around this date were included in the composite. Clouds and the associated shadows in each image were detected using Fmask (Zhu and Woodcock 2012). The best available observation for each pixel within our study area was selected based on the scoring functions defined in White et al. (2014) and refined in Hermosilla et al. (2016), which ranks pixels based on presence and distance to clouds and cloud shadows, proximity to target date, and acquisition sensor.

Following the development of annual BAP composites, spectral trend analysis was carried out at the pixel level to remove noise (e.g. haze, smoke) in the trajectory and infill data gaps with interpolated values enabling generation of gap-free, annual, surface reflectance composites over our time frame of interest (Hermosilla et al. 2015a, 2016). From the time series of surface reflectance pixel composites, we derived Tasseled Cap brightness (TCB), greenness (TCG), and wetness (TCW) (Crist, 1985). Trajectory metrics, typically derived from trend analysis using time series data as slope of the trend, have proven to be useful for characterizing forest structure (Pflugmacher et al. 2012; Frazier et al. 2014). The nature and characteristics of the metrics that describe a given multi-year trend offer insights on successional processes operating, for instance, after disturbances (Bolton et al. 2018). Following this reasoning, we derived temporal trajectory metrics to characterize historic forest dynamics and related trends using the slope of the Tasseled Cap indices from the pre-fire time series, which encompasses from 1984 (first year of the time series) to the year before the year analyzed.

In addition to the annual multispectral surface reflectance composites, the C2C approach also enables the generation of change metrics as a part of the spectral trend analysis (i.e. change magnitude, change persistence, year of greatest change) from Landsat time series (Hermosilla et al. 2016). For example, change magnitude represents the difference in normalized burn ratio (NBR, Key and Benson 1999) before and after a disturbance, while change persistence measures the duration of the disturbance in years. Subsequently, the pixel-level spectral information from the annual composites was combined with temporal and geometry metrics to attribute disturbance objects to a type (e.g. harvest, wildfire, and non-stand replacing disturbance) (Hermosilla et al. 2015b).

Based upon annual Landsat surface reflectance pixel composites and disturbance information, Hermosilla et al. (2018) demonstrated a time series-informed framework for land cover mapping termed Virtual Land Cover Engine (VLCE). VLCE integrates post-classification probabilities, forest disturbance information and forest successional knowledge to mitigate classification errors associated with spurious class transitions through time. The establishment of this time series informed land cover mapping framework allowed for the production of Canada-wide land cover maps with an annual periodicity. Among the 12 land cover classes covered by the VLCE framework, four are treed classes, including wetland-treed, coniferous, broadleaf, and mixedwood.

Forest structure data

Matasci, Hermosilla, Wulder, White, Coops, Hobart, and Zald (2018); Matasci, Hermosilla, Wulder, White, Coops, Hobart, and Bolton et al. (2018) combined annual Landsat surface-reflectance composites and airborne LiDAR samples (Wulder et al. 2012) to seamlessly map forest structure variables annually over Canada’s forested ecosystems using a Nearest Neighbor imputation approach. Herein, we used forest structure information from six modeled structural layers: the mean, standard deviation, and coefficient of variation of canopy heights, the 95th percentile of canopy heights, canopy cover, and gross stem volume.

Topographic and geographic information

Topographic variables were derived from the digital elevation model (DEM) acquired by Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) 30 m version 2 Global Digital Elevation Model (GDEM) product (Tachikawa et al., 2011). Topographic variables included elevation, slope, topographic solar radiation index (TSRI) (Roberts and Cooper, 1989), and topographic wetness index (TWI) (Kirkby and Beven 1979). Finally, latitude and longitude of the pixel centroids (in decimal degrees) were included in this study to account for the association between geographic location and fire disturbance.

Canadian national fire database (CNFDB)

The Canadian National Fire Database (CNFDB) is a spatially explicit database representing wildfire events from 1959 to present time. The CNFDB compiles data from several jurisdictional fire management agencies and from Parks Canada, and replaces the former Large Fire Database (Stocks et al. 2002). The mapping methods used for the generation of the database have varied through time as technology has evolved and new data sources have become available. We used fire year in the CNFDB to inform on wildfires that occurred prior to the start of the C2C time series in 1985 (1959–1984). Information on historical fires allowed us to obtain a more comprehensive understanding of the influence of long-term historical fire regimes on current burn probability over large areas.

Sample selection for model calibration and validation

Figure 1B presents the spatial distribution of forest fires detected between 1985 and 2014 in two time frames (1985 to 1994 and 1995 to 2014), as well as undisturbed treed areas and waterbodies for Saskatchewan. We selected samples (pixels) from burned and undisturbed forests for each year in the time series using several selection criteria. First, only treed pixels (Hermosilla et al. 2018) within Saskatchewan’s forested ecosystems were considered. Second, guided by the forest disturbance dataset (Hermosilla et al. 2015b), the burned samples were selected randomly from areas that were burned in each year, and undisturbed samples were selected randomly from pixels that were not disturbed in that same year, such that areas disturbed by other agents would not be selected (e.g. harvesting). For simplicity, undisturbed samples are referred to as unburned for the remainder of the paper. To ensure representative spectral trajectory metrics of sufficient length (i.e. TCB, TCG, and TCW slopes), a minimum of ten years of observations was deemed necessary; therefore, only fires that occurred between 1995 and 2014 were considered in our model, while samples from 1985 to 1995 were reserved as a temporal baseline for computation of the trajectory metrics. Therefore, the time interval associated with the trajectory variables was at least 10 years, and varies depending on the year each burned or unburned sample is selected for and is representative of.

Following a stratified random sampling strategy, 2500 burned and unburned samples were extracted annually in the time frame of 1995 to 2014, such that the sampling process was not biased by size of fires or numbers of fires in a given year. As a result, 50,000 samples were obtained for burned and unburned classes respectively, representing 0.01% of total burned areas (8.9 Mha) and unburned areas (8.73 Mha) in our study area respectively between 1995 and 2014.

Predictor variables

Predictor variables were derived from Landsat and ancillary data sources to model burn probability (Table 1). Those derived from Landsat time series primarily characterized pre-fire biotic conditions of forests, including land cover, forest structure, and spectral trajectories. The information provided by these layers was extracted for the year prior to each fire occurrence year. We also incorporated topographic and geographic metrics to examine the influence of abiotic factors on burn probability.

Table 1. Predictor variables used to model burn probability.

Characterization domain	Data type	Description	Variable name	Reference
Biotic	Pre-fire land cover (VLCE)	Land cover	Land cover	Hermosilla et al. (2018)
	Pre-fire spectral trajectory metrics	Slope of time series TCB	TCB slope	Bolton et al. (2018)
		Slope of time series TCG	TCG slope
		Slope of time series TCW	TCW slope
	Pre-fire forest structure	Coefficient of variation of vegetation height	Canopy height cv	Matasci, Hermosilla, Wulder, White, Coops, Hobart, and Zald et al. 2018; Matasci, Hermosilla, Wulder, White, Coops, Hobart, and Bolton et al. 2018)
		Average vegetation height	Canopy height mean
		95th percentile of vegetation height	Canopy height p95
		Standard deviation of vegetation height	Canopy height stddev
		Gross stem volume	Gross stem volume
		Percentage of first returns above 2m	CANOPY cover
	Canadian National Fire Database	Fire disturbance year before 1985	CNFDB	Canadian Forest Service (2017)
Abiotic	Geographic coordinates	Latitude	Latitude	–
		Longitude	longitude
	Topography	Elevation	Elevation	GDEM V2 (Tachikawa et al. 2011)
		Topographic wetness index	TWI
		Topographic solar radiation index	TSRI
		Slope	Slope

Burn probability modeling

Predictor variables representing 1994–2013 conditions and corresponding historical trends were used to predict the probability of an area being either burned or unburned in the following year (1995–2014). For burn probability analysis, Random Forests (RF; Breiman 2001) was run in classification mode, where the binary classes of burned and unburned were predicted. As an ensemble learning algorithm, RF consists of a number of decision trees, with final class membership based upon the majority vote from the individual decision trees. These votes are commensurate with class membership likelihoods in the context of soft classification (Du et al. 2012) and therefore can provide a measure of burn probability from 0 to 100% for the subsequent year. Given that decision trees in RF are calibrated with bootstrap samples, approximately one third of the original samples were not used by each decision tree for model calibration, or out-of-bag (OOB) samples. As a result, OOB samples allow for internal cross-validation, without relying on an independent set of validation data (Belgiu and Drăguţ 2016). The output of the internal cross-validation is referred to as OOB error, which enables the calculation of variable importance scores through random permutations. As predictor variables are permuted one at a time, the more important a variable is, a greater decrease will be observed in model accuracy (Breiman 2001).

The developed model was first assessed with OOB samples using confusion matrix based on estimated area proportions, from which overall accuracy and user’s and producer’s accuracies were computed, as well as 95% confidence intervals. The variable importance scores derived therein were used to examine the contribution of biotic and abiotic predictors to burn probability modeling. Partial dependence plots were derived to further evaluate the relationship between these predictor variables and burn probability as captured by the model, in an attempt to uncover the nature between their interactions. Specifically, a partial dependence plot depicts the correspondence between the response variable and a small set of predictor variables, while averaging out the effects of other predictors in a model (Greenwell 2017).

Using the same abiotic variables and biotic predictors representing 2014 conditions, the model was then applied to treed pixels in 2014 to generate a wall-to-wall burn probability (0−100%) map for 2015 at 30 m spatial resolution across Saskatchewan’s forested ecosystems. The burn probability map was then compared to actual burned areas in 2015 that were independently detected using Landsat data (Hermosilla et al. 2016, 2017).

Results

Prediction of burn probability

The RF OOB assessment of the burn probability model was used to evaluate model performance (Table 2). Burned and unburned pixels were discriminated with an overall accuracy of 92.5% and balanced user’s and producer’s accuracies, indicating that the biotic and abiotic predictors were useful for discriminating burned and unburned forest areas for the training data (i.e. 1995–2014; Table 1). Among the biotic factors considered, structural forest descriptors (i.e. canopy cover and mean canopy height) exhibited the highest importance to our burn probability model. Land cover information also ranked among the most important predictors. Abiotic factors describing geographical position (i.e. latitude) and topography (i.e. elevation) were among the top five predictors for modeling burn probability (Figure 3).

Figure 3. Variable importance for the burn probability Random Forests model.

Table 2 Confusion matrix of the estimated proportions of area for the burn probability prediction model for calibration data 1995–2014.

	Reference
	Unburned	Burned	Total	User’s accuracy (%)
Predicted
Unburned	0.421	0.043	0.464	90.7 ± 0.6
Burned	0.032	0.504	0.536	94.0 ± 0.5
Total	0.453	0.547
Producer’s accuracy (%)	92.9 ± 0.4	92.1 ± 0.3
			Overall accuracy: 92.5% ± 0.4%

Biotic and abiotic influence on burn probability

The partial dependence plots in Figures 4 and 5 show the relationship between various biotic factors used to predict burn probability. Overall, forest stands with lower volume and relatively simple structure had higher probability to burn (Figure 4). Wildfires were more prevalent among open forest stands with lower height variability and volume. Interestingly, stands with average height between 3 and 5 m exhibited low burn probability, while those with lower or higher stand height were more susceptible to wildfires (Figure 4A). Results on the influence of land cover on predicting burn probability (Figure 5A) showed that coniferous forests were markedly more likely to burn than wetland-treed, broadleaf, or mixedwood forests. Regarding the influence of spectral trajectory metrics on burn probability, our results indicated that forest areas with declining brightness and wetness trends corresponded to higher burn probability (Figure 5B).

Figure 4. Partial dependence plots for biotic variables on burn probability: (A) canopy height mean and canopy height standard deviation, (B) canopy cover and gross stem volume. The inner ticks along the x-axes represent deciles of the predictor variables.

Figure 5. Partial dependence plots for (A) land cover and (B) temporal trajectory metrics on burn probability, depicting the relationship between burn probability and long-term spectral trends. The inner ticks along the x-axes represent deciles of the predictor variables.

In terms of latitudinal gradient (Figure 6A), northern Saskatchewan is generally more prone to wildfires than southern Saskatchewan, with the high burn probability hotspots mostly located between 55° and 60° N. Although western Saskatchewan exhibited higher burn probability than the eastern region, the longitudinal trend was not as strong as the latitudinal gradient. In terms of topographic drivers (Figure 6B), lowland forests were less likely to burn compared to areas with higher elevation (Figure 6).

Figure 6. Partial dependence plots for (A) geographic and (B) topographic controls on burn probability. The inner ticks along the x-axes represent deciles of the predictor variables.

Wall-to-wall burn probability map assessment in 2015

The wall-to-wall map of predicted burn probability for Saskatchewan’s forested ecosystems in 2015 indicates a marked difference between the low burn probability values found in southern (below latitude 56°N) compared with northern areas (Figure 7A). Similarly, the western forests of the province presented higher burn probability values compared with eastern forests. Areas predicted to have a high burn probability exhibited strong spatial agreement with actual wildfires that occurred in 2015 at a range of spatial scales, due to the capacity of our 30-m map to capture fine spatial details of burn probability (Figure 7). Compared to the relatively compact burned patches in 2015, the predicted probability map is less spatially continuous at a local scale (Figure 7C); however, the burned area overall exhibited higher burn probability than surrounding unburned forest stands. Indeed, pixels that burned in 2015 had a much higher burn probability (mean = 69.7%), which was markedly skewed toward the highest burn probability values (median = 81%) compared to those pixels that did not burn in 2015 (mean = 27.3%), which were also markedly skewed toward the lowest probabilities (median = 19%) (Figure 8). Only 6.8% of areas that burned presented burn probabilities lower than 19%. Overall, 71.5% of area burned had burn probabilities that exceeded 50%.

Figure 7. (A) Predicted forest burn probability in 2015 in Saskatchewan’s forested ecosystems, in comparison to actual burned forest areas for 2015, as detected independently by the C2C approach. Zoom-in examples showing spatially detailed agreement between predicted burn probability and detected fires in (B) Boreal Shield West and (C) Boreal Plains.

Figure 8. Distribution of predicted burn probability values for burned and unburned treed areas in our study area detected from Landsat time series data. Vertical solid lines represent the median probability for the burned and unburned class. The vertical dotted line represents the median probability for the entire study area.

Discussion

Biotic and abiotic influence on burn probability

We assessed the driving factors (biotic and abiotic) on burn probability of Saskatchewan’s forested ecosystems using time series of Landsat data. Our objective was to determine the importance of these driving factors, identify the nature of their relative influence on fire occurrence, and map burn probability. Based on fires that occurred between 1995 and 2014, our results suggested that biotic predictors describing forest structure are key to modeling burn probability. This aligns with recent findings identifying the availability of fuels as the most significant driving force behind high-severity fires across western United States (Parks et al. 2018). Herein we found that canopy cover was the most important predictor (Figure 3), exhibiting a strong negative correlation with burn probability (Figure 4B), so areas with high canopy cover had a lower probability to burn. Considering that taller forest stands exhibited a higher burn probability (Figure 7), burn probability could correspond with forest succession (Whitman et al. 2019). Thus, as forest stands progress from stand initiation to stem exclusion and the understory re-initiation stage (Oliver et al. 1996), stand height generally increases while competition among dominant trees might lead to increased canopy openness. As a result, development in the understory and overall vertical distribution would result in higher loading in forest floor fuels (de Groot et al. 2009) and ladder fuels (Menning and Stephens 2007), thus increasing burn probability. Erni et al. (2018) demonstrated that stand age could impose significant influence on the spatial patterns and magnitude of burn probability due to the stand age resistance effect, which illustrates that young naturally regenerated stands are more resistant to wildfire than mature stands owing to lack of fuels. Although we incorporated the historical information of wildfires provided by the CNFDB as a proxy of stand age into our model of burn probability, this information did not contribute strongly to our burn probability predictions (Figure 3). This may suggest that forest structural metrics themselves provide an adequate description regarding forest age for the identification of fire prone stands. In addition to forest structure, land cover and information on the pre-disturbance spectral trends also contributed to the burn probability model (Figure 3). Our results indicated that coniferous stands were more prone to burn compared to other stand types, concurring with the trends identified by the Canadian Fire Danger Rating System (CFDRS; Stocks et al. 1989) and the land cover investigation of Castilla et al. (2016). As we did not incorporate climate or weather data into our burn probability model, it is challenging to determine the magnitude of influence that can be uniquely attributed to these biotic variables, apart from what may be correlated with climatic patterns.

Pre-fire vegetation dynamics were characterized in this study by analyzing the historical spectral trends on wetness and brightness observed by Landsat time series. These spectral trajectory predictors exerted a moderate influence on burn probability. In particular, forest stands with declining trends in brightness and wetness over time exhibited higher burn probability (Figure 5), potentially indicating the influence of non-stand replacing disturbances (e.g., drought, defoliating, insects, vegetation decline; Cohen et al. 2016; Bell et al. 2018), which have been found to precede wildfire events (Hermosilla et al. 2019).

Among the abiotic factors analyzed, latitude and elevation were ranked within the most important predictors for modeling burn probability (Figure 3). In particular, forest stands above 600 m were much more likely to burn than lowland forests. This finding agrees with the relationships between burn severity and elevation found by Parks et al. (2014) in wildfires in New Mexico and Idaho. Geographically, northwestern Saskatchewan was more likely to be disturbed by fire than southeastern Saskatchewan (Figure 4). Saskatchewan’s managed forest areas is found primarily in the southern extent of the province’s forested ecosystems. These managed forest areas comprise forest tenures for harvesting (Stinson et al. 2011) and have intensive fire suppression practices compared to the unmanaged forest areas to the north (White et al. 2017), resulting in a suite of landscape characteristics with a lower probability of burning. Overall, the results indicate that the fine-scale forest structural information enabled by Landsat time series can improve our understanding of the relationship between various structural characteristics and susceptibility to wildfire.

Burn probability mapping at a 30-m spatial resolution and implications for Forest management

In this study we built a model to predict forest burn probability in 2015, calibrated using burned and unburned pixels from the period 1995 to 2014. RF votes were used as class membership likelihoods to characterize burn probability in Saskatchewan’s forested ecosystems (Figure 7A). We found strong spatial agreement between wildfires that actually occurred in 2015 and area with high burn probability (Figures 7 and 8). When assessing spatial agreement at a local scale, it is worth noting that not all pixels in a burned stand in 2015 exhibited high burn probability based on their pre-fire vegetation and topographic characteristics. Accordingly, unburned stands might include patches with relatively high predicted burn probability, which is common in the context of burn probability modeling. As highlighted by Parisien, Ager et al. (2020), the stochastic nature of wildfire occurrence leads to fires burning areas predicted with low probability by mere chance.

While burn probability predictions are able to capture fine scale variability in fuel conditions with a spatial resolution of 30 m, the pixel-based model did not take the spatial context, such as fire patch perimeter or area, into consideration in the probability model (Parisien et al. 2013; Erni et al. 2018; Whitman et al. 2018). For example, a patch with high fuel loadings would have lower burn probability if it was surrounded by stands with low fuel loadings due to the fire shadow effect (Parisien et al. 2011; Erni et al. 2018). Additionally, many other factors controlling forest ignition and burn probability (e.g. lighting strike, seasonal precipitation, and weather) remain stochastic in nature, which are not accounted for in our model. Nonetheless, more than 71.5% of area burned in 2015 had burn probability values exceeding 50%. To further improve model performance, annual climatic variables could be incorporated to capture factors such as ignition probability in a spatially explicit manner (Scott et al. 2012).

Compared to simulation-based burn probability models, such as FlamMap (Finney 2006), FSim (Finney et al. 2011), and Burn-P3 (Parisien et al. 2005), our model shares the functionality of generating burn probability maps, yet differs from these simulation-based models in a number of key ways. In terms of overall model design, simulation-based burn probability models generally follow a mechanistic approach to model fire growth through fire ignition and spread, while accounting for effects of topography, fuels, and fire weather (Parisien, Dawe et al. 2020). In contrast, we aim to establish the connection between pre-fire biotic and abiotic factors and burned and unburned samples mapped by time-series Landsat spectral data following a statistical approach, which generates burn susceptibility estimates indirectly. While both approaches rely heavily on input data, especially those that represent vegetation conditions, fuel grids used in existing burn probability models are typically converted from land cover data, which may not be sufficiently detailed, up-to-date, or accurate (Parisien et al., 2013). In comparison, the forest structural layers used in our model are derived at an annual time step using airborne LiDAR transect data (Matasci, Hermosilla, Wulder, White, Coops, Hobart, and Zald et al. 2018; Matasci, Hermosilla, Wulder, White, Coops, Hobart, and Bolton et al. 2018) covering a range of forest structural attributes, offering a capability to characterize vegetation conditions in a timely manner. As a result, the statistical burn probability model presented in this research estimates burn probability for a specific year based on forest conditions, as highlighted with the 2015 prediction. In contrast, simulation-based models operate in a temporally holistic fashion over numerous iterations and calculate burn probability as the number of times a cell burned over the total number of iterations (Parisien, Dawe, et al. 2020). This difference in model output is evident from the range of predicted burn probability from our model (0–100%), compared to that of simulation-based burn probability outputs (e.g. Burn-P3), which tends to be much narrower and skewed to lower values (e.g. Parisien et al. 2005; Beverly and McLoughlin 2019). Therefore, burn probability derived from our model should be interpreted as a relative susceptibility to wildfire (Davis et al. 2017), thus requiring additional normalization prior to any direct comparison with burn probability estimates from simulation-based models. Lastly, due to the mechanistic nature of simulation-based models, they offer users the flexibility to explore burn probability and fire behavior under alternative scenarios, such as fuel treatment (Finney 2007), wildlife conservation (Ager et al. 2010), and climate projection (Wang et al. 2016). As ours is a statistical-based model, it is less suitable for answering such “what if” questions. Even though modifications to the input variables can be made to mimic effects such as fuel treatment, it is not currently designed for exploring alternative scenarios. It is evident that there are diverse approaches to gather insights on where wildfire may occur, all of which have strengths and weaknesses. Given the potential effects of wildfires, any relevant information should be considered (Parisien, Ager, et al. 2020).

The spatial agreement between areas with high predicted burn probabilities and realized fire locations in 2015 suggests that burn probability maps could be used as additional input information in wildfire risk models (Taylor et al. 2013, Xi et al. 2019) and decision support tools. Burn probability maps can contribute meaningfully to strategic fire and fuel management decisions, including fire risk evaluation (Finney et al. 2011; Alcasena et al. 2017), fuel treatment optimization (Salis et al. 2018), and prevention planning (Miller et al. 2008). The continuous nature of burn probability outputs is useful for both managers and modelers. Forest managers and decision makers can establish priority areas of interest by determining thresholds of burn probability that meet their information needs, to prioritize resources or harden at-risk communities, while also being mindful of omission and commission error (Figure 9). From a forest management perspective, the selection of the threshold value would have implications on extent of the area of interest, and hence on the material and economic requirements needed to undertake the prevention, treatment, or monitoring of areas with high probability of burning. In many ways, the use of the full range of probability values remains of greatest utility, as these provide more detailed information, while avoiding possibly arbitrary cutoffs. As indicated in Figure 9, overall agreement between predicted burn probability and actual burned areas in 2015 increases as the burn probability increases. Overall, approximately 21% of Saskatchewan’s forested ecosystems had a burn probability greater than 50% in 2015. From our analysis herein, we determined that 75% of this area actually did burn in 2015 and captured the majority of the actual fires (burned area underestimation = 28.8%) but with an evident overestimation regarding the burned area (fire overestimation 83.1%). Liu et al. (2013) found that the application of a relatively high burn probability threshold was useful for identifying focus areas for fuel treatment, which in turn can aid in reducing the number of fires and burned area. Similar opportunities exist to examine the potential influence of wildfire on urban interfaces (Massada et al. 2009; Salis et al. 2013) and wildlife conservation (Thompson et al. 2013; Bentley and Penman 2017).

Figure 9. (A) Predicted burned area, and (B) predictive agreement with actual fires in 2015 as function of the threshold applied to burn probability to predict burned areas.

Conclusion

In this research, we presented a data-driven approach to predict burn probability across forested ecosystems of Saskatchewan, Canada, using a suite of Landsat time series based forest structural, spectral, and land cover metrics, together with ancillary geographic and topographic information. Herein we generated a burn probability model and further analyzed the relationship between these predictive metrics representing biotic and abiotic factors and the model outcome. Subsequently, we applied the burn probability model to data relating 2014 conditions to map burn probability in 2015, which exhibited strong correspondence with actual fires occurring that year captured from an independent source. Among burned forests in 2015, only 6.8% had burn probabilities that were lower than 19%, in contrast, 71.5% had burn probabilities greater than 50%. Descriptors of forest structure were the most useful predictors, highlighting the importance of data sources that are capable of characterizing forest structure over large areas with sufficient accuracy, spatial and temporal resolution. Short and open stands with low height variability and volume exhibited relatively high burn probability, understood as a function of the forest characteristics present in the northern extent of the study area that is also most prone to wildfire. Among various cover types, coniferous stands were the most likely to burn, while broadleaf forests were the least likely to burn. From a stand development perspective characterized by temporal trajectory metrics, stands with significant negative long-term trajectories, likely under the influence of non-stand replacing disturbances, were the most susceptible to wildfires. In terms of topography, stands located on gentle slopes and higher elevations tended to be more susceptible to wildfires. The approach presented herein demonstrates how detailed, spatially-explicit burn probability maps may be generated from information characterizing forest structure and past disturbances derived from Landsat time series. Value is demonstrated in that these burn probabilities are well aligned in a predictive sense with areas that actually burned, representing a new and useful source of information for wildland fire risk modeling, management, and mitigation.

Acknowledgements

This research was undertaken as part of the “Earth Observation to Inform Canada’s Climate Change Agenda (EO3C)” project jointly funded by the Canadian Space Agency (CSA), Government Related Initiatives Program (GRIP), and the Canadian Forest Service (CFS) of Natural Resources Canada (NRCan). Support has also been provided by a NSERC Discovery grant to Coops. This research was enabled in part by computational resources provided by WestGrid (www.westgrid.ca) and Compute Canada (www.computecanada.ca). Open access is supported by the Government of Canada. Geordie Hobart and Dr. Giona Matasci are thanked for their contributions to generating the image composites and the forest structural layers used herein.