Effect of Topography on Persistent Fire Refugia of the Canadian Rocky Mountains

Persistent fire refugia, which are forest stands that have survived multiple fires, play an important ecological role in the resilience of mountainous forest ecosystems following disturbances. The loss of numerous refugia patches to large, high-severity fires in recent years is prompting the need to better understand drivers of fire refugia endurance. We investigate the role of topographic features on fire refugia survivorship based on pre-1950 fire regime conditions. Mapped refugia patches (n = 557) covering 28% of the forested landscape were used to develop three predictive models based on patch size (all sizes, <30 ha, <10 ha), as a function of explanatory variables describing several components of topography. Five topographic variables consistently favoured persistent fire refugia occurrence, though the ranking of explanatory variable importance varied among patch-size models. For the all-refugia model, elevation (23.7%), proportion of non-fuel at a 5000-m scale (20.3%), solar radiation (14.6%), Topographic Position Index at a 2000-m scale (10.1%), and distance from rivers (10.1%) were the top variables. The models’ predictive abilities were high, but decreased with larger patch sizes. We conclude that many suitable areas are currently unoccupied by fire refugia; that random elements affect their survivorship; and that additional environmental factors not considered in this study may contribute to their persistence. With changing climate and fire-regime conditions, careful fire and forest management considerations will be needed to limit future losses of persistent fire refugia forests.


Introduction
The simple definition of a forest fire refugium is a forested area, big or small, that has escaped burning.However, various studies define the state of what qualifies as refugia, or fire refugia, differently.Franklin et al. [1] define refugia as areas persisting through a disturbance that contain biological legacies, such as organisms, organic materials, and organically generated patterns, which can be incorporated into the recovering ecosystem.Generally speaking, fire refugia are the result of site characteristics and spatial locations favouring lower fire susceptibilities [2,3].Some consider residual island remnants within a recent fire to be refugia [4,5], whereas others are more concerned with persistences in which some patches of forests can survive multiple fires [2,3].In the context of this study, persistent fire refugia are examined.Persistent fire refugia patches are common in the Canadian Rocky Mountains of Alberta (hereafter 'Alberta Rockies'), which is remarkable considering how active fire has been on the landscape during the last millennium [6][7][8].In this fire-prone environment, fire refugia are forest stands that have escaped wildfires for centuries.They represent old-age forest patches embedded in or on the periphery of a historically active fire mosaic.Embedded patches tend to be the outcome of remnant patches having survived lower-intensity fires from short fire intervals and having built resistance over time through fuel development characteristics.In contrast, peripheral fire refugia patches do not tend to show evidence of surviving multiple fires, because fires do not repeatedly spread along the same footprint.In the context of this study, we defined and identified persistent fire refugia based on the guideline that stand ages had to be three times greater than the median fire return interval ("FRI") [9].In this case, fire refugia correspond to forests older than approximately 300 years given the historical fire regime of the Subalpine Natural Subregion of southern Alberta is regulated by a median FRI of 90 years [8].
The role of fire refugia in overall ecosystem functions and health is significant [10].Refugias' biological legacies contribute to faster ecosystem recovery following large and intense fire disturbances, as they form small ecological pools from which seeds can propagate [11].This ecological function is particularly important for mesic sites in high-elevation terrain where common pioneer species, such as lodgepole pine (Pinus contorta subsp.latifolia Loudon) and aspen (Populus tremuloides Michx.), are outside of their ecological range and cannot effectively contribute to forest recolonisation following fire [12].Fire refugia add to the overall biodiversity of a watershed by offering a diversified age-class mosaic [13,14] and by providing ecological niches specific to life organisms living in old-growth forests.In addition, forests older than 250 years are superior carbon sinks, due to increased carbon storage capacity through large wood volume, large coarse woody debris, and soil and ground organic matter [15].The well-developed moss and shrub layers of old forests also contribute to stream bank stability [16] and stream water quality through sediment retention, water filtering and cooling [17].The destruction of these stands by wildfire negatively impacts water quality for several years, due to the profuse release of sediments [18], nutrients [19,20], and heavy metals [21,22] that have accumulated over centuries.
Few studies have focused specifically on factors explaining the survivorship of fire residual patches, either ephemeral or persistent [2,4,5,33].The most recent studies examined numerous terrain metrics, including slope and aspect, and characteristics of the catchment basin, as well as topographic indices of wetness, convergence and wind shelter, along with climate and fire weather.They found that catchment slope, catchment aspect and topographic wetness index were the most significant model components for predicting fire residual patches density in northern New Mexico, USA [5].In contrast, flat and moderate topographic complexity within a rugged mountain environment produced a greater proportion of unburned area under benign and moderate fire weather condition scenarios [4].We expect that certain combinations of topographic elements can greatly increase the probabilities of survival and development of fire refugia zones in the Alberta Rockies, notably fire barriers.Patches of forests adjacent to features that do not support burning, or that reduce fire spread significantly, such as rocky outcrops, water bodies, or a vegetation type with low flammability, have better chances of escaping fire [9,32,34,35] (Figure A1).
In this study, we seek to identify topographic components favouring the persistence of fire refugia in the Alberta Rockies.Specific objectives consisted of: (1) assessing the contribution of topographic variables to fire refugia occurrence; (2) evaluating the functional form of the response between refugia occurrence and top explanatory variables; and (3) predicting and mapping the likelihood of refugia.Models were produced for different size thresholds to determine if controls of fire refugia occurrence vary Forests 2018, 9, 285 3 of 20 by refugia patch size.A large mountain landscape in the Alberta Rockies was chosen for its available fire refugia data set, which was produced from aerial photographs taken ca.1950, capturing historical burning patterns before climate change and significant improvement to the effectiveness of fire suppression.

