Burning Potential of Fire Refuges in the Boreal Mixedwood Forest

: In boreal ecosystems, wildﬁre severity (i.e., the extent of ﬁre-related tree mortality) is affected by environmental conditions and ﬁre intensity. A burned area usually includes tree patches that partially or entirely escaped ﬁre. There are two types of post-ﬁre residual patches: (1) patches that only escaped the last ﬁre; and (2) patches with lower ﬁre susceptibility, also called ﬁre refuges, that escaped several consecutive ﬁres, likely due to particular site characteristics. The main objective of this study was to test if particular environmental conditions and stand characteristics could explain the presence of ﬁre refuges in the mixedwood boreal forest. The FlamMap3 ﬁre behavior model running at the landscape scale was used on the present-day Lake Duparquet forest mosaic and on four other experimental scenarios. FlamMap3 was ﬁrst calibrated using BehavePlus and realistic rates of ﬁre spread obtained from the Canadian Fire Behavior Prediction system. The results, based on thousands of runs, exclude the effects of ﬁrebreaks, topography, fuel type, and microtopography to explain the presence of ﬁre refuges, but rather highlight the important role of moisture conditions in the fuel beds. Moist conditions are likely attributed to former small depressions having been ﬁlled with organic matter rather than present-day variations in ground surface topography.


Introduction
Fire is one of the dominant ecological drivers affecting vegetation patterns and dynamics in the circumboreal region [1,2]. Fire effects vary spatially depending on fire behavior [3,4] and the fire return interval. A burned area usually includes residual patches that partially or entirely escaped fire [5][6][7][8]. Two types of post-fire residual patches have been distinguished in eastern North American boreal mixedwood forests [9]: (1) "transient residual patches" that only escaped the last fire, probably due to peculiar but temporary unsuitable conditions for fire propagation; and (2) "fire refuges" that escaped several consecutive fires, likely due to specific site conditions. Although fire refuges represent a small proportion of the total area burned, they could provide unique habitats in post-fire successional landscapes [8]. Indeed, the ecological continuity recorded in fire refuges (unlike in transient residual patches, which only escaped the last fire) [9] could provide refuges for species with specific biodiversity

Study Area
The reference forest mosaic used in this study is an 11,000-ha natural forest mosaic encompassing the Lake Duparquet Research and Teaching Forest (Figure 1). The study area is located in the eastern Canadian boreal mixedwood forest [26], and was previously used to test the effect of landscape composition on fire size distribution [27]. The studied landscape is characterized by balsam fir (Abies balsamea (L.) Mill.), paper birch (Betula papyrifera Marsh.), white spruce (Picea glauca (Moench)), trembling aspen (Populus tremuloides Michx.), and eastern white cedar (Thuja occidentalis L.) as the main tree species [28]. The geomorphology is characterized by the presence of a massive clay deposit left by pro-glacial lakes Barlow and Ojibway [29]. The climate is cold temperate with a mean annual temperature of 0.7 • C and mean annual precipitation of 890 mm [30]. The closest meteorological station is located at La Sarre, 42 km north of the study area.

Site Selection
Typical post-fire succession on mesic clay deposits in the eastern Canadian boreal mixedwood forest involves a gradual change from post-fire pioneering stands dominated by shade-intolerant deciduous tree species (trembling aspen or paper birch) during the first ca. 75 years, to mixed stands with an important white spruce component in the next ca. 75 years, to coniferous stands dominated by balsam fir and eastern white cedar after ca. 150 years [28]. It was thus possible, using ecoforestry maps produced by the Quebec Ministry of Forests, Wildlife and Parks [31], to perform an exhaustive fire refuge census by distinguishing, within an otherwise relatively homogeneous forest matrix, stands considered as post-fire residual patches due to contrasting composition and structure representative of late-successional (older) stands. From this preliminary selection [9], thirteen post-fire residual patches were found in areas where the last known fire, reconstructed from dendrochronology, occurred in 1944 or 1923 (depending on site location) [28,32], and the second-to-last fire in 1717 or 1760 [28,32]. These stands have been previously identified [9] as coniferous old-growth forest patches (with balsam fir and eastern white cedar) embedded in a matrix of younger deciduous forests (with trembling aspen or white birch).
From the palaeoecological reconstruction of fire activity based on radiocarbon dating of macroscopic soil charcoal peaks in stratigraphic sections sampled from these post-fire residual patches [9], eight of the 13 patches were identified as fire refuges that had escaped two or more consecutive fires (1923 or 1944 and 1760 or 1717). The five other patches only escaped the most recent fire (1923 or 1944). These were therefore recorded as transient residual patches (Table 1), not used in this study focusing on refuges, and therefore considered as regular coniferous stands thereafter.

