Are fire refugia less predictable due to climate change?

Fire refugia—unburnt habitat within a wildfire’s perimeter—play a key role in wildlife persistence and recovery. While studies have shown that the location of refugia is influenced by local topographic factors, growing evidence points to extreme fire weather becoming the dominant factor driving high-severity wildfires that result in the location of fire refugia being less predictable. Between September 2019 and February 2020, a series of mega-fires in eastern Australia burned largely in broadleaf forest. We assessed burned and unburned areas of forest in eastern Australia using Sentinel-2 satellite data, aggregated monthly over the fire season to calculate a fire severity layer at a 20 m pixel resolution. We found that fires burned 5.7 × 106 ha−1 of forest and woodland. The total percentage area of unburned forest within the wildfire footprint was approximately 10%. The majority (94%) of the unburnt forest and woodland patches within the fire perimeter occurred as patches <1 ha (n = 842 622 and 111 707 ha) with far fewer large unburnt patches (>100 ha) (n = 575 and 286 080 ha). Boosted regression tree analyses of the relationships between fire severity and potential explanatory variables revealed that 63%–78% of the variable importance in the models were climatic and weather-related factors. Fire weather index was the single most important variable for analyses, accounting for 40%–52% of modelled results. Our results reinforce mounting evidence that a shift is underway in the balance between deterministic and contingent factors in the occurrence of fire refugia with local topographic controls being increasingly overridden by severe fire weather conditions, and declining topographic effects as fire severity increases. Further studies are needed over a longer time frame, inclusive of prior forest management impacts, to confirm that the ability to predict fire refugia is permanently declining.


Introduction
The role of fire refugia-unburnt habitat within the perimeter of a wildfire-in the persistence and recovery of wildlife populations is well recognized [1,2]. Their significance for conservation management is escalating with the increasing prevalence and extent of forest fires, particularly in temperate and boreal biomes subject to mega-wildfires (>1.0 × 10 6 ha) [3][4][5][6]. In addition to protecting animals from being burned and killed, refugia provide vegetation-based habitat resources for food, shelter and reproduction that can support population recovery. Vegetation structure is critical in providing suitable habitat for many forest vertebrates and invertebrates and fire is therefore a major driver of variation in habitat structure, resource availability and wildlife demography in many terrestrial ecosystems [7][8][9].
The capacity to manage fire refugia is largely constrained by the extent to which deterministic factors enable their occurrence in the landscape to be accurately predicted and therefore integrated into conservation management plans. However, if contingent factors related to fire weather conditions dominate, then the occurrence of refugia will be random or haphazard, with implications for wildlife recovery and conservation. A likely key driving factor in this shift is more severe fire weather conditions due to climate change [6,10,11], increasing the likelihood of contingent factors dominating.
In addition to fire weather conditions, the occurrence of forest fire refugia has been found to be related to factors such as topography and fuel characteristics (types, loads and dryness) [12]. Topographic factors include topographic shading, as well as surface and subsurface water run-on areas, which result in cooler and more moist site conditions, enhanced plant water availability and less combustible fuel loads [2]. Where the occurrence of refugia is strongly influenced by these topographic factors, refugia can persist through multiple fires and their location can be predicted in space and time [13]. Fuel characteristics vary with major vegetation types [14], with the canopy of closed forest generating micro-climates that can limit the movement of fires. Deterministic factors affecting the location of fire refugia can include topographic position and forest ecosystem types [12,15]. However, extreme fire weather conditions can drive high-severity wildfires with impacts that override the local topographic-related determinants [16], potentially resulting in the location of refugia being more contingent, ephemeral and unpredictable.
Between September 2019 and February 2020, approximately 5.7 × 10 6 ha of mainly temperate broadleaf forest was burned in eastern Australia by a series of mega-fires [17]. To date, research and commentary on the impacts of these fires have focused on the area of forest burned, fire severity, the immediate impacts on wildlife, attribution to climate change and extreme fire weather, and wider societal impacts [10,18,19]. However, information on the location of unburned forest within the fire perimeter, together with understanding their potential to serve as refugia and their role in post fire recovery, is critically important for conservation.
The 2019 and 2020 mega-fires presented a unique opportunity to test whether the drivers of refugia have changed from being dominated by topographic factors to being driven by weather and climate. Our underlying null hypothesis was that no factors affect the location of refugia, i.e. their distribution is random. We also tested two hypotheses. Hypothesis 1 was that topography remains the dominant determinate of the location of refugia as reported previously [2,13]. Hypothesis 2 was that weather and climate now dominate the location of fire refugia. If the null hypothesis and hypothesis 1 are rejected, and hypothesis 2 accepted, this would mean that the location of fire refugia are more likely than not to be contingent on fire weather conditions with major implications for conservation planning and management.