Study Area
The study area forms a total land base of 911,951 ha, 47% of which is forested (425,011 ha), and represents a collection of adjacent land jurisdictions located in the Alberta Rockies.It includes: Banff National Park (BNP), the Whitegoat and Siffleur Wilderness Areas (WSWA) including the North Saskatchewan Valley from BNP east boundary to the middle of Abraham Lake, the Spray Lake Recreation Area (SLRA), and the Elbow watershed (Figure 1).The majority of the area is in the Subalpine Natural Subregion, with the exception of main drainage valley bottoms (i.e., Bow and North Saskatchewan valleys) between 1400 and 1600 m, which are part of the Montane Natural Subregion.Nearly 50% of the subalpine is composed of rock and ice, heavily dissecting the forest cover.The treeline is found at around 2200 m.Serotinous lodgepole pine is the main fire-adapted species colonizing the Montane, the lower Subalpine and it can also be found in the upper Subalpine on warm aspects [36].Engelmann spruce (Picea engelmannii Parry ex Engelm)-subalpine fir (Abies lasiocarpa (Hook.)Nutt.) forest communities favour mesic conditions and are found in wetter areas of the Montane, and on cool facing slopes and high elevations of the Subalpine [36].Alpine larch (Larix lyalii Parl.) are found along the treeline of the Subalpine, while pockets of aspen trees tend to be restricted to the Montane or the lowest elevations of the Subalpine.
Forests 2018, 9, x FOR PEER REVIEW 3 of 21 was chosen for its available fire refugia data set, which was produced from aerial photographs taken ca.1950, capturing historical burning patterns before climate change and significant improvement to the effectiveness of fire suppression.

Study Area
The study area forms a total land base of 911,951 ha, 47% of which is forested (425,011 ha), and represents a collection of adjacent land jurisdictions located in the Alberta Rockies.It includes: Banff National Park (BNP), the Whitegoat and Siffleur Wilderness Areas (WSWA) including the North Saskatchewan Valley from BNP east boundary to the middle of Abraham Lake, the Spray Lake Recreation Area (SLRA), and the Elbow watershed (Figure 1).The majority of the area is in the Subalpine Natural Subregion, with the exception of main drainage valley bottoms (i.e., Bow and North Saskatchewan valleys) between 1400 and 1600 m, which are part of the Montane Natural Subregion.Nearly 50% of the subalpine is composed of rock and ice, heavily dissecting the forest cover.The treeline is found at around 2200 m.Serotinous lodgepole pine is the main fire-adapted species colonizing the Montane, the lower Subalpine and it can also be found in the upper Subalpine on warm aspects [36].Engelmann spruce (Picea engelmannii Parry ex Engelm)-subalpine fir (Abies lasiocarpa (Hook.)Nutt.) forest communities favour mesic conditions and are found in wetter areas of the Montane, and on cool facing slopes and high elevations of the Subalpine [36].Alpine larch (Larix lyalii Parl.) are found along the treeline of the Subalpine, while pockets of aspen trees tend to be restricted to the Montane or the lowest elevations of the Subalpine.As a result of ecological succession, persistent fire refugia patches are largely composed of Engelmann spruce-subalpine fir [36].Lodgepole pine trees, which have a shorter life expectancy than the long fire-free intervals associated with fire refugia, tend to be sparse.However, they can still be an important component of refugia stands located on warm aspects and at lower elevations.Refugia represent decadent old-growth forests renewing themselves through small gap dynamics as a result of individual tree or tree cluster mortality [37][38][39].Stands are made up of scattered large-diameter trees, standing live and dead, with in-filling of multi-size trees forming a layered canopy [13,14,40].Large-diameter down woody debris in various degrees of decomposition and with various states of moss cover are scattered throughout [40].The stands often include open pockets that can be vegetated by a shrub layer and small trees, as well as dense pockets of in-filling subalpine fir trees [41].
The study area is located on the east side of the Continental Divide and is subject to a lightning strike shadow resulting in a considerably lower risk of lightning fires (Figure A2) [42].The Montane fire regime is almost exclusively driven by anthropogenic fires [8,43].In spite of being lightning-ignition limited, there is abundant evidence of historical fires having contributed to shaping the forest mosaic of subalpine landscapes for thousands of years [6].Fire frequency is significantly greater in the Montane, with a median FRI of about 30 years [8].In the Subalpine, the median FRI is 90 years [8], but important spatial variations in FRIs based on topographic locations have been recorded [9,44].