Site Selection
Typical post-fire succession on mesic clay deposits in the eastern Canadian boreal mixedwood forest involves a gradual change from post-fire pioneering stands dominated by shade-intolerant deciduous tree species (trembling aspen or paper birch) during the first ca. 75 years, to mixed stands with an important white spruce component in the next ca. 75 years, to coniferous stands dominated by balsam fir and eastern white cedar after ca. 150 years [28]. It was thus possible, using ecoforestry maps produced by the Quebec Ministry of Forests, Wildlife and Parks [31], to perform an exhaustive fire refuge census by distinguishing, within an otherwise relatively homogeneous forest matrix, stands considered as post-fire residual patches due to contrasting composition and structure representative of late-successional (older) stands. From this preliminary selection [9], thirteen post-fire residual patches were found in areas where the last known fire, reconstructed from dendrochronology, occurred in 1944 or 1923 (depending on site location) [28,32], and the second-to-last fire in 1717 or 1760 [28,32]. These stands have been previously identified [9] as coniferous old-growth forest patches (with balsam fir and eastern white cedar) embedded in a matrix of younger deciduous forests (with trembling aspen or white birch).
From the palaeoecological reconstruction of fire activity based on radiocarbon dating of macroscopic soil charcoal peaks in stratigraphic sections sampled from these post-fire residual patches [9], eight of the 13 patches were identified as fire refuges that had escaped two or more consecutive fires (1923 or 1944 and 1760 or 1717). The five other patches only escaped the most recent fire (1923 or 1944). These were therefore recorded as transient residual patches (Table 1), not used in this study focusing on refuges, and therefore considered as regular coniferous stands thereafter.

Field Sampling
Stand composition, structure, and fuel loads were measured in each fire refuge according to the sampling design set by Hély et al. [33], along and in the vicinity of a single 30-m sided equilateral triangle [34]. The same stand sampling methodologies were used to avoid false significant differences in fuel or structures between fire refuges and other matrix stands (previously sampled by Hély et al. [33]) due to changes in sampling method. In each stand, forest structure and canopy characteristics (tree species and diameter at breast height, total height, and canopy base height) were therefore estimated from 24 trees (12 dominant and 12 suppressed) that were selected in the triangle vicinity based on the point-centered quadrant method [34].
Loads of all fuel types defined in the BehavePlus system [22] were measured within each stand along or apart from the triangle sides depending on fuel type [33]. Woody debris from the three American time-lag classes, representative of desiccation times (1 h, 10 h, and 100 h time lags) and corresponding approximately to diameters <0.6 cm, 0.6-2.5, cm and 2.7-7.6 cm, respectively [22], were measured using the line intersect method [35,36]. The same species coefficients and equations as in Hély et al. [37] were used to estimate fuel loads from twig and branch numbers as they had been adapted to the boreal mixedwood forest. Shrub and herbaceous loads were measured in six 1-m 2 quadrats [38], evenly spaced along the 90-m triangle transect. Shrub loads were estimated from basal stem diameter measurements using species dry weight-basal diameter relationships set from shrub samples previously collected in the Duparquet area [39][40][41]. Herbs were collected to obtain oven-dried weight. Litter and duff layer depths were measured in six quadrats (25 cm × 25 cm) and total litter and duff material was separately collected to obtain oven-dried weights.
Stand characteristics of young (deciduous), intermediate (mixed), and old (coniferous) stands representative of the mixedwood boreal forest matrix of the studied landscape were obtained from Hély [42]. They were merged to the stand characteristics of fire refuges ( Table 2) sampled in the present study to create a stand fuel load dataset to be applied to the ecoforestry map polygons ( Figure 1). The few Pinus banksiana stands present in the landscape were classified separately, as these stands represent young but not deciduous stands and were not considered in the analyses.

