Exceptionally large areas burned in 2014 in central Northwest Territories (Canada), leading members of the Tłchǫ First Nation to characterize this year as ‘extreme’. Top-down climatic and bottom-up environmental drivers of fire behavior and areas burned in the boreal forest are relatively well understood, but not the drivers of extreme wildfire years (EWY). We investigated the temporal and spatial distributions of fire regime components (fire occurrence, size, cause, fire season length) on the Tłchǫ territory from 1965 to 2019. We used BioSIM and data from weather stations to interpolate mean weather conditions, fuel moisture content and fire-weather indices for each fire season, and we described the environmental characteristics of burned areas. We identified and characterized EWY, i.e., years exceeding the 80th percentile of annual area burned for the study period. Temperature and fuel moisture were the main drivers of areas burned. Nine EWY occurred from 1965 to 2019, including 2014. Compared to non-EWY, EWY had significantly higher mean temperature (>14.7°C) and exceeded threshold values of Drought Code (>514), Initial Spread Index (>7), and Fire Weather Index (>19). Our results will help limit the effects of EWY on human safety, health and Indigenous livelihoods and lifestyles.
Des superficies exceptionnellement vastes ont brûlé en 2014 dans le centre des Territoires du Nord-Ouest (Canada), ce qui a conduit les membres de la Première Nation Tłchǫ à qualifier cette année d'« extrême». Les facteurs climatiques descendants et les facteurs environnementaux ascendants expliquant le comportement du feu et les superficies brûlées en forêt boréale sont relativement bien connus, mais pas les facteurs associés aux années de feu extrêmes (EWY). Nous avons étudié les distributions temporelles et spatiales des composantes du régime des feux (occurrence, taille, cause, durée de la saison de feu) sur le territoire Tłchǫ de 1965 à 2019. Nous avons utilisé BioSIM et des données de stations météo pour interpoler les conditions météorologiques moyennes, la teneur en humidité du combustible et des indices forêt-météo pour chaque saison de feu et nous avons décrit les caractéristiques environnementales des zones brûlées. Nous avons identifié et caractérisé les EWY, c.à-d., les années dépassant le 80e percentile des superficies annuelles brûlées de la période d'étude. La température et la teneur en humidité du combustible étaient les principaux déterminants des superficies brûlées. Neuf EWY sont survenues entre 1965 et 2019, incluant 2014. Les EWY se différencient des autres années par une température moyenne significativement plus élevée (> 14,7°C) et par le dépassement de valeurs seuils de l'indice de sécheresse (> 514), de l'indice de propagation initiale (> 7) et de l'indice forêt-météo (> 19). Nos résultats aideront à limiter les effets négatifs des EWY sur la sécurité et la santé des personnes ainsi que sur le maintien des moyens de subsistance et du mode de vie autochtones.
Introduction
Fire regimes have become more active in Canada's boreal forests over the past 50 years (Kasischke and Turetsky 2006; Hanes et al. 2019). In northwestern Canada, rising temperatures associated with climate change have resulted in shorter fire intervals, larger areas burned, and more severe fire effects on ecosystems (Whitman et al. 2018, 2019). For example, in 2014, more than 3.4 million hectares burned in the Northwest Territories (NT), compared to 600,000 hectares burned on average each year from 2009 to 2019 (NTENR 2015; Supplemental Material S1 (teco_a_2070342_sm6945.docx)). In central NT, the 2014 fires affected more than 20% of the territory of the Tłchǫ First Nation (∼800,000 ha), leading community members to characterize this year as ‘extreme’, both in terms of magnitude and consequences of the fires (Morarin 2020).
Large and very large fires, from several hundred to several thousand hectares, respectively, can alter ecosystem ecological functions, increase the risk for human safety and socioeconomic interests, and threaten human communities and values (Moritz et al. 2014), including cultural and subsistence activities of Indigenous peoples (Adams 2013; Morarin 2020). Large fire sizes may reduce fire return intervals in some areas, thereby affecting species composition, succession and stand structure (Bergeron 1998; Johnstone and Chapin 2006; Soja et al. 2007; Whitman et al. 2018; Baltzer et al. 2021). Shorter intervals between fires can reduce ecosystem resilience by limiting tree regeneration (Asselin et al. 2006; Johnstone and Chapin 2006; Burton et al. 2008), and alter the long-term net carbon balance of boreal forests (Walker et al. 2019). Alteration of the forest mosaic by large fires can also modify surface albedo with consequences for soils and climate, altering the energy balance between the atmosphere and the Earth's surface (Schimel and Baker 2002; Bond-Lamberty et al. 2007; Bonan 2008; Oris et al. 2014). By altering the physical landscape, large fires affect wildlife and their habitat (Barrier and Johnson 2012; Anderson and Johnson 2014). Some studies have also highlighted destruction of infrastructure due to large fires (Williams 2013; Gauthier et al. 2015), as well as negative effects on human health through plumes of smoke transporting fine particles, which decrease air quality and cause respiratory diseases in local populations (Kochtubajda et al. 2016; Dodd et al. 2018). Large fires limit the capacity of Indigenous people to engage in cultural and subsistence activities on the land (Morarin 2020), reducing their access to ecosystem services (Berkes and Davidson-Hunt 2006; Royer and Herrmann 2013; Bélisle and Asselin 2021). These fires can jeopardise some cultural keystone places crucial to land-based Indigenous livelihoods and lifestyles (Woo et al. 2007; Christianson 2011; Legat 2012; Cuerrier et al. 2015).
Whatever their size, wildfires are natural disturbances that contribute to forest ecosystem processes and determine forest community structure, health, productivity and diversity (Weber and Stocks 1998; Wright and Heinselman 2014; Guindon et al. 2018). Large and severe crown fires are thus part of the natural functioning of boreal forests, although they are relatively infrequent (Stocks et al. 2002; Boonstra et al. 2016). In the NT, the year 2014 was notable, both because of the number of fire ignitions related to lightning activity (Kochtubajda et al. 2016; Veraverbeke et al. 2017) and the large areas burned covering several hundred thousand hectares. At northern latitudes, fires are commonly left to burn due to low population density, few industrial activities and slow response time (Tymstra et al. 2020). Hence, when weather, fuel moisture content and environmental conditions are favourable, coalescing intense crown fires spread rapidly, sometimes reaching 50 m.min–1 (Damoah et al. 2006), which makes their control difficult (Peters et al. 2004; Stephens et al. 2014; Erni et al. 2020). In the NT, only high priority areas are protected, i.e., residential, resource-rich and culturally significant areas (Government of the Northwest Territories 2016).
Fire regimes describe ‘the characteristic patterns of wildfires over large spatial and temporal scales and, they are sensitive to changes in climate, vegetation, and ignitions’ (Higuera 2015, p. 13,137). The spatial components of a fire regime include fire type, extent, and severity, while the temporal components describe the recurrence of fire at various timescales (frequency, return interval, seasonality, season length). Fire regimes are controlled by the fire environment, notably fuel type and abundance, which affects fire behavior and determines fire effects through time (Crutzen and Goldammer 1993; Morgan et al. 2001). In northern Canadian boreal forests, fire behavior is influenced by ‘top-down’ drivers such as weather and climate, and ‘bottom-up’ drivers such as topography, soil and land cover type (Whelan 2006; Falk et al. 2007; Gaboriau et al. 2020). Top-down drivers determine weather and climate conditions, fuel moisture and flammability at hourly, seasonal, and successional timescales (Macias Fauria et al. 2011). Bottom-up drivers are landscape characteristics corresponding to elevation, slope, fuel type and availability, and land cover (Mansuy et al. 2014; Walker et al. 2020). The interaction of top-down and bottom-up drivers influences fire occurrence and spread. The main top-down drivers of northern boreal fire regimes are fairly well documented (Macias Fauria and Johnson 2006, 2008). Previous work has shown that exceptionally high summer temperatures and dry conditions increase the fire season length and the number of days of potential fire spread (Jain et al. 2017), in turn increasing the number of large fires. The effects of bottom-up drivers such as snow cover (Euskirchen et al. 2016), permafrost thawing (Schuur et al. 2008), topography, land cover and vegetation changes (Cumming 2001; Beck et al. 2011; Myers-Smith et al. 2011) have also been reported. However, while long-term processes related to fire behavior and area burned in the boreal forest are relatively well understood, the drivers of large fires and extreme wildfire years (i.e., years with the largest area burned in a study area and time period; hereafter EWY) are still underreported.
We aimed to reconstruct the contemporary wildfire regime in central NT. For that, we characterized some temporal and spatial fire regime components (fire ignition cause, fire occurrence and size, fire season length) on the territory of the Tłchǫ First Nation (NT) from 1965 to 2019. We explored the associations between climate, environmental conditions and fire regime components using multivariate methods and statistical modelling. In addition, we identified years with the largest annual area burned to characterize EWY over the study period. We discriminated EWY and non-EWY based on seasonal weather, fuel moisture and fire-weather conditions, hence establishing thresholds associated with EWY occurrence. In view of previous studies on crown fires and their drivers, we hypothesized that seasonal warm temperature and drought (Balshi et al. 2009; Parisien et al. 2011; Whitman et al. 2019), as well as fuel availability and forest species composition (described as land cover types), would be the main drivers determining the annual area burned (Bernier et al. 2016; Thompson et al. 2017; Gaboriau et al. 2020). Our findings could inform management actions to limit the negative effects of large fires on cultural keystone places and ecosystem services that benefit the Tłchǫ people and neighboring First Nations.
Material and methods
Study area
The study area covers 39,400 km2 and corresponds to the Tłchǫ First Nation territory, Northwest Territories, Canada (Figure 1a). In 2003, an agreement was signed between the Tłchǫ First Nation and the Canadian and NT governments to recognize Tłchǫ ownership and management of land and resources within the territory (Tłchǫ Government 2003). As a result, timber harvesting in the study area is subject to strict regulations to conserve Tłchǫ heritage, identity and traditional land uses, for example by maintaining ancestral trails, caribou migration corridors, and watersheds (Government of the Northwest Territories 2016). The four Tłchǫ communities have a total population of 2,901 members (Northwest Territories Bureau of Statistics 2021).
The elevation in the study area ranges from 156 meters above the sea level (masl) south of Behchokǫ̀ to 479 masl east of Wekweètì (Figure 1b; Natural Resources Canada 2016). Climate is continental and dry with long cold winters and short warm summers (Beck et al. 2018). The Yellowknife weather station south of the study area (62°27′ N, 114°26′ W, 206 masl; Environment Canada 2021) recorded a mean annual temperature of –4.3°C during the period 1981–2010, with mean monthly temperatures of 17°C and –25.6°C in July and January, respectively. The mean annual precipitation averaged 289 mm, with 170 mm (59%) occurring as rain from April to October (Environment Canada 2021). The western part of the study area is in the Taiga Plains ecoregion, which is warmer and wetter (Figures 1c and 1d) than the eastern part located in the Taiga Boreal Shield ecoregion (Figure 1e). The Taiga Plains are characterized by a thick layer of glacial sediments, and their soils are mostly gelic histosols, while the Taiga Boreal Shield is characterized by a Precambrian bedrock with thin dystric cambisols (Ecological Stratification Working Group 1996).
The study area includes many lakes formed after glacier retreat ca. 10,000–9,000 years before present (BP) (Dyke 2005; Latifovic et al. 2017) and large peatlands dominated by sphagnum mosses established after 6,000 BP (MacDonald 1995). Except for the northeast, which is mainly composed of grassland, lichen and moss barrens, sub-polar and taiga needleleaf forests dominate the landscape (Figure 1f). Forest density is higher in the southwest than in the northeast, following elevation, temperature and precipitation gradients (Figures 1b, 1c and 1d). Within the study area, black spruce (Picea mariana (Mill.) B.S.P.) dominates forest composition, whereas jack pine (Pinus banksiana Lamb.) is more frequent in the south where tree species diversity, temperature and fires are higher (Beaudoin et al. 2014). Broadleaf tree species, such as trembling aspen (Populus tremuloides Michx.) and paper birch (Betula papyrifera Marsh.), are commonly found in the southern part of the study area and near wetlands (Beaudoin et al. 2014).
Fire data
We used the Canadian National Fire Database (CNFD; Natural Resources Canada 2021) to extract information on fires having occurred from 1965 to 2019 on the Tłchǫ First Nation territory. We removed very small fires (<1 ha) that are usually human-caused and located near settlements (Stocks et al. 2002). We then calculated the number of fires per year or annual fire occurrence (AFO), annual area burned (AAB), and the fire size distribution (Table 1). We log-transformed the AAB data set (LogAAB) to meet the normality assumption, and we used Pearson's Rho test to measure the relationship between AFO and LogAAB. We measured the direction, strength and significance in time series of AFO and AAB from 1965 to 2019 by conducting a Mann–Kendall test using the ‘Kendall’ package in R (McLeod 2015).
Figure 1.
(A) Location of the study area encompassing the four communities of the Tłchǫ First Nation in the Northwest Territories, Canada (GA = Gamèti, WH = Whatì, BE = Behchokǫ̀, WE = Wekweèti). (B) Elevation at 10-m spatial resolution (in meters). (C) 1965–2019 mean temperature gradient (southwest to northeast) at noon (degrees Celsius) for the fire season period extracted from four weather stations and interpolated with BioSIM in 10-km grid cells and (D) 1965–2019 sum of precipitation (mm) gradient (northeast to southwest) for the fire season period interpolated with BioSIM. (E) Soil-based ecoregion classification in 10-km grid cells. (F) Main land cover types in 2010 per 10-km grid cell.