Fire Refugia Data Set
All jurisdictions had existing field-based fire history data, from which either complete-stand-origin maps or fire-refugia-only maps were produced [45][46][47][48][49]. Fire history sampling methods followed those developed for stand-replacing fires [50,51].While the sampling design at the time was not meant to study fire refugia, the large number of sample plots (1879) and tree cross-sections collected (6101) ensured all age cohorts present on the landscape were well sampled.Sampling plots were positioned on the edge of stand boundaries and four cross-sections on average were collected from dominant trees on each side of stand ecotones.The majority of fire refugia patches were identified from stand-origin dates, as few fire scars, if any, remained for dating old-growth forests in this stand-replacing fire regime environment.Not all fire refugia patches were visited due to difficulty of access on mountain flanks.Aerial photo interpretation using photography from ca. 1950, in combination with age data from nearby similar-looking stands, were used to assign stand-origin dates.Due to the distinctive photo tone and texture associated with old forests within a fire-generated forest mosaic, we are confident that fire refugia stands were correctly identified.
By using the rule of thumb of three times the length of the median FRI, forest stand-origin pre-dating 1700 were used to form the fire refugia data set.The mapping resolution ranged from 2 to 9 ha due to the variable map scales used.The pre-2000 stand-origin mapping was done at a coarser scale of 100 m, whereas fire refugia patches for the Elbow watershed were digitised on-screen from 2007, 0.5 m resolution, ortho-rectified imagery using 1949 aerial photography and fire history plots to guide fire refugia mapping.In total, 557 fire refugia polygons were identified within the study area.Collectively, they amount to 117,515 ha, representing 27.7% of the forested land base (425,011 ha).Fifty percent of patches were less than 30 ha in size, while 37% were less than 10 ha.In spite of the important number of small-sized refugia patches, they amount to only 1.8% and 0.7% of the total fire refugia area, respectively.Due to the large number of small-sized refugia patches available, we used these two patch size thresholds to verify whether small size refugia responded differently to topographic features than those of greater size.The spatial distribution of fire refugia is strongly associated with the Subalpine Natural Subregion, as few patches were recorded at valley bottom in the Montane.A visual assessment of the data showed a fairly even distribution of refugia patches throughout most of the subalpine watersheds within the study area.
Forest species composition is often a controlling element of fire regime characteristics [52,53], and fire refugia could potentially respond differently to topographic features based on vegetation.In our study, we decided not to address this component, as our refugia patches showed an overwhelming dominance of spruce fir stands, which would make a topographic species response analysis difficult to interpret.We were also mindful that in rugged mountains, both aspect and elevation govern the micro-climate and ensuing species assemblage.In that regard, high correlations of FRI response between spruce and high elevation and cool aspects have been identified in the Alberta Rockies [9].

Topographic Variables Data Set
We considered topographic variables that could ecologically explain consistent changes in fire behaviour and fire spread patterns to the benefit of long-term forest survival during wildfire events.Following a literature review of studies having examined topographic elements in relation to fire residual patches [2,4,5,32], we shortlisted ten topographic variables believed to strongly influence persistent fire refugia: elevation, slope, solar radiation as a surrogate for aspect, topographic position index (TPI) [54], maximum rate of slope curvature change (maxcurv), which detects convergent and divergent terrain, proportion of non-fuel (NF), and distance from rivers (Strahler stream order ≥5) as a mean to capture headwaters of side drainages.Total annual solar radiation (Mwh/m 2 ) was calculated using the Area Solar Radiation tool from the Spatial Analyst Extension [55] in ArcGIS.The proportion of non-fuel (rocks and water) was calculated from the Ecological Land Classification Map for Banff National Park [56], from the Alberta Ground Cover Classification Mosaic [57,58] for the Elbow Watershed, and from 1:50,000 National Topographic Survey maps for the rest of the landscape.The remainder of the topographic variables were derived from a 30-m digital elevation model (DEM) produced from Canadian Digital Elevation Data at a scale of 1:250,000 [59].
Many area-based environmental variables are scale dependent [60].Consequently, we assessed various scales as circular windows for TPI, NF, and maxcurv.TPI was tested at a scale of 300 m, 500 m, 1000 m, 2000 m, and 4000 m; NF was tested at a scale of 300 m, 500 m, 1000 m, 2000 m, and 5000 m; whereas maxcurv was tested at a scale of 300 m, 500 m, 1000 m, and 2000 m.This brought the total number of variables examined to 17.All topographic variables were mapped at 100-m resolution using Idrisi TerreSet [55,61,62].
The following variables are commonly assessed in topographic analysis studies, but we chose to exclude them.Topographic ruggedness index (TRI) had a strong association with rocky terrain above treeline, whereas less rugged values largely corresponded to the forested area, making it a poor indicator of refugia on the landscape.Similarly, topographic wetness index (TWI) had low scores across all angled terrain and only showed high scores for valley bottoms of main rivers.We also opted against any catchment metric variables aiming to capture cool-air-pool (CAP) drainages [63].The Canadian Rockies do not have canyon-type ecosystems, and the few narrow canyons existing do not contain much treed vegetation that could become climate refugia.Hanging catchments in the upper Subalpine do not tend to be cold-air sinks during fire spread events associated with hot days.During the day, heated air rises along slopes, an effect that is exacerbated on warm aspects.Thus, the likelihood of a catchment metric being associated with persistent fire refugia is unlikely.

Topographic Variable Selection
The exploratory variable selection was performed through a series of steps to minimise the number of variables used in the prediction model.The Spearman rank correlation test served to eliminate highly correlated variables with a |r| ≥ 0.7.In order to select the best variable from each correlated group of variables, we built single-variable generalised additive models (GAMs) predicting refugia presence for each independent variable.The variable with the higher deviance explained in each correlated group was retained.This resulted in the retention of nine independent variables to be considered for analysis (Table 1).