Weather Data and Associated Fuel Moisture Scenarios
To simulate fire ignition and early fire behavior under different weather conditions, two fire weather indices from the Canadian Fire Weather system [43] were selected. They represent moderate and high fire danger (Fire Weather Index (FWI) = 5 and 15, respectively), as used by the SOPFEU (Quebec Society for the protection of forests against fire) in the studied region. Two days representative of each fire danger were used in order to take into account wind or dry air effects and to partially capture intrinsic weather variability (Table 3). Table 3. Details of the fire weather indices [33]. Note: FFMC, Fine Fuel Moisture Content, which is a numerical rating of the moisture content of litter and other cured fine fuels. This code is an indicator of the relative ease of ignition and flammability of fine fuel. ISI, Initial Spread Index, which is a rating of the expected rate of fire spread. It combines the effects of wind and FFMC on rate of spread without the influence of fuel quantity. BUI, Buildup index, which is a numerical rating of the total amount of fuel available for combustion. FWI, Fire Weather Index, which is a rating of fire intensity that combines ISI and BUI. It is suitable as a general index of fire danger throughout the forested areas of Canada [44].

FWI
These weather characteristics were transformed into fuel moisture content, based on the range of values from scenarios provided by BehavePlus ( Table 4). The fine fuel moisture content (FFMC), duff moisture code (DMC), and drought code (DC) of the Canadian Fire Weather Index [43] were assumed to match with 1 h, 10 h, and 100 h moisture contents, respectively, used in both BehavePlus and FlamMap3 fire models, knowing that dead fuel moisture of extinction (maximum fuel moisture content, which limits fire propagation) is usually set at 30% in BehavePlus. For the calibration process, all stands (broadleaved, mixed, coniferous, and fire refuge stands) were attributed the same water content (1 h, 10 h, and 100 h moistures). Wind direction used in the FlamMap3 model corresponded to the main wind direction (from the south-southwest) over the studied area and is representative of the fire-season wind [45].

Parameterization of Fuel Models at the Stand Level
First, the FBP model was used under the four different weather conditions (Table 3) in order to determine fire behavior variability for each stand type in the landscape mosaic. The FBP System is indeed regarded as producing good predictions of fire behavior when compared to natural or experimental fire records [46][47][48].
Slope and elevation were maintained constant (0 • and 300 m above sea-level, respectively) in FBP and BehavePlus runs to ease comparison and because the study area has an overall flat topography. In the FBP model, stand types are characterized by tree species composition and density, and fuels are qualitatively described [20]. Fire behavior variables (e.g., rate of spread (ROS) and head fire intensity (HFI)) are predicted from empirical relationships computed from many fire measurements recorded during both wildfires and experimental fires, and covering a large range of weather conditions [49]. Moreover, different fire behavior relationships exist for spring (without leaves) and summer (with leaves) stand types that include broadleaved species. Spring fire behavior relationships were selected, because springs sustain faster and more intense fires than summers [33,50]. Spring broadleaved and mixed fuel types (i.e., Dl and Ml, respectively) were therefore chosen, as they were deemed representative of the present boreal mixedwood forest mosaic. The coniferous percentage was increased in M1 stands to age stands toward late successional states dominated by coniferous trees.
The second methodological step involved the BehavePlus model. BehavePlus is a non-spatially explicit deterministic model based on the physical properties and combustion properties of fuel types (see below). The generic fuel model "Moderate load broadleaf litter" available in the BehavePlus system [23] was used for all stand types, replacing fuel load values (for dead (1 h, 10 h, 100 h time lag) and living (herbaceous and woody shrubs) fuels, respectively), fuel bed depth, and canopy structure by those measured in situ ( Table 2). Based on the same topographical (slope and elevation) and weather conditions as for the FBP model, BehavePlus fuel models were calibrated by comparing BehavePlus simulated rates of spread (ROS ) with those obtained from the FBP. Hély et al. [48] have concluded that systematically slower ROS were simulated by Behave (earlier version of BehavePlus) as compared to the ROS simulated by FBP, and that these differences were likely due to the exclusion of the duff layer in Behave simulations. We addressed and solved this problem by adding the duff load to the 1-h fuel load. Consequently, the fuel bed depth was also increased in the calibration process until BehavePlus ROS predictions were in the same range as those from the FBP predictions with still realistic fuel bed depths in the range of observed values (Table 1).