Climate data and fire season length
To estimate average climate and fire-weather conditions during the fire season from 1965 to 2019, we extracted daily weather data from the four weather stations located nearest the study area (Lac La Martre, Lower Carp Lake, Rae Lakes, Yellowknife Hydro; between 62° 42′ N, 113°54′ W, and 64°05′ N, 117°17′ W; Figure 1a) using the BioSIM software (Régnière and Saint-Amant 2014). For each weather station location, we extracted daily weather data including temperature (°C), wind speed (km h–1), relative humidity (%) and total daily precipitation (mm, as water equivalent) measured at noon. We then used the BioSIM software to produce time series of daily weather data for each centroid of 10-km grid cells covering the study area (n = 475 grid cells). BioSIM interpolated data from the four closest weather stations using inverse distance weighting output, while adjusting for differences in elevation and location differentials with regional gradients (Régnière and Saint-Amant 2014). Next, we used interpolated daily weather data as input to the Canadian Fire Weather Index (FWI) daily model implemented within the BioSIM software for calculation of the daily fuel moisture content and fire-weather conditions for each centroid. The FWI is an internationally used fire danger prediction system based on daily weather observations affecting the potential of fires to ignite and spread (Van Wagner 1987). Fuel moisture content is described by the Fine Fuel Moisture Code (FFMC), the Duff Moisture Code (DMC) and the Drought Code (DC) (Table 1). The FFMC, DMC, and DC indicate the moisture content in surface fine fuels, organic and deep organic layers, respectively, with higher values of these indices representing lower moisture content (Van Wagner 1987; Wotton 2009). The other three components are indices of fire behavior, and their values are proportional to fire danger (Table 1): rate of fire spread (Initial Spread Index; ISI), fuel availability (Buildup Index; BUI), and fire intensity at the fire front (FWI) (Van Wagner 1987; Wotton 2009).
Table 1.
Characteristics of the fire regime components, mean weather conditions, and averages of fuel moisture content and fire-weather indices in the study area for the fire seasons from 1965 to 2019.