Model Building and Analysis
We constructed boosted regression tree models (BRTs) of refugia presence-absence for three refugia size classes: <10 ha, <30 ha, and all refugia.BRT is a machine-learning modelling approach that fits multiple regression trees to attain high predictive performance that can account for nonlinear relationships [64].BRTs do not require any a priori model specification or test of hypothesis [65], and can account for complex interactions among variables.The technique builds regression trees through "boosting", using a stage-wise model-fitting process to selectively build each subsequent regression tree to minimise the residuals of the previous tree.At each individual tree, the training data is split by a ratio (i.e., the bag fraction), which defines the proportion of the training data to use to build a tree.
To ensure that each refugia size class model was generated in an unbiased fashion, all three models used the same nine variables.BRT model-building methods followed [66].We used a stratified random sampling approach to extract 10,000 sample points for each of the three refugia size classes, and 10,000 absence (i.e., non-refugia) points in other forested areas.Each refugia size model was built from a random subset of 400 points from the presence pool of points and 400 points from the absence pool of points.We repeated this process 100 times and averaged the resulting 100 BRT models to create an ensemble model for each size class.This approach allowed for selecting different sample points between each individual sub-model in order to capture variation in fine-scale topographic variables and to limit the stochasticity in model outcomes.
All models were fitted using R gbm package [67].We assumed a Bernoulli squared error loss function, given that our response variable was binomial.We built the models using 3-fold cross validation, a bag fraction of 0.5, and a tree complexity of 5. Tree complexity refers to the number of nodes in a tree, and we deemed a 5-node tree to be sufficiently deep to capture moderately complex interactions between topographic variables.We used a learning rate of 0.02, 0.015, and 0.005 for the <10 ha, <30 ha, and all-refugia models, respectively.We adjusted this learning rate for the three models such that each fire refugia class model built between 500 and 1500 regression trees for each submodel of the ensemble.The number of trees in each submodel corresponded to the point at which adding more trees would reduce cross-validated model performance.We performed model validation using k-fold cross-validation, which randomly divides the sample into k subsamples, of which a single subsample is retained for model validation, and the remaining folds are used for model training.Model performance was assessed using cross-validated area under the receiver operating characteristic curve (AUC), a measure of how well the model can distinguish between refugia and non-refugia.All test subsets are considered and the results from all k folds are then averaged.Variable importance was assessed for each variable of the three refugia-size models.As computed from the gbm package [67], variable importance consists of the number of times the variable was selected as a tree node and the improvement resulting from including those nodes.We plotted the relative contribution for each variable as a percentage of the total deviance explained.
As a complement to the variable importance, we evaluated the relationships of likelihood of fire refugia to the top explanatory variables for each refugia-size model.The top five variables of the all-refugia model were retained for this analysis.These relationships were visualised by plotting the partial dependence to each explanatory variable of interest.Partial dependence plots represent the marginal effect of a given variable by holding the values of all other variables of the model constant at their mean value, thereby statistically controlling for the effect of these variables.The variable importance and partial dependence plots can be interpreted together to provide both the relative importance and the relationship of each explanatory variable to fire-refugia presence.
Finally, model predictions were mapped for each refugia-size model using the raster package in R [68].The modelled likelihood of fire refugia was also mapped for each observed refugia patch, in order to assess model predictions on a per-patch basis.

Results
Mean model cross-validated AUC were 0.71, 0.77, and 0.84 for the all-refugia, <30 ha, and <10 ha models, respectively.BRTs consistently identified the same four topographic variables, in different orders of relative influence, as significant drivers of small, medium and all sizes of persistent fire refugia patches: elevation, the proportion of non-fuel in a 5000-m moving window, solar radiation (i.e., aspect) and distance from rivers (Figure 2).They represent a combined relative influence of 69%, 58% and 64% for the all-refugia, <30 ha, and <10 ha models, respectively.The other lead topographic variables, ranking fifth and sixth, according to the order of their relative influence among the three models, included TPI in a 2000-m moving window, the maximum rate of slope curvature change in a 1000-m moving window and proportion of non-fuel in a 1000-m moving window.While distance from rivers had less than 10% influence in the all-refugia model, its influence was 20% and 15% in the <10 ha and <30 ha models, respectively.Solar radiation had a consistent influence of 12 to 15% across all models, whereas some of the highest divergence in influence was proportion of non-fuel in a 1000-m window.At 12% influence, it ranked fifth for the <30 ha model, while it received the lowest influential score in the other models.
Forests 2018, 9, x FOR PEER REVIEW 7 of 21 variable was selected as a tree node and the improvement resulting from including those nodes.We plotted the relative contribution for each variable as a percentage of the total deviance explained.As a complement to the variable importance, we evaluated the relationships of likelihood of fire refugia to the top explanatory variables for each refugia-size model.The top five variables of the allrefugia model were retained for this analysis.These relationships were visualised by plotting the partial dependence to each explanatory variable of interest.Partial dependence plots represent the marginal effect of a given variable by holding the values of all other variables of the model constant at their mean value, thereby statistically controlling for the effect of these variables.The variable importance and partial dependence plots can be interpreted together to provide both the relative importance and the relationship of each explanatory variable to fire-refugia presence.
Finally, model predictions were mapped for each refugia-size model using the raster package in R [68].The modelled likelihood of fire refugia was also mapped for each observed refugia patch, in order to assess model predictions on a per-patch basis.