Settings for Simulations at the Landscape Level
Spatially explicit fire simulations at the landscape mosaic level were performed using the FlamMap3 model [24,25]. FlamMap3 is based on BehavePlus fuel models and fire behavior at the pixel level. FlamMap3 propagates fire from an ignition location (randomly selected or not) through neighboring pixels over a given simulation time (see below). Potential fire behavior calculations include surface ROS [51], crown fire initiation [52], and crown ROS [53]. The ignition-propagation scheme is repeated a given number of times set by the user (with random ignition locations). FlamMap3 provides an output map for each fire behavior variable (e.g., ROS and Head Fire Intensity (HFI)). It also provides a burn probability map reporting for each pixel the number of times fire went through it compared to the total number of ignitions.
Slope ( • ), elevation (m above sea-level), aspect ( • N), and composition maps for the Lake Duparquet landscape, including fire refuges, were rasterized (225 m 2 spatial resolution, i.e., 15 × 15 m) using the ArcGIS software to produce Ascii grids required as input data in FlamMap3. A landscape fuel map was created by applying to the stand composition ecoforestry map the four fuel models representative of mean broadleaved, mixed, coniferous, and Pinus banksiana stands, respectively. A fifth fuel model, representative of mean fuel conditions measured in fire refuges, was specifically defined and projected on the fire refuge locations. Water bodies and recently logged stands were represented as a generic non-fuel model and considered as firebreaks. Moreover, as for BehavePlus at the stand level, for each fuel model (stand type), fuel load (kg/ha) ( Table 2), canopy cover (%), tree height (m), crown bulk density (kg/m 3 ), and height-to-live crown base (m) were used to create the FlamMap3 required maps based on stand attributes.
In order to simulate fires, whose size distribution would express the natural variability recorded in the boreal mixedwood forest [27], four simulation times in FlamMap3 (63, 105, 408, and 4835 min) were selected, corresponding to the four quartiles (25%, 50%, 75%, 100%) of the fire size distribution (Figure 2) computed from Quebec public archives (1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007) provided by the SOPFEU [54]. Simulation times were calculated using the mean pixel-based ROS computed from preliminary FlamMap3 output run over the entire landscape. However, as the smallest simulated fires (63-min fires) could not cross several pixels or reach fire refuge stands (except when ignition occurred within), results were only reported for simulations with fires whose size was at least equal to the median archived fire size in Quebec. Hereafter, the median size fires were called small fires, the 75% size fires were called medium fires, and the 100% size fires were called large fires.