For each year and each centroid, we averaged the daily weather data, daily fuel moisture content and fire-weather indices, and summed daily precipitation over the fire season length, which we estimated following the standard protocols applied by fire managers in Canada to estimate the start and end of each fire season. The annual time series of climate variables and FWI indices respect the assumption of normality. We assumed that the fire season started on the third day after snowmelt or on the third consecutive day with noon temperature above 12°C, whichever occurred first (Lawson and Armitage 2008). We also considered the overwintering effect that is representative of the previous year's legacy index values at the end of the fire season instead of a simple ‘reset’ of indices computation. Similarly, we assumed the fire season ended with the accumulation of 2 cm of snow on the ground for seven consecutive days or after three consecutive days with daily minimum temperature below 0°C, whichever occurred first (Lawson and Armitage 2008). To assess the accuracy of our fire season length estimation with regard to snowmelt, we compared fire season length estimates with snowmelt timing extracted from Moderate Resolution Imaging Spectroradiometer (MODIS) images for 2001–2015, after having resampled grid cell sizes from 250-m to 10-km resolution. We found a mean difference of 22 days between data sets, with MODIS images showing a later snowmelt than fire season length estimates. However, we also found that data sets were significantly and positively correlated (Pearson’s r = 0.72, p = 0.002; Supplemental Material S2 (teco_a_2070342_sm6945.docx)). We measured the direction, strength and significance of temporal trends in the fire season length and climate variables over the 1965–2019 period, using Mann-Kendall tests.
Fire-climate associations
First, we linearly detrended climatic variables that showed significant linear trends according to Mann-Kendall tests, by extracting the residuals from linear regressions. Next, we measured the first-order autocorrelation on the detrended series. The AR1 (autoregressive model of order 1) for fire season length and climatic variables was non-significant. Then, we examined the associations (Pearson's r correlations) between two fire regime components (annual fire occurrence (AFO) and log-transformed annual area burned (LogAAB)), and annual averages for climate conditions (sum for annual precipitation), fuel moisture content and fire-weather indices during fire seasons from 1965 to 2019.
We applied two modelling approaches based on Generalized Additive Mixed Models (GAMM; Wood 2017) and Multivariate Adaptive Regression Splines (MARS; Friedman 1991). The comparison of GAMM and MARS could provide insights into their relative strengths. GAMMs are an extension of generalized additive models that use smoothing functions that allow flexible description of complex responses to predictor relationships. GAMMs require normality of residuals. On the other hand, MARS combines the strengths of regression trees and spline fitting by replacing the step functions normally associated with regression trees with piecewise linear basis functions. Both methods allow modelling of complex non-linear relationships between a response variable and its predictors. However, the computational complexity of GAMM makes cumbersome the generation of predictions for independent data sets. The simple rule-based basis functions in MARS greatly facilitate prediction.
We tested for multicollinearity of explanatory variables by calculating variance-inflation factors (VIF) for linear models using the ‘car’ and ‘olsrr’ packages in R (Hebbali 2017; Fox and Weisberg 2019). For a given predictor (p), the VIF measures how much the variance of a regression coefficient is inflated due to multicollinearity in the model. We removed, one by one, variables that yielded VIF > 5. By retaining a single variable from a subset of correlated variables, we assume that some drivers explaining fire activity could be excluded, despite their strength as an explanatory variable. We determined that seven of the explanatory variables, representing climate conditions and fuel moisture content, were non-collinear, namely fire season averages of temperature, relative humidity, wind speed, DMC, DC, ISI and the sum of precipitation during the fire season. For both GAMM and MARS approaches, we assumed that these seven variables had the same weight and we compared explanatory variables selected by each model, looking at generalized cross-validation statistics. More information about the parameters used for each model is provided in Supplemental Material S3 (teco_a_2070342_sm6945.docx). To determine the best model to predict fire metrics (AFO and LogAAB) through time, we calculated the Pearson correlation between observed and predicted values by each model from 1965 to 2019. The null hypothesis of no significant correlation was rejected when p < 0.05.
Associations between elevation, land cover types and fires
We examined the spatial associations between elevation, land cover types, and fire regime components during the study period. First, we summed the number of fire ignitions and area burned [log(total area burned by cell) +1] between 1965 and 2019 in each of the 475 cells of 10-km resolution covering the study area. For each grid cell, we extracted the mean elevation from a 10-m Canadian Digital Elevation Model (Figure 1b; Natural Resources Canada 2016) and we resampled the data, calculating the mean at 10-km resolution. We then measured the spatial correlations between elevation and fire regime components in each cell covering the study area for the total study period.
We extracted the proportion of land cover types in 2010 at a 30-m resolution and resampled it for each 10-km grid cell within the study area using Landsat data from the Canada Centre for Remote Sensing (NALCMS 2017) (Figure 1f). Seven land cover types were distinguished: sub-polar needleleaf forest, taiga needleleaf forest, sub-polar shrubland, sub-polar grassland, grassland-lichen-moss, barren lands, and water ( Supplemental Material S4 (teco_a_2070342_sm6945.docx)). We restricted the period of analysis for associations with land cover types (2010–2019) because we did not have information before 2010, and land cover types likely changed since 1965, particularly as a result of fire. To assess spatial associations, we used the ‘SpatialPackv0.3–8’ package in R (Osorio et al. 2019) and conducted a modified version of the t-test, which accounts for the spatial structure (latitude and longitude) of predictor variables (Dutilleul et al. 1993) and prevents type I error inflation.
Detection of extreme wildfire years
As suggested by the definition of extreme events (IPCC 2012), we defined extreme wildfire years (EWY) as years at the upper end of the AAB data set from 1965 to 2019. We applied Bartlett's goodness-of-fit test using the ‘Renext’ package in R (Deville et al. 2016) to confirm the exponential distribution of AAB (p < 0.001). We estimated the parameters of the AAB data set, retaining only fire years (i.e., 48 out of 55 years), by fitting a model including the 20th-80th percentile range corresponding to the bulk of observations. We used the ‘getOutliers’ function in the ‘extremevalues’ package in R (Van der Loo 2020) to identify the residuals of the estimated AAB data set higher than the 80th percentile, i.e., years within the highest 20% of the distribution and with a low probability of occurrence (p = 0.01 confidence limit; Supplemental Material S5 (teco_a_2070342_sm6945.docx)) that we called EWY.
Top-down and bottom-up drivers of extreme wildfire years
We applied Wilcoxon-Mann-Whitney tests to assess significant differences between extreme wildfire years (EWY; n = 9) and non-EWY (n = 46) based on their average climate conditions during the fire season and the length and start date of the fire season from 1965 to 2019. We applied the same tests to assess significant differences between EWY and non-EWY based on fuel moisture content and fire-weather indices from 1965 to 2019. For each index, we counted the number of days during the fire season when the index value was ≥75th percentile of all daily data of the 1965–2019 period.
Figure 2.
(A) Annual fire occurrence (blue line CNFD – points) and annual area burned (bars – CNFD polygons) on the territory of the Tłchǫ First Nation over the 1965–2019 period, with extreme wildfire years (EWY) highlighted in red. Ten-year intervals are separated by black dotted lines. The number of asterisks represents the number of fires > 20,000 ha recorded each year in the study area. (B) Monthly total area burned (grey bars – CNFD polygons) and total fire occurrences (blue line – CNFD points) during the fire season from 1965 to 2019 on the Tłchǫ First Nation territory.

To test for significant differences between EWY and non-EWY based on elevation and land cover types, we applied Wilcoxon-Mann-Whitney tests. We distinguished the mean elevation of fire perimeters for 103 fires >1 ha that burned during EWY, and 229 fires >1 ha that burned during non-EWY from 1965 to 2019. We distinguished the land cover types burned for 23 fires >1 ha that burned during EWY and 28 fires >1 ha that burned during non-EWY from 2010 to 2019.
Results
Fire regime components
In all, 752 fires >1 ha have burned ∼ 18,000 km2 of forests or 45% of the Tłchǫ territory from 1965 to 2019. On average 14 fires (AFO) and 326 km2 of forests burned (AAB) per year (sd = 10 and 1111 km2; Figure 2a). During the single 2014 year, 27 fires burned 8,032 km2 of forests in the study area (Figure 2a). No significant trend in AFO and AAB was observed during the study period ( Supplemental Material S6 (teco_a_2070342_sm6945.docx)). AFO and LogAAB were positively correlated (Pearson's r= 0.67, p < 0.001, n = 55), but time series had interannual variability over the past five decades (Figure 2a). The 1970s and 1990s had the highest AFO of the study period, and the 1970s and 2010s recorded largest AAB (Figure 2a). Although the number of fires per year was important during the 1990s, the annual area burned was relatively small, implying many small fires, whereas during the 2010s fires were less numerous but larger. The fire season mainly spans from June to August, and the peak in fire occurrence and area burned occurs in June and July (respectively, 78% of the total number of fires and 97% of the total area burned during both months; Figure 2b). From 1965 to 2019, fires >1 ha were mostly ignited by lightning (79%), compared to human causes (17%); 4% of the fires had an unknown cause. Lightning-caused fires were responsible for 96% of the total area burned and affected the whole territory ( Supplemental Material S7a (teco_a_2070342_sm6945.docx)), whereas human-caused fires were mainly located in the southeastern part of the territory and around the four Tłchǫ communities ( Supplemental Material S7b (teco_a_2070342_sm6945.docx)).
Fire-climate-environment associations
The fire season length was negatively correlated with the date of snowmelt estimated with the FWI daily model from 1965 to 2019 (Pearson's r = –0.77, p < 0.001). Thus, spring snow depth is a determinant of fire season length in the study area. The estimated mean fire season length was 159 days and varied from 127 days in 2004 to 180 days in 2010 ( Supplemental Material S8 (teco_a_2070342_sm6945.docx)). The mean starting and ending dates of the fire season were April 27 and October 2. The results showed significant trends from 1965 to 2019 for temperature (positive) and wind speed (negative) (Figure 3). Both LogAAB and AFO were positively correlated with seasonal temperature, fuel moisture content and fire-weather indices, while negatively correlated with seasonal precipitation (Table 2). It is worth noting that seasonal wind speed and relative humidity were neither correlated with LogAAB nor with AFO.
Table 2.
Pearson's r correlations between fire regime components (LogAAB and AFO) and potential top-down drivers (mean weather conditions, sum of precipitation, average of fuel moisture content and fire-weather indices) computed for the fire seasons from 1965 to 2019. See Table 1 for signification of acronyms.

Table 3.
Results from the modified version of the t-test on correlations between the spatial characteristics of the fire regime components (total fire occurrence and log-scaled area burned) and potential top-down climatic and bottom-up environmental (elevation) drivers. Fire regime components and explanatory variables were compiled for each 10-km grid cell and from 1965 to 2019. * p < 0.05.

The MARS and GAMM models predicted LogAAB with similar accuracy (adjusted R2 = 0.32 and 0.33, respectively, Supplemental Material S9a (teco_a_2070342_sm6945.docx)). The Pearson correlations between the observed and predicted values of LogAAB resulting from the MARS and GAMM models were 0.57 (p < 0.001) and 0.60 (p < 0.001), respectively ( Supplemental Material S9b (teco_a_2070342_sm6945.docx)). The variable selection procedure for the MARS model highlighted DMC as the best predictor of LogAAB, while the GAMM model selected both DMC and DC. MARS predicted AFO with higher accuracy than GAMM (adjusted R2 = 0.58 and 0.30, respectively, Supplemental Material S9c (teco_a_2070342_sm6945.docx)). The Pearson correlations between the observed and predicted AFO values resulting from the MARS and GAMM models were 0.76 (p < 0.001) and 0.57 (p < 0.001), respectively ( Supplemental Material S9d (teco_a_2070342_sm6945.docx)). Both models selected RH and ISI as the best predictors of AFO, but the MARS model also selected temperature, wind speed and precipitation.
Table 4.
Results from the modified version of the t-test on correlations between the spatial characteristics of the fire regime components (total fire occurrence and log-scaled area burned by 10-km grid cell for the period 2010–2019) and potential bottom-up environmental drivers (land cover types). Range and median correspond to the proportion of each 2010 land cover type by cell. ** p < 0.01 and * p < 0.05.

Spatial associations showed that temperature was positively correlated with fire occurrence on the study territory (Table 3). The presence of needleleaf forest was positively correlated with fire occurrence and area burned from 2010 to 2019, whereas the presence of grassland was negatively correlated with area burned (Table 4).
Identification and characterization of extreme wildfire years
We differentiated EWY from other years based on AAB. Nine years were in the top 20% (>80th percentile) of the AAB data set (1966, 1973 1976, 1979, 1998, 2008, 2013, 2014 and 2016) and were thus characterized as EWY (Figure 2a and Supplemental Material S5 (teco_a_2070342_sm6945.docx)). The nine EWY had AAB > 365 km2 (Figure 2a). The remaining 39 years having recorded fires had small to moderate AAB, each between 0.14 and 365 km2 depending on the year. Eight of the nine EWY (all but 1973) recorded fires >20,000 ha (Figure 2a), sometimes exceeding hundreds of thousands of hectares such as in 2014 ( Supplemental Material S10 (teco_a_2070342_sm6945.docx)). A total of 17 fires >20,000 ha were recorded from 1965 to 2019. These very large fires essentially occurred from the end of June to the middle of July (Figure 2b and Supplemental Material S11 (teco_a_2070342_sm6945.docx)). The northeastern parts of the study area have been less affected by very large fires (Figure 4). It is worth noting that EWY did not necessarily have high AFO; only three of the nine EWY were in the top 20% of the AFO distribution from 1965 to 2019 ( Supplemental Material S12 (teco_a_2070342_sm6945.docx)).
Thresholds associated with EWY occurrence
We observed that 75% of the EWY recorded an average temperature >14.7°C during the fire season compared to only 30% for non-EWY (Figure 5). We also noted significant differences in fuel moisture content and fire-weather conditions between EWY and non-EWY. The DC > 514 threshold was exceeded at least 36 days during the fire season for 75% of the EWY, compared to only 48% for non-EWY (Figure 5). The ISI > 7 threshold was exceeded at least 41 days during the fire season for 75% of the EWY compared to only 26% for non-EWY (Figure 5). The FWI > 19 threshold was exceeded at least 36 days during the fire season for 75% of the EWY compared to only 35% for non-EWY (Figure 5).
Environmental conditions (elevation and land cover types) of areas burned during EWY and non-EWY were not significantly different ( Supplemental Material S13 (teco_a_2070342_sm6945.docx)). Only one significant difference was observed for burned land cover occupied by water. We can explain this difference by the 30-m spatial resolution of the information used to describe land cover types and because fires during EWY were stopped on lakeshores (Figure 4).
Discussion
Characteristics of extreme wildfire years
Years with the largest annual areas burned (that we define as EWY) occurred essentially during the 1970s and 2010s on the territory of the Tłchǫ First Nation. EWY typically had lightning-induced fires >20,000 ha, and one fire even exceeded 300,000 ha in 2014. These very large fires occurred in the summer, during the annual peak of fire activity, more precisely from the end of June to the middle of July. This period of the year has previously been shown to have presented extreme fire risk from 1950 to 2001 in Alaska (Kasischke et al. 2002) and from 1980 to 1989 in the Canadian boreal forest (Stocks et al. 1998).
Figure 3.
Temporal changes and Mann-Kendall correlation tests (tau = Kendall rank correlation coefficient) from 1965 to 2019 on the Tłchǫ First Nation territory for (A) fire season length, (B-E) weather conditions, and (F-K) fuel moisture content and fire-weather indices. Blue lines represent significant trends built with a linear regression, while the shaded area represents the 95% confidence interval.

Top-down and bottom-up drivers of extreme wildfire years
While FWI components and seasonal precipitation did not vary significantly from 1965 to 2019, there were a decrease in seasonal average wind speed and an increase in temperature (Figure 3). The occurrence of several EWY during the last decade is related to high temperatures. Typical EWY fires occurred during periods with warm temperature (>14.7°C) and dry deep organic soil layers (DC > 514). These results confirm that top-down drivers have contributed to the occurrence of EWY during the last five decades in the boreal forest of central NT. Hence, as suggested in previous studies, our results emphasize the role of heat and drought in EWY occurrence (Kochtubajda et al. 2016) and more globally in areas burned in the Canadian boreal forest (Balshi et al. 2009; Parisien et al. 2011). Fires >20,000 ha were lightning-ignited, emphasizing the role of thunderstorm activity in the occurrence of EWY (Veraverbeke et al. 2017). Once ignited, EWY fires were mainly driven by fuel moisture content, particularly in deep organic layers (DC > 514), and by fire-weather conditions in terms of fire spread (ISI > 7) and fire risk (FWI > 19). Counterintuitively, the annual sum of precipitation did not differ significantly between EWY and non-EWY. Because annual precipitation is always relatively low in the study area, fire-conducive drought conditions are rather caused by above-average temperatures. Based on previous research on the sensitivity of fuel moisture to temperature and precipitation changes in the boreal forest (Girardin et al. 2009; Flannigan et al. 2016), we suggest that evapotranspiration exacerbated by warmer temperatures considerably reduces fuel moisture during EWY.
Figure 4.
Distribution of burned areas and year of ignition of fires > 20,000 ha from 1965 to 2019 on the territory of the Tłchǫ First Nation.

Increased temperature and drought conditions prevailing since the 1990s in the study area (Figure 3) contributed to a higher occurrence of large fires, especially during the last decade (Stephens et al. 2014). Our results corroborate other studies having demonstrated the influence of temperature and drought on burn rates, in Alaska for the period 1950–2003 (Duffy et al. 2005), and in the North American boreal forest for the period 1921–2008 (Flannigan et al. 2005; Girardin and Sauchyn 2008; Balshi et al. 2009; Senici et al. 2010; Parisien et al. 2011). Contrary to findings from previous studies having shown the role of wind speed in fire spread (Hirsch et al. 1998; Parisien et al. 2011), we observed that average wind speed during the fire season was not a significant driver of EWY occurrence, probably because this variable is highly dynamic at a (sub)daily scale and was smoothed on a season-wide basis in our analyses.
Figure 5.
Wilcoxon t-test between EWY and non-EWY from 1965 to 2019 for mean weather conditions, total precipitation, number of days with fire-weather indices exceeding the 75th percentile of the distribution of all fire season days, and length and start date of the fire season. *** p < 0.001; * p < 0.05. Red lines represent climatic or FWI thresholds associated with EWY occurrence. Note that fire-weather indices are unitless.

Top-down drivers had a more important influence on the occurrence of EWY than bottom-up drivers, for which, however, our analysis was based on a shorter time period (10 years) because of data availability constraints. Elevation and land cover type were not significant drivers differentiating burned areas during EWY and non-EWY, particularly because the elevation range within the study area is relatively small (323 m) and because the study area is dominated by needleleaf forests (Walker et al. 2018). Nevertheless, lakeshores were more affected by fires during EWY ( Supplemental Material S13 (teco_a_2070342_sm6945.docx)), suggesting the intensity and spread of very large fires, and the role of lakes as firebreaks, such as observed in the southwestern part of the study area in 2014 (Figure 4).
Because they were obtained from a relatively small area comprising the territory of the Tłchǫ First Nation, our results should not be extrapolated to the total extent of the NT. Nonetheless, our results are in line with previous studies having revealed that high-intensity wildfires in the NT largely affected organic soil layers of black spruce stands (Walker et al. 2018). In addition, we recognize that the available fire archives only cover the last five decades and that the early years of this time period were likely incomplete due to limitations in fire detection and mapping.
Climate change and risk of very large fires
Climate change has already had an effect on fire season length and annual area burned in the Canadian boreal forest (Gillett et al. 2004; Jain et al. 2017; Hanes et al. 2019). It can be anticipated that climate change will result in even larger fires and more frequent or severe EWY in the future. Indeed, the thresholds for top-down drivers of EWY could be exceeded more frequently in the next few decades if climate change trends continue in the same direction (Wang et al. 2014).
Projections from 24 climate models ( https://climatedata.ca/) indicate that the mean annual temperature in the study area could reach 2°C by 2100 compared to –6°C currently. Consequently, area burned is expected to increase in the future (Flannigan et al. 2005), and large fires could become more recurrent (Balshi et al. 2009). Warmer temperature could reduce the return interval of large fires and increase the risk of rapid reburning, possibly causing vegetation shifts from forest to woodland states, especially in black spruce stands (Barrett et al. 2011; Whitman et al. 2019). However, if the precipitation rate increases in the future, the fire risk could decrease, as observed in the 1980s in the study area (Flannigan et al. 2013; Kitzberger et al. 2017). Furthermore, warm conditions could increase the abundance of deciduous tree species across the study area (Boulanger et al. 2014; Chaste et al. 2019), which are less fire-prone than needleleaf species such as black spruce, which currently dominates the landscape. Conversely, increased fire activity could rejuvenate the forest, limiting fire-prone fuel (Parks et al. 2016; Baltzer et al. 2021). Hence, the effects of climate changes on future fire and vegetation remain uncertain in the northwestern boreal forest in Canada.
EWY occurrence thresholds for fire management on Indigenous territories
Higher recurrence of EWY and very large fires in the future in the Canadian boreal forest will directly affect Indigenous people's livelihood, security, health and traditional activities (Berkes and Davidson-Hunt 2006; Royer and Herrmann 2013; Dodd et al. 2018; Morarin 2020). By examining the fire regime over the past five decades on the Tłchǫ First Nation territory, we highlighted the extreme nature of the 2014 fire season in a short-term perspective and we have shown that other years have been comparable to 2014 even though they had lower annual areas burned. We also contributed to a better understanding of the drivers associated with EWY occurrence. Threshold values of temperature, fuel moisture content (DC) and fire behaviour (ISI, FWI) associated with EWY occurrence could be included in climate monitoring systems to more accurately predict EWY and to prevent large fire risk in a context of climate change. Thresholds can also inform management actions to limit the negative consequences of large wildfires on cultural keystone places and ecosystem services that benefit to the Tłchǫ people and neighboring First Nations.
Acknowledgments
The authors would like to thank Xiao Jing Guo (Laurentian Forestry Centre, QC, Canada) for helping with modelling, and Raphaël Chavardès (UQAT, QC, Canada) for providing a friendly review of the manuscript. We also thank the two reviewers who made comments that helped us improve the manuscript. Permission to work in the area was granted by the Tłchǫ government.
Funding
This work was supported by a grant to HA from Polar Knowledge Canada (# NST-1718-0014), a grant to DMG from the National Geographic Society (# EC-386R-18), and a Discovery Grant to MPG from the Natural Sciences and Engineering Research Council of Canada (NSERC).
Data availability statement
The data and R code used in this study are available at the following address: 10.5281/zenodo.6440362.