Results
Mean model cross-validated AUC were 0.71, 0.77, and 0.84 for the all-refugia, <30 ha, and <10 ha models, respectively.BRTs consistently identified the same four topographic variables, in different orders of relative influence, as significant drivers of small, medium and all sizes of persistent fire refugia patches: elevation, the proportion of non-fuel in a 5000-m moving window, solar radiation (i.e., aspect) and distance from rivers (Figure 2).They represent a combined relative influence of 69%, 58% and 64% for the all-refugia, <30 ha, and <10 ha models, respectively.The other lead topographic variables, ranking fifth and sixth, according to the order of their relative influence among the three models, included TPI in a 2000-m moving window, the maximum rate of slope curvature change in a 1000-m moving window and proportion of non-fuel in a 1000-m moving window.While distance from rivers had less than 10% influence in the all-refugia model, its influence was 20% and 15% in the <10 ha and <30 ha models, respectively.Solar radiation had a consistent influence of 12 to 15% across all models, whereas some of the highest divergence in influence was proportion of non-fuel in a 1000-m window.At 12% influence, it ranked fifth for the <30 ha model, while it received the lowest influential score in the other models.The relationships, expressed as partial dependence plots, of the likelihood of fire refugia to individual explanatory variables differed substantially among variable and among refugia-size models (Figure 3).For the all-refugia model, which is largely weighted towards patches larger than 1 km 2 , trends tend to be linear, whereas patterns depicted in small-sized refugia models fluctuate more as patch size decreases.Generally, there are more refugia as elevation increases, and with an increasing proportion of rocks in a 5000-m window.As solar radiation increases from cool to warm constant in the all-refugia model whereas the tendency is for increased number of small refugia away from main river valleys, towards headwaters.Refugia response to TPI in a 2000-m window is substantially variable, based on patch size.The all-refugia model suggests there are more refugia on towering landforms, whereas the <30 ha model indicates that refugia can be higher than the surrounding areas, as well as lower in depressed landforms.The <10 ha model points to a prevalence of refugia associated with uniform terrain and slightly lower than the surrounding land.The relationships, expressed as partial dependence plots, of the likelihood of fire refugia to individual explanatory variables differed substantially among variable and among refugia-size models (Figure 3).For the all-refugia model, which is largely weighted towards patches larger than 1 km 2 , trends tend to be linear, whereas patterns depicted in small-sized refugia models fluctuate more as Forests 2018, 9, 285 9 of 20 patch size decreases.Generally, there are more refugia as elevation increases, and with an increasing proportion of rocks in a 5000-m window.As solar radiation increases from cool to warm aspects, the presence of refugia diminishes.Fire refugia response to distance from rivers is relatively constant in the all-refugia model whereas the tendency is for increased number of small refugia away from main river valleys, towards headwaters.Refugia response to TPI in a 2000-m window is substantially variable, based on patch size.The all-refugia model suggests there are more refugia on towering landforms, whereas the <30 ha model indicates that refugia can be higher than the surrounding areas, as well as lower in depressed landforms.The <10 ha model points to a prevalence of refugia associated with uniform terrain and slightly lower than the surrounding land.The spatial prediction maps of fire refugia probability show an area reduction of high probabilities (>0.7) with decreasing patch size (Figure 4).A calculation of the proportion of forested landscape falling in refugia probabilities >70% amounted to about 10% across all models, whereas the proportion of landscape showing 50% likelihood of refugia occurrence varied from 46, to 33 and 26% for the all-refugia, <30 ha and <10 ha models, respectively.When we assessed how well the models performed against the actual position of fire refugia patches in the study area, refugia falling in the >0.7 probability zone amounted to 23, 38 and 57%, for the all-refugia, <30 ha and <10 ha models, respectively.However, the same assessment using probability zones >0.5, which is a positive likelihood of refugia occurrence, showed the models performed much better, with 71, 74 and 84% of refugia observations being well classified (Table A1, Figure A3).The spatial prediction maps of fire refugia probability show an area reduction of high probabilities (>0.7) with decreasing patch size (Figure 4).A calculation of the proportion of forested landscape falling in refugia probabilities >70% amounted to about 10% across all models, whereas the proportion of landscape showing 50% likelihood of refugia occurrence varied from 46, to 33 and 26% for the all-refugia, <30 ha and <10 ha models, respectively.When we assessed how well the models performed against the actual position of fire refugia patches in the study area, refugia falling in the >0.7 probability Forests 2018, 9, 285 10 of 20 zone amounted to 23, 38 and 57%, for the all-refugia, <30 ha and <10 ha models, respectively.However, the same assessment using probability zones >0.5, which is a positive likelihood of refugia occurrence, showed the models performed much better, with 71, 74 and 84% of refugia observations being well classified (Table A1, Figure A3).