Modeling Experiments
Five modeling experiments were performed to test the effect of different stand characteristics and their combinations on the occurrence of fire refuges. Each experimental design produced a burn probability map.
The first experiment (Table A1) tested the effect of observed stand attributes (fuel type proportions, loads, and stand structure) combined with the indirect effect of stand location (the occurrence of fire breaks in the vicinity, topography). A given fuel moisture scenario was applied to all stands throughout the landscape. This experiment included a total 120,000 simulations representing 12 runs of 10,000 ignitions, each run being a combination of one simulation time (over three) and one weather/moisture scenario (over four) with 10,000 random ignition locations.
The second experiment (Table A1) tested the effect of low occurrence of fire refuges in the landscape. To proceed, the number of fire refuges was artificially increased to 800 stands and randomly placed in the landscape to represent 2% of the forest mosaic area. To save time, computer space, as well as to keep high computer performance, only 40,000 randomly ignited fires were simulated, based on the four weather-moisture scenarios and the 408-min simulation time (75th percentile) only. The potential differences in burn probabilities computed with this second experiment and those extracted from the previous simulations with fire refuges in their real location but only with the corresponding 408-min simulation time could be attributed to their low occurrence in the landscape.
The third experiment (Table A1) tested the combined effect of firebreaks and topography on the occurrence and persistence of fire refuges. Fire refuges were manually placed in the vicinity of water bodies and in particular topographical situations (depressions) to test their effect on the probability of burning. In this experiment, as for the previous one, only 40,000 fire ignitions were performed using only the 408 min simulation time.
In the fourth experiment (Table A1), the fire refuge microtopography was tested (stand scale), as the flat macrotopography (landscape scale) was not expected to explain the location of fire refuges. Indeed, the thick organic matter layer previously measured in most of the fire refuges [9] could reflect the presence of shallow topographical depressions created just after the drainage of Holocene proglacial lake Obijway-Barlow. These depressions would have been filled with organic matter during the Holocene. This assumption comes from the fact that topography at ground surface in fire refuges did not differ from the surrounding forest matrix while fire refuges had significantly thicker organic matter accumulation [55]. Therefore, in order to reproduce the initial conditions in fire refuges (without organic matter accumulation), the presence of small depressions was tested by artificially changing the altitudinal conditions of fire refuges. As for the two last experiments, only 40,000 new randomly ignited fires were simulated (with 10,000 per moisture scenario based on the 408 min simulation time).

Modeling Experiments
Five modeling experiments were performed to test the effect of different stand characteristics and their combinations on the occurrence of fire refuges. Each experimental design produced a burn probability map.
The first experiment (Table A1) tested the effect of observed stand attributes (fuel type proportions, loads, and stand structure) combined with the indirect effect of stand location (the occurrence of fire breaks in the vicinity, topography). A given fuel moisture scenario was applied to all stands throughout the landscape. This experiment included a total 120,000 simulations representing 12 runs of 10,000 ignitions, each run being a combination of one simulation time (over three) and one weather/moisture scenario (over four) with 10,000 random ignition locations.
The second experiment (Table A1) tested the effect of low occurrence of fire refuges in the landscape. To proceed, the number of fire refuges was artificially increased to 800 stands and randomly placed in the landscape to represent 2% of the forest mosaic area. To save time, computer space, as well as to keep high computer performance, only 40,000 randomly ignited fires were simulated, based on the four weather-moisture scenarios and the 408-min simulation time (75th percentile) only. The potential differences in burn probabilities computed with this second experiment and those extracted from the previous simulations with fire refuges in their real location but only with the corresponding 408-min simulation time could be attributed to their low occurrence in the landscape.
The third experiment (Table A1) tested the combined effect of firebreaks and topography on the occurrence and persistence of fire refuges. Fire refuges were manually placed in the vicinity of water bodies and in particular topographical situations (depressions) to test their effect on the probability of burning. In this experiment, as for the previous one, only 40,000 fire ignitions were performed using only the 408 min simulation time.
In the fourth experiment (Table A1), the fire refuge microtopography was tested (stand scale), as the flat macrotopography (landscape scale) was not expected to explain the location of fire refuges. Indeed, the thick organic matter layer previously measured in most of the fire refuges [9] could reflect the presence of shallow topographical depressions created just after the drainage of Holocene proglacial lake Obijway-Barlow. These depressions would have been filled with organic matter during the Holocene. This assumption comes from the fact that topography at ground surface in fire refuges did not differ from the surrounding forest matrix while fire refuges had significantly thicker organic matter accumulation [55]. Therefore, in order to reproduce the initial conditions in fire refuges (without organic matter accumulation), the presence of small depressions was tested by artificially changing the altitudinal conditions of fire refuges. As for the two last experiments, only 40,000 new randomly ignited fires were simulated (with 10,000 per moisture scenario based on the 408 min simulation time).
Finally, as the thick organic matter in fire refuges has been assimilated to peat [55] and therefore very moist ground conditions, the fifth experiment (Table A1) tested the effect of fuel moisture on the probability of burning. To proceed, the moisture content was artificially increased for fire refuges only, without exceeding the fuel moisture of extinction, and without changing fuel moisture in the forest matrix (kept constant for a given weather/moisture scenario as in the other experiments). We ended with fire refuge fuel moistures of 21%, 24%, and 26% (for 1 h, 10 h, and 100 h dead fuel, respectively). Once again, only 40,000 randomly ignited fires were simulated based on the median 408 min simulation time and the four general weather-moisture scenarios used in the other experiments.