Methods
We provide the workflow for the calculation of the burn severity index classes and related statistical analyses in figure 1. Further details on workflow steps are provided in supplementary information (available online at stacks.iop.org/ERL/16/114028/ mmedia). We defined fire refugia as unburned patches within the fire perimeter. We assessed burned and unburned areas of forests using Sentinel-2 satellite data for the period 1 September 2019 to 28 February 2020, aggregated monthly over the fire season using Google Earth Engine [20]. We created monthly fire perimeters using the differenced normalized burn ratio (∆NBR) and cross validated with data from visible infrared imaging radiometer suite (VIIRS) and moderate resolution imaging spectroradiometer (MODIS) fire detection products to remove incorrectly allocated areas [21,22].
We mapped burn severity within the fire perimeter using previously established differenced normalized burn ratio classes thresholds [23] which have been found to be accurate within temperate Australia [24] and are similar to those established by Key and Benson [25] and employed globally to assess fire severity across varying ecosystem and forest types [24,[26][27][28]. The lower threshold was assessed for each month and region and increased to between 0.11 and 0.13 to minimize the impact of rapidly drying vegetation to increase the accuracy unburnt/low burn severity threshold. We classified areas of burnt forest into five burn severity classes: unburnt (<0.11-0.13), low (>0.11-0.13), moderate-low (>0.27), moderate-high (>0.27) and high (>0.44).
We evaluated the burn severity layers using GIS spatial analysis tools to map forest patches of contiguous cells within each of the five severity classes. The burn severity classification from each month was derived from the highest burn severity per cell within the perimeter. We calculated patch statistics based on their area and location and with respect to forest type and land tenure. We stratified our analyses by the Interim Biogeographic Regionalisation for Australia (IBRA), a nationally recognized classification and map of bioregions and subregions (IBRA) [29] (figure S1).
We used boosted regression trees using the gbm package in R to examine the relationships between fire severity class and a candidate set of potential explanatory variables. We constructed a series of boosted regression tree models to test which environmental factors had the greatest influence on fire severity. This machine learning technique has been found to cope with complex non-linear relationships through creating large numbers of regression trees. Numerous studies assessing the factors influencing fire and fire severity have found that boosted regression tree analysis to be a useful method [30][31][32]. To complete the analysis, we generated a one-hectare resolution grid across the entire fire perimeter, within forests and woodlands, resulting in a total of 71 32 452 cells. Additionally, we created polygons of unburnt and high severity patches using contiguous cells within the fire perimeter. Mean and mode values were then calculated for each variable per grid cell and patch using zonal statistics which calculates statistics of a raster layer for each feature of an overlapping polygon vector layer. We implemented the boosted regression tree model using the method established by Elith et al [33], with a learning rate of 0.05 and a tree complexity of 'five' . We conducted model training on 5% (356 620 grid cells) of the data and we tested twice from random 5% subsets of the remaining 95% of the data.
We chose the candidate set of potential explanatory variables based on prior research that had established their significance in understanding forest wildfire behaviour and severity [15,31,34]. We had to generate, filter or adapt, a number of the variables for the analysis (supplementary information). The candidate explanatory variables were: fire weather index (FWI); forest fire danger index (FFDI); mean monthly total shortwave radiation; aspect; topographic wetness index (TWI); topographic position index (TPI); terrain ruggedness index; forest types; visible atmospheric resistant index (VARI); and normalized difference infrared index band 6 (NDIIb6). The FWI and FFDI are generated from interpolated weather station data on humidity, temperature, wind speed and precipitation. The FFDI also including a drought factor based accumulated soil moisture deficit. Due to the similarities between the two indices, we tested them in separate models. To avoid collinearity in potential explanatory variables, we tested all using a correlation matrix, ensuring variables had a low or moderate Pearson's pairwise correlation (r < 0.7) [35]. We provide definitions and data sources for these variables in table 1. We completed a set of analyses where the spatial unit of analysis were either 10 4 m 2 grid cells or polygons. For the grid-based analysis, the independent variable was continuous values of the differenced ∆NBR. For the polygon-based analyses, the independent variables were patches of contiguous cells assigned to either unburnt or high fire severity class.

Results
Our analysis of the 2019 and 2020 wildfires showed that the mega-fires burned 7.04 × 10 6 ha of land spanning a 1160 km stretch from south east Queensland to eastern Victoria, of which 5.7 × 10 6 ha was forest and woodland (figure 2 and table 2). This is approximately 30% of the total forest area within the bioregions impacted (table 3). The wildfires burned in all forest types, bioregions and land tenures, including protected areas and forests managed for commodity (wood) production, located within the fire perimeter (figure 2, tables 2 and 3). There were substantial Table 1. The independent variables in the boosted regression tree analysis. Also included is the source dataset (which was either directly used or modified for the analysis; see supplementary information), as well as a description of each data layer.

Variables
Source dataset The National Vegetation Information System is a standardized spatial layer of vegetation types throughout the country. The major vegetation groups layer was chosen for the analysis. VARI -500 m VARI used to assess vegetation moisture and has been found to be important for measuring fuel load (34). The index was created using MODIS satellite imagery and Google Earth Engine and with imagery 1 month prior to the fire being normalized to data from 5 years prior. NDIIb6 -500 m The NDIIb6 is a vegetation moisture index used to measure drought. The index was created using MODIS satellite imagery and Google Earth Engine and with imagery 1 month prior to the fire being normalized to data from 5 years prior.
areas (777 378 ha) of high severity fire in all bioregions, forest types and tenures (figure 3, tables 2 and 3). Natural World Heritages sites in these bioregions also were impacted, with 16% of the 151 662 ha of rainforest in the Gondwanan Rainforest of Australia World Heritage Area burnt. Approximately 34% of impacted forest was burnt at a medium-high to high level of severity (table 3). The total percentage area of unburned forest within the wildfire perimeter was approximately 10% (tables 2 and 3; additional statistics provided in table S1). Eucalypt woodlands (forest with 10%-30% canopy foliage cover) had the smallest percentage remaining unburnt area of forest types (6%) and experienced the largest proportion of high severity fires (18%; 201 469 ha). While rainforest was the least impacted of the broad forest types assessed, the impact was still considerable with 163 620 ha burnt (51% at a low severity, 22% of medium-high and high severity) and only 30% unburnt. The majority (94%) of the unburnt forest and woodland patches within the fire perimeter occurred as patches <1 ha (n = 842 622) (figure 4) with far fewer large unburnt patches (>100 ha) (n = 575) (table 4). We found that three variables-FWI, forest type and the vegetation moisture index (NDIIb6)-most strongly influenced burn severity (∆NBR), collectively accounting for 70%-79% of all three boosted regression tree models. Higher burn severity was associated with increasing FWI values and decreasing NDIIb6 values. Low fire severity was also associated with rainforest forest type. FWI was the single most important variable for analyses, accounting for 42%-53% of the relative importance in the model of each variable. The next important variable was forest type (15%-26%) and NDIIb6 with 8%-11%. Even when The location of fire refugia patches was not associated with local topographic factors as represented by the candidate set of potential explanatory variables in the boosted regression tree models. Neither was there any discernible variation from visual inspection of the kernel densities in the distribution of fire refugia in relation to the topographic gradients, with unburnt, low, moderate and higher severities found across both the TWI and the TPI gradients ( figure 6). The results suggest that both small (<1 ha) and large (>100 ha) fire refugia patches occurred across topographic gradients and land tenures.

Discussion
The results of our boosted regression tree analyses revealed that the occurrence of refugia was strongly correlated with climatic and weather-related factors (the FWI, the FFDI, and the NDIIb6 vegetation moisture index) and that the topographically-related factors played only a minor role in influencing where refugia occurred in the landscape. The FWI was the single most important variable for analyses, accounting for 42%-53% of modelled results. Previous studies found strong correlations between the spatial distribution of wildfire refugia and local topographicrelated factors [2,12,15,36].
Relative to the scale of the wildfires, we found that large fire refugia were uncommon, and the occurrence of both small and large patches were spread across bioregions, forest types, and land tenures and did not appear to be correlated with local topographic gradients. The dangerous fire weather conditions driving the south-eastern Australia mega-fires [11] resulted in the location of fire refugia within the fire perimeter being more contingent on stochastic daily fire weather conditions than deterministically based on local environmental conditions largely related to topography. Forest type was one of three important explanatory variables in the boosted regression tree models which helps account for the range in the proportion of unburnt areas (15%-26%). Forest type can exert a strong control on the location of wildfire refugia due to variation in fuel loads and moisture levels [37]. Severe fires in 2019 and 2020 occurred in all forest types, including in rainforest and riparian forest that is typically fire-resistant.
When the location of fire refugia in a forest landscape is largely determined by local environmental factors-topography and forest type-then it is feasible to plan for their conservation by protecting them from inappropriate land management practices such as post-fire logging [38] and inappropriate applications of fuel reduction burning [39,40]. Should the south-eastern Australia mega-fires prove to be a single event, unlikely to be repeated, then possibly no change is required in forest management beyond dealing with the current impacts. However, if the fire regime has undergone a phase shift to a new state with permanently enhanced fire risk with a haphazard distribution of fire refugia, then new forest management strategies would be required including for wildlife conservation. For fire refugia contingent on prevailing weather conditions and largely de-coupled from locally scaled environmental determinants, it is not possible to plan for their conservation and provide guidance by prescribing for specific locations what forest management practices are permitted or should be excluded to help maintain their wildlife habitat value.
The 2019/20 Black Summer bushfire season in south-eastern Australia was unprecedented in terms of the area of forest burnt, the types of forest burnt, and the number of fires that developed into extreme pyroconvective events [11]. A dry winter followed by hot dry windy conditions in September, with litter moisture content across eastern Australian forests at record low levels, providing the preconditions for such major wildfires [10,41]. The evidence for how these extreme fire weather conditions are related to climate change is mounting. The 1 • C of warming since pre-industrial levels is associated with longterm change in fire weather conditions since the 1950s, including an increase in the FFDI during spring and summer in southern Australia [42]. Attribution studies focused on the eastern Australia megafires have established that some drivers of FFDI show an imprint of anthropogenic climate change Table 3.     See table 1 for details about each variable. Note that the partial dependency plot does not consider interactions each variable has with another. Table 5. Diagnostics from the boosted regression trees of relationships between fire severity class and potential explanatory (table 1) [11,43]. South-eastern Australia is one of a number of regions globally that have experienced an increase in dangerous fire weather conditions since the 1970's [4,6,42,44] and these elevated conditions are projected to either be maintained or continue to escalate further, depending on the world's success or failure to mitigate greenhouse gas emissions [45,46]. Two critical forest policy and management issues arise in the wake of the mega-fires and projected climate impacts. These are: (a) the ongoing management of the wildfire refugia from the 2019 and 2020 fires, and (b) planning for future fire refugia. The fires had a massive impact on wildlife with 119 animal species requiring urgent conservation management intervention [47] and some ecosystems, particularly relictual Gondwanan rainforests, susceptible to regeneration failure [48]. Careful management of fire refugia will be critical given their importance for post-fire recovery especially given climate change trends and land use pressures [40,[49][50][51]. Both small and large fire refugia patches will be important, as noted in other studies in temperate forests [12,15,16].
The refugial values of burned areas must also be considered in supporting in situ survival of residual populations. For example, a study of two native mammal species-agile antechinus (Antechinus agilis) and bush rat (Rattus fuscipes)-in eucalypt forests of south eastern Australia found there were still present in burnt sites though at 30% and 12% of the density observed in unburnt sites, with in situ survival, and not recolonisation, being the most plausible explanation [52]. Less severely burnt areas will provide greater protection for a wider range of species and retain more habitat resources and recover more quickly, enhancing species dispersal and recolonisation [52]. Re-consideration is therefore warranted of current forest management approaches on all land tenures including state forests and national parks. Existing refugia should not be subject to further human caused disturbances to help ensure they maintain their wildlife habitat value. Given the trend in the increasing difficulty of predicting the location of future refugia, other disturbances in the landscape that are predictable, and which disrupt natural processes of forest regeneration and reduce forest resilience should be limited or entirely withdrawn wherever possible [49].
A limitation of this study was the relatively coarse temporal resolution of the fire weather data compared with previous studies which had access to hourly fire weather data and spatially fine grained topographic data for analysing the fire-front [53]. We used the best available data for the 182-d wildfire event that burned 5.7 million ha of forest. Further research is need into the accuracy of the ∆NBR burn severity index and derived classes given that for many species, low and unburnt areas may not be functionally different. Validation analysis from a study of wildfire impacts on southern Australian Eucalypt forest found an overall accuracy of 85% with 98% for the high burn severity class, 80% for the low class and 75% for the unburnt class [54]. Similar studies are need to confirm if these results hold across a wider range of forest types, including rainforest and vine thicket, and how the error for low and unburnt classes can be reduced. Research is also needed to address the potential uncertainty associated with the high end of ∆BNR if the signal saturates and is unable to distinguish ecological significant thresholds. In particular, the difference between a severe 'crown burn' and total 'crown scorch' is important for keystone forest canopy species that regenerate only from seed, with adults usually being killed by the latter [55].

Conclusion
Our results reinforce mounting evidence that a shift is underway in the balance between deterministic and contingent factors in the occurrence of fire refugia with local topographic controls being increasingly overridden by extreme fire weather conditions and declining topographic effects where burn severity increases [13,24]. Given the currently exacerbated and projected increases in dangerous fire weather conditions [11,42,45], it is more likely than not that fire refugia will increasingly be the result of temporally contingent rather than spatially deterministic factors, placing further stress on these forests' biodiversity. Further studies are needed over a longer time frame, inclusive of prior forest management impacts, to confirm that the ability to predict fire refugia is permanently declining.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.

Source code
The computer code associated with all main analyses carried out during the study can be found in the GitHub repository (https://github.com/patrickm-norman/Fire_refugia_study).