Discussion
Given the unique perspective stemming from a rich data set covering a large-sized territory, our study provides a valuable contribution to understanding persistent fire refugia in the Canadian Rockies.The study design made it possible to define and validate spatial predictions, as well as having the flexibility to assess the role of topography on different sizes of refugia patches.The robust and meaningful model results not only identified the topographic variables that had the strongest response on persistent fire refugia, but also captured the response trends when lead variables were combined together.The validation of all three spatial prediction models of refugia size probabilities on the landscape allowed us to better understand and capture where fire refugia patches are most likely to persist.When we compared observed refugia patches against the probability of refugia occurrence map, understanding that the prediction map was built from those same refugia patches, we readily discerned which patches may truly be persistent as a result of their topographic location, versus those which may have remained more by chance (green and blue zones), given that they do not have a particular topographic setting that limits fire ignition and spread.
Elevation, solar radiation (i.e., aspect), proportion of non-fuel in a 5000-m window, and distance from rivers were consistently the most influential topographic variables across all models, but responses varied based on fire refugia size.The all-refugia model, which is heavily weighted towards large-sized patches, showed gentle linear response trends for all lead variables.This is explained by the large extent of some refugia patches, which can fill entire side drainages, and thus individual patches can cover a wide range of topographic variations.In contrast, patches <30 ha showed more sensitivity to topographic elements-notably, elevation and distance to rivers-but the best model response was obtained from refugia patches <10 ha.The small-sized model points to a strong positive linear relationship with elevation, but the other variables tested showed fluctuating (i.e., non-linear) refugia responses, which are common with environmental variables [69].The environmental conditions in which fires burn can be highly variable due to important changes in the weather, vegetation types, and terrain, or from subtle changes in the relative humidity, temperature, wind direction and fuel moisture content, which ultimately dictate the rate of fire spread and the direction of movement of the fire.In this context, it is important to acknowledge that an element of randomness exists in the survivorship of forest stands to fire, even for persisting fire refugia [4].
As expected, elevation was a significant driver of persistent refugia in the Alberta Rockies, given that the longest FRIs are found at higher elevations [9,70,71].This is explained by the cooler air, higher relative humidity conditions, greater precipitation amounts and prolonged snow cover that are associated with high elevations (>1800 m).It has been determined that for every metre of increased elevation, there is a 0.33% chance of decreased probability of burning [9]; thus, stand ages can easily exceed three hundred years nearing treeline and in proximity to headwaters.
Cooler facing slopes were found to be more likely to host fire refugia, but this was not a steadfast rule, notably for the smallest refugia patches.Receiving less solar radiation, especially during the early spring and late autumn, means that fuels from north-and east-facing slopes retain more moisture and require a greater number of warm and rain-free days to be available for burning [72].Similar to high elevations, cool aspects have a longer lasting snow cover and hold a higher relative humidity than warm-facing slopes [73].In the Alberta Rockies, high elevations combined with cool aspects produce the longest FRIs and allow for the development of old-growth forests [9].Thus, it is logical that this synergy would also favour the persistence of fire refugia.
Similar to studies in other fire-prone areas, forest patch isolation as a result of rocky ridges, lakes or wetlands is a strong determinant of old forests [32,35,74].The proportion of land without vegetation (non-fuel) in a 5000-m moving window captured the high proportion of rocks associated with small valleys bound by rocky ridges on three sides.This represents the most rugged regions of the mountains, where rocky ridges above treeline heavily dissect the forest cover and inhibit fire spread from one drainage to the other.After elevation, proportion of nonfuel was the most influential driver for the persistence of refugia patches in the all-refugia model, as many of the largest fire refugia patches (>100 ha) are located in narrow valleys surrounded by rocky ridges.
Distance from rivers was an important driver in all refugia-size models, but more so for the <30 ha and <10 ha patches.While the response is not linear, there is generally an increase in fire refugia occurrence as a function of the distance from rivers.In our study area, a first wave of refugia tends to occur in river valleys at distances less than 5 km and includes small patches along the treeline, as well as larger patches along the river banks of north-facing slopes.The second upsurge in refugia probability represents patches located in secondary and tertiary drainages (Strahler stream orders 1, 2, 3), which are furthest from the main valleys and where fire is unlikely to spread deeply into these small side drainages.The influence of distance from rivers largely resides in the probabilities of ignition, which are the greatest at valley bottom of important drainages where human travel corridors are found (roads, trails).
Historically, anthropogenic fires had a tendency to originate from main valleys in the Montane [43] and propagate into smaller drainages, provided that burning conditions were favourable and winds shifted in such a manner that flames could be pushed into small drainages [8].In combination with favourable topographic features, areas with low probability of ignitions, from both anthropogenic and lightning sources, are the likely explanations for the persistence of the largest patches of fire refugia.We presume that smaller-sized patches are the result of a greater exposure to fire over centuries, where larger refugia stands were dissected by one or several fires, or where fire residuals within fire perimeters built resistance to burning over time.
When we compare our results to those of similar studies, we corroborate elevation [2], aspect [4,5,75,76], fuel breaks [32], relative topographic position [4,5], convergence [4,5] and headwaters [2] as also being important topographic drivers of fire refugia occurrence.However, due to the different landscapes, vegetation and fire regimes, slope did not make the top nine variables in our study, although it ranked highly in others [4,5].This outcome may also be due to the use of different variables, as interactions between some covariates can mask the true effect of weaker variables, as well as using different spatial scales of analysis [60].Another reason for differences between studies could be from the type of fire refugia assessed.Our refugia were largely part of the overall forest fabric where patches tend to be in periphery of fires rather than residuals within recent burns.These discrepancies in results among studies may suggest that topographic drivers are somewhat different between ephemeral and persistent refugia (Figure 5).
The assessment of land cover from the spatial predicted maps of fire refugia revealed interesting facts.Only about 10% of the forested landscape offers ideal topographic conditions in which the probabilities of refugia occurrence are estimated to be greater than 70%.This value was similar across all three refugia-size models.It highlights how vulnerable fire refugia may become in a warming climate if only 10% of the forested Subalpine offers the best sheltering topographic conditions from burning.However, given that the current cover of refugia patches already amount to 28% of the forested area, and that the majority of refugia patches were classified in zones with predicted probabilities of fire refugia greater than 50%, we deduce there could be up to 46% of the forested landscape capable of hosting large-sized persistent refugia.This number diminishes to 25% of the forested landscape for small-sized refugia patches.When we compare the current refugia cover to that from the predicted probabilities, it is clear that forests growing on topographic features amenable to the formation of fire refugia are not necessarily protected from burning.Keeping in mind that our data set reflects pre-1950 fire regime conditions, when there were more fires burning and before the current climate warming, this outcome puts forward additional questions with regard to other driving factors that may influence the sustainability of persistent fire refugia.
our study, although it ranked highly in others [4,5].This outcome may also be due to the use of different variables, as interactions between some covariates can mask the true effect of weaker variables, as well as using different spatial scales of analysis [60].Another reason for differences between studies could be from the type of fire refugia assessed.Our refugia were largely part of the overall forest fabric where patches tend to be in periphery of fires rather than residuals within recent burns.These discrepancies in results among studies may suggest that topographic drivers are somewhat different between ephemeral and persistent refugia (Figure 5).Given that the extent of the area with favourable topographic features is 39% greater than the current fire refugia cover, we conclude that other factors must be at play.First, there is an element of randomness to the distribution of fire refugia and their length of persistence (Figure A4), and additionally, topography itself may play a smaller role in the overall fire environment than anticipated, notably during large fire events [77]; second, other factors in combination with topography likely assist the ability of refugia to persist for centuries.Such factors could include the age of a fire refugium and its associated stage of fuel development.In subalpine old-growth forests, moist fuel conditions normally arise from a deep moss layer, large-diameter tree boles and a thick understory cover [40].Such fuel characteristics have been found to be an important contributing factor of fire resistance of boreal forest fire refugia [78].The amount and large diameter of snags and coarse woody debris associated with decadent forests require a drying time lag of 52 days following rain or snowmelt [72], making these stands more difficult to ignite unless there have been prolonged drought conditions.
The presence of recent burns and stands aged less than 50 years have been shown to limit the size of subsequent fires and influence future fire spread patterns [35,[79][80][81].The historical fire regime, which was largely anthropogenic on the east slopes of the Canadian Rocky Mountains, documented shorter pre-1950 FRIs at the valley bottom of primary drainages, which produced mixed-severity fire effects (i.e., fire intensities were variable, and fires were not entirely lethal) [8,43,72].The stand age-resistance factor to fire ignition and spread was thus an additional contributing element enabling the development and persistence of fire refugia, notably when mean forest ages were younger within a watershed, or immediately adjacent to existing fire refugia.Coupled with key topographic conditions conducive to maintaining moister conditions such as high elevation and cool aspects [9], or being protected by fuel breaks [32], some forest stands were able to escape fire over decades and centuries during the historical fire regime.These older stands essentially built resistance against moderate-intensity fires, and even against the intensities produced from young lodgepole pine stands crowing; a typical fire behaviour associated with stands less than 60 to 70 years old, due to their high tree densities and low branches [82,83].
In this contemporary era of very large fires following decades of fire exclusion policies and a warming climate [84], we are questioning whether combinations of factors that promoted fire refugia in the past may no longer effectively do so in the future.In recent years, we have observed burning conditions in the Canadian Rockies that are conducive to producing large, high-severity fires with few island remnants, where fires burn at higher elevations, and frequently up to the treeline.Similar observations have been made in the Sierra Nevada [85] and the Northern Rockies [86].Fuel development in tandem with sustained droughts and erratic wind events during burning are now making it exceptionally difficult for established fire refugia to survive [33,35,87].These severe burning conditions are believed to be the results of combined factors directly affecting the fire environment once a fire has ignited.First, the warming climate trend is extending the fire season globally, both at the start and end [88]; second, long-lasting blocking ridges of high pressure have been leading to sustained drought conditions allowing for large-diameter standing fuels to dry out and burn more easily [89]; and third, the combination of fire exclusion and effective fire suppression over the last several decades has made the landscape forest mosaic more uniform from maturing forest conditions (e.g., densification, development of laddering fuels), which lead to higher intensity fires [85,90].
Following a fire hiatus of 80 years in many subalpine valleys of the Canadian Rockies, a homogeneous cover of dense mature lodgepole pine forests has developed [90] and offers little difference in canopy height or stand structure with nearby old forests or fire refugia stands.Under such fuel conditions, we believe fire suppression will have little effect, notably during extreme fire weather events.As a result, we have observed the decimation of very large old-growth patches from recent wildfires in the mountain national parks (e.g., Waterton Lake Kenow 2017 fire, Banff Spreading Creek 2015 fire).A troublesome post-fire outcome of these extreme-severity fires is the consumption of all soil organic materials, including coniferous seed banks (pers.obs.Spreading Creek Fire, Banff National Park).Under such circumstances, ecosystem resilience may be compromised considering ephemeral or persistent fire refugia are almost absent and vast expanses of area burned are far from seed sources, such as a refugium or the unburned forest mosaic.

Conclusions
Fire refugia are increasingly facing adverse external pressure that could eradicate these biological legacies.There is an even greater cause for concern, as results show topographic components, while having significant influence on the persistence of fire refugia, may not be as effective or as relevant as first thought under a warming climate and current fuel conditions.We are likely seeing the cascading effect of unintended consequences from decades of fire exclusion at lower elevations, exacerbated by severe droughts, and which are predicted to become more frequent in this climate warming trend.It will be necessary to further document the location and relative age of recently burned fire refugia patches to determine if randomness or a changing fire regime was the likely cause of destruction.If lost refugia are estimated to be more than 400 to 500 years old, there is a greater need for concern due to sensitive and very unique habitats being lost.With evolving climate and fuel conditions, contemporary burning regimes are changing, and the sustainability of persistent fire refugia over the next century will be challenged.The outcome of this research, derived from historical fire regime conditions, offers a benchmark from which to compare future trends.We advocate a careful management of landscapes containing persistent fire refugia (i.e., old-growth forests) moving forward.

Figure 1 .
Figure 1.(A) Study area located in the Canadian Rocky Mountains of Alberta, Canada, composed of four jurisdictions; (B) Distribution of fire refugia patches (green), within the forested zone (beige).Refugia smaller than 30 ha are highlighted using pink dots.

Figure 1 .
Figure 1.(A) Study area located in the Canadian Rocky Mountains of Alberta, Canada, composed of four jurisdictions; (B) Distribution of fire refugia patches (green), within the forested zone (beige).Refugia smaller than 30 ha are highlighted using pink dots.

Figure 2 .
Figure 2. Relative importance of explanatory variables for each refugia size class, as percent of total deviance explained.

Figure 2 .
Figure 2. Relative importance of explanatory variables for each refugia size class, as percent of total deviance explained.

Figure 3 .
Figure 3. Partial dependence plots of fire refugia probability (y-axis) on five selected exploratory variables shown for the three fire refugia size-class models.The variables represent the five most important, as determined by the all-refugia model.Error bands represent 95% confidence intervals.

Forests 2018, 9 , 21 Figure 3 .
Figure 3. Partial dependence plots of fire refugia probability (y-axis) on five selected exploratory variables shown for the three fire refugia size-class models.The variables represent the five most important, as determined by the all-refugia model.Error bands represent 95% confidence intervals.

Figure 5 .
Figure 5.This photograph shows ephemeral fire refugia patches at lower elevations (in the foreground outlined in yellow) that survived the most recent fire.Ephemeral patches had burned during the previous two fires and are maintaining lodgepole pine stand communities, whereas the persistent fire refugia patches (pink outline) tucked at the headwall of a side drainage have avoided multiple fires and largely consist of spruce trees.

Figure 5 .
Figure 5.This photograph shows ephemeral fire refugia patches at lower elevations (in the foreground outlined in yellow) that survived the most recent fire.Ephemeral patches had burned during the previous two fires and are maintaining lodgepole pine stand communities, whereas the persistent fire refugia patches (pink outline) tucked at the headwall of a side drainage have avoided multiple fires and largely consist of spruce trees.

Figure A1 .
Figure A1.Examples of subalpine environment in which fire refugia patches are found: headwaters, near treeline, against fuel breaks.(A) 1949 aerial photo showing persistent fire refugia patches outlined in red following a 1910 fire that burned most of the valley; (B) Fuel breaks such as rock slides and treeline inhibit fire spread and enhances survivorship to fire; (C) Small side valley surrounded with rocks is completely composed of fire refugia forest; (D) Persistent fire refugia at high elevations avoided a prescribed burn.

Figure A1 .
Figure A1.Examples of subalpine environment in which fire refugia patches are found: headwaters, near treeline, against fuel breaks.(A) 1949 aerial photo showing persistent fire refugia patches outlined in red following a 1910 fire that burned most of the valley; (B) Fuel breaks such as rock slides and treeline inhibit fire spread and enhances survivorship to fire; (C) Small side valley surrounded with rocks is completely composed of fire refugia forest; (D) Persistent fire refugia at high elevations avoided a prescribed burn.

Figure A2 .
Figure A2.Average annual lightning strike density recorded over 5 km × 5 km grid cells.The lightning strike shadow corresponds to the zone with less than approximately 10 strikes/year.Fire refugia patches are shown as white outlines.Lightning data source: Government of Alberta, 2006 to 2015.

Figure A2 .
Figure A2.Average annual lightning strike density recorded over 5 km × 5 km grid cells.The lightning strike shadow corresponds to the zone with less than approximately 10 strikes/year.Fire refugia patches are shown as white outlines.Lightning data source: Government of Alberta, 2006 to 2015.

Figure A3 .
Figure A3.Model spatial prediction weighted means summarised by individual refugia for the allrefugia model (A); refugia <30 ha model (B); and refugia <10 ha model (C); refugia smaller than 30 ha are represented as points.

Figure A4 .
Figure A4.Persistent fire refugia patches (pink outline) are near treeline and at the headwaters of the side drainage.Whereas a refugium located mid-slope (blue outline) that was surrounded by forest may also be persistent, but owes its survival to random effects or factors other than topography.

Figure A3 .
Figure A3.Model spatial prediction weighted means summarised by individual refugia for the all-refugia model (A); refugia <30 ha model (B); and refugia <10 ha model (C); refugia smaller than 30 ha are represented as points.

Figure A4 .
Figure A4.Persistent fire refugia patches (pink outline) are near treeline and at the headwaters of the side drainage.Whereas a refugium located mid-slope (blue outline) that was surrounded by forest may also be persistent, but owes its survival to random effects or factors other than topography.

Figure A4 .
Figure A4.Persistent fire refugia patches (pink outline) are near treeline and at the headwaters of the side drainage.Whereas a refugium located mid-slope (blue outline) that was surrounded by forest may also be persistent, but owes its survival to random effects or factors other than topography.

Table 1 .
Description of lead topographic variables considered.Value minimums, maximums, means and standard deviations are presented.