Statistical Analyses
All analyses were performed using the R software [56]. The analyses were performed not directly on the burn probability output from FlamMap3, but rather on stand propensity to burn as defined as the ratio of cumulative area burned probability for each stand type (burned proportion) on its representativeness (stand proportion in the forest mosaic matrix). For a given stand type, a ratio higher (lower) than 1 highlights a relative propensity to burn or to propagate fire (to escape fire or to slow down fire spread). This was necessary due to the fact that fire refuges represent only a few stands in the overall studied landscape and likely are even less abundant in the simulated area burned.
For each modeling experiment, stand type propensities to burn were tested using Chi-square and contrast tests in order to compare their cumulated burned proportion (observed values) to their cover representativeness in the forest matrix (theoretical percentage values). Second, we looked for significant differences in burn probabilities among forest stand types (i.e., among fuel model types) using one-way non-parametric analyses of variance on rank scores (Kruskal-Wallis non-parametric test) followed by a Tukey multiple comparison test.
Using the simulation dataset from the first experiment, we tested the effect of stand attributes (fuel type proportions, loads, and stand structure) including as well the combined effects of macro-topography and firebreak locations (experiment 1). The similar overall response through weather scenarios and through simulation times helped us in restricting the last four simulation experiments to the 75 percentile simulation time of 408 min.
The number of fire refuges (experiment 2), their location as compared to firebreaks (experiment 3), their maintenance over micro-topographic depression buildup (experiment 4), as well as their fuel moisture content (experiment 5) were tested using the same statistical approach, but only applied to the series of 408 min runs with the four fuel moisture scenarios.

Results
At the end of the parameterization procedure, ROS predicted by the FBP and BehavePlus systems were in good agreement (R 2 = 0.8), ranging from 0.2 to 9 m/min. A slight overestimation from BehavePlus in conifer stands was noted, with the heaviest fuel loads under highest fire risk with windy conditions. This satisfactory agreement between both fire behavior models running at the stand level confirmed that the BehavePlus model could be adapted to the boreal mixedwood forest if the fermentation layer (deep litter) was taken into account in fuel load characterization. This adjustment of fuel load yielded more realistic values for all stand types.
While the large number of iterations (10,000 fire ignitions) for each scenario and experiment almost always resulted in rejection of the null hypothesis, we chose to interpret only the most significant signals in stand propensity to burn with a special focus on fire refuge type. We also chose to interpret only the general pattern among stand types in order to evaluate the ecological relevance of the landscape fire model, discarding minor differences which were difficult to interpret from an ecological standpoint.
Simulation comparisons from the first modeling experiment based on actual forest mosaic stand characteristics (location, number, and fuel types) but with set weather/fuel moisture scenarios, showed, for each fire size class, that fire refuges seem to have a high propensity to burn (up to 3.46), except for scenario 4 (Figure 3). Moreover, the relative propensity to burn of fire refuges seems to be higher than those of any of the other stand types (Table 5). Between-stand differences decreased with increased fire danger (scenarios 3 and 4) and fire duration. . Propensity of stand types to burn as a function of fire danger (scenarios) and fire duration (a proxy for small, medium, and large fires based on the 50%, 75%, and 100% quartiles of the burned area distribution from the 1994-2004 SOPFEU fire database [54]). In each of the twelve panels, we tested the stand type propensity to burn as compared to its landscape representativeness using Chi-square and contrast tests. * for significant departures from 1 (burning equal to representativeness, dashed line).

Figure 3.
Propensity of stand types to burn as a function of fire danger (scenarios) and fire duration (a proxy for small, medium, and large fires based on the 50%, 75%, and 100% quartiles of the burned area distribution from the 1994-2004 SOPFEU fire database [54]). In each of the twelve panels, we tested the stand type propensity to burn as compared to its landscape representativeness using Chi-square and contrast tests. * for significant departures from 1 (burning equal to representativeness, dashed line). Table 5. Comparisons of burn probabilities among stand types, based on actual mosaic composition and spatial arrangement, as a function of fire danger (weather/fuel moisture scenarios) and fire duration (a proxy for small, medium, and large fires based on the 50%, 75%, and 100% quartiles of the burned area distribution from the 1994-2004 SOPFEU fire database [54] 408 min (75%) The fifth modeling experiment ( Figure 4, Table 6), testing the effect of very moist fuels in fire refuges compared to other stands, was the only experiment that systematically showed significant reduction in fire refuge propensity to burn through all fire danger scenarios (Figure 4, right panel). The reduction was so strong that the fire refuge propensity to burn shifted to a propensity not to burn (ratio < 1) and fire refuge stands were systematically ranked in third position, just above deciduous stands but below mixed and coniferous stands. For all other modeling experiments, propensity to burn generally increased from broadleaved to coniferous stands, and those from fire refuges were among the highest (Figure 4, first three columns and Table 6). Differences in propensity to burn of fire refuges compared to those from their actual locations (second experiment in Figure 4 versus Figure 3) were likely due to the random ignition effect. Table 6. Comparisons of burn probabilities among stand types as a function of experiments using only the medium fire duration (408 min runs) representing the 75% quartile of the burned area distribution from the 1994-2004 SOPFEU fire database [54]. Significant differences from Kruskal-Wallis and Tukey tests (with p < 0.05) among stand types are represented by different letters.

Number of refuges
Topography and firebreaks High fuel moisture Number of fire refuges experiment Topography-firebreaks experiment Small depression experiment Higher fire refuge fuel moisture experiment Figure 4. Propensity of stand types to burn as a function of experiments using only the medium fire duration (408 min runs) representing the 75% quartile of the burned area distribution from the 1994-2004 SOPFEU fire database [54]. In each of the twelve panels, we tested the stand type propensity to burn as compared to its landscape representativeness using Chi-square and contrast tests. * for significant departures from 1 (burning equal to representativeness dashed line). The propensity to burn of fire refuges in scenarios 1 and 3 of the topography-firebreak experiment reached the value of 8 (although we topped the y-axis at 4 to ease readability).

Discussion
Results from the present study are based on a large number of fire behavior simulations at the landscape level, combining different fire sizes and weather/moisture scenarios, and testing different environmental factors possibly explaining the occurrence of fire refuges, a peculiar type of post-fire residual patches [18,33,57,58]. Moreover, due to the influence of environmental features such as fire severity, location in the landscape, and fuel characteristics, fire refuges, unlike other transient residual patches, are assumed to not be randomly distributed [59].
Although the large number of fire ignitions almost systematically pushed the statistical interpretation toward rejection of the null hypothesis, we clearly and logically found an increase in propensity to burn from broadleaved to coniferous stands, suggesting that FlamMap3 is an efficient model to study the relationship between vegetation and fire in the mixedwood boreal forest. Our analysis also suggests that stand fire hazard increases through the successional sequence in canopy tree species replacement [33]. This increase in propensity to burn only results here from differential fuel load accumulation and fuel spatial arrangement (input in FlamMap3) [27,33], but not from changes in surface fuel quality (composition) since it has not been directly considered as input in the FlamMap system. The late successional coniferous forest stands accumulate fuels, leading to increased forest combustibility known as build-up [4,60,61]. This trend was confirmed in most simulation results except for three simulations whose patterns could only be explained by the randomness of ignition locations (Figure 4), which were likely overrepresented in some types (broadleaved, refuges, or even in water bodies and human disturbed areas).
A second logical result from the real forest mosaic composition (i.e., stand proportions) was the effect of weather conditions, which induced a decrease in stand propensity to burn with increasing fire weather risk, as well as with increasing fire duration (i.e., fire size). Indeed, both changes allowed fires to spread more easily and therefore increased the probability to reach many stand types.
Regarding the main objective of the present study, we showed that fire refuges in the eastern Canadian mixedwood boreal forest burned more than the surrounding forest matrix when fuel moisture per fuel type (1 h, 10 h, 100 h, litter, and live fuels, respectively) was held constant throughout the landscape. From this first analysis, it is clear that in their current location, the association of fuel types and loads, topography at ground level, and firebreaks in their vicinity cannot explain the presence of fire refuges nor their spatial distribution in the forest mosaic. This was confirmed by the results of the second experiment with more numerous fire refuges randomly located in the landscape, as fire refuges maintained a higher propensity to burn than the surrounding forest matrix. Hence, fire refuges in the Lake Duparquet area actually present fuel characteristics favorable to fire ignition and fire spread in terms of shape, size, density, loading, chemical properties, and spatial configuration [62], but not moisture.
Similarly, our results excluded the predominant effects of firebreaks and topography in the occurrence of refuges based on the first, third, and the fourth modeling experiments. This result contrasts with those from previous studies explaining that the spatial occurrence of post-fire residual stands was influenced by firebreaks, such as rock outcrops and water bodies, which disrupt horizontal fuel continuity, thus preventing fire propagation [3,17,18,55]. Despite the relatively high proportion of water bodies in the study area, their role as firebreaks does not seem to explain the lower propensity to burn of fire refuges. However, this study cannot exclude the effect of firebreaks on other residual patches (transient residual patches) in the landscape [63]. Indeed, the role of lakes intercepting and stopping fire spread more likely depends on wind direction during fire and on ignition location in the landscape in comparison with residual patches location. The importance of protected topographical positions was also reported as an important factor in interior forests of eastern Washington [64] and south central Wyoming [65], where old forests are found. The generally flat topography in our study area as compared to the abovementioned regions likely explains why no link between the occurrence of fire refuges and landforms was found here. Moreover, while we tested for the presence of Holocene small depressions as initial conditions favorable to fire refuge creation, they appeared to have no negative effect on fire behavior and therefore cannot fully explain the occurrence of these persistent unburned stands in the landscape. However, fire refuges burned less than other stand types (except pure broadleaved stands) when they were assigned higher fuel moisture content in the simulations as observed and measured in the field. Because fire refuges escaped multiple stand-replacing fires that have occurred in the surrounding forest matrix, high fuel moisture appears to be the critical factor reducing fire intensity, as also reported in previous studies [63,66]. Moreover, it is worth noting that under past extreme fire weather conditions, even fire refuges sometimes have burned [9]. Beyond weather condition, fuel moisture is usually related to topography and aspect [67]. The aspect determines the solar-flux (cool, wet north and east facing aspects) and can impact soil moisture, which has important influences on fire behavior [48,68,69]. In the present study, due to the flat topography, the high moisture content of the litter in fire refuges seems to be rather related to the indirect effect of shallow depressions that have been filled with organic matter though time (centuries to millennia [9]), in which water tends to accumulate and peat to develop.
While the results of this study seem to be in accordance with palaeoecological analyses performed in situ, comparisons with other studies or regions are difficult. First, differences between fire refuges and other transient residual patches have never been tested before. Second, the role of firebreaks and topography can vary regionally based on the macro-and micro-topography of the study area. Therefore, to understand the general pattern of fire refuge occurrences in a given landscape, more investigations are needed locally to take into account the present day mosaic specificities and potential past changes as revealed by palaeoecological analyses. Hence, given the potential importance of fire refuges in the landscape (biodiversity hot spots [10,11]), they should be subjected to special conservation efforts.

Conclusions
To conclude, fuel moisture content appears to be the most important factor influencing the distribution of fire refuges at the landscape scale in the Eastern Canadian mixedwood boreal forest where topography on clay belt is relatively flat. This result is in good agreement with palaeoecological analyses performed in the same stands, which showed the occurrence of aquatic taxa and moisture tolerant tree species in fire refuges [70]. Hence, fire refuges could be considered as "soaked powder kegs", having well enough fuel to burn, but too much humidity. As dryer climate conditions are expected in the northeastern North American boreal forest over the next decades [71], fire conditions leading to the burning of fire refuges could become more frequent. Simulation studies using various climate change scenarios are necessary to evaluate potential effects on the persistence of fire refuges.