Influences of climate fluctuations on northeastern North America’s burned areas largely outweigh those of European settlement since AD 1850

There is a pressing need for a better understanding of changing forest fire regimes worldwide, especially to separate the relative effects of potential drivers that control burned areas. Here we present a meta-analysis of the impacts of climate fluctuation and Euro-Canadian settlement on burned areas from 1850 to 1990 in a large zone (>100 000 km2) in northern temperate and boreal forests of eastern Canada. Using Cox regression models, we tested for potential statistical relationships between historical burned areas in 12 large landscapes (reconstructed with dendrochronological data) with climate reconstructions, changes in the Euro-Canadian population, and active suppression (all reconstructed at the decadal scale). Our results revealed a dominant impact of climate fluctuations on forest burned areas, with the driest decades showing fire hazards between 5 to 15 times higher than the average decades. Comparatively, the Euro-Canadian settlement had a much weaker effect, having increased burned areas significantly only during less fire-prone climate conditions. During periods of fire-prone climate, burned areas were maximum independent of fluctuations in Euro-Canadian populations. Moreover, the development of active fire suppression did not appear to reduce burned areas. These results suggest that a potential increase in climate moisture deficit and drought may trigger unprecedented burned areas and extreme fire events no matter the effects of anthropogenic ignition or suppression.


Introduction
Burned areas and extreme fire events are increasing in many forested biomes of the world, including, the tropical forests of Amazonia (Barlow et al 2020), the temperate forests of Western North America (Holden et al 2018) and Australia (Boer et al 2020), and the northern boreal forests of Canada and Russia (Hanes et al 2019, Feurdean et al 2020. Fire significantly impacts the global earth system (e.g. climate regulation; Bowman et al 2009), ecosystem biodiversity, dynamics and services (Kelly et al 2020), as well as human lives and infrastructures (Bowman et al 2017). There is thus a pressing need for a better understanding of changing fire regimes worldwide.
Changes in fire regimes are generally assumed to be driven by interactions between climate and direct human control (Marlon et al 2008, Pechony and Shindell 2010, Kelley et al 2019, Rogers et al 2020. Climate mediates fire activity through temperature and moisture changes , Pechony and Shindell 2010, Ali et al 2012, Jolly et al 2015. Warmer temperatures lengthens the fire season and increase the probability of periods of Figure 1. Location of the 12 landscapes for which burned areas were reconstructed in eastern Canada (information about individual landscapes is provided with their abbreviation in table 1). The five landscapes that have been settled by Euro-Canadian populations are shown with red contours (see figure 2). Forest zones were modified from Baldwin et al (2018) and the FAO global ecological zones mapping.  Drever et al (2006) extreme heat, both increasing the potential burned area (Coogan et al 2021). More intense or frequent drought increases vegetation flammability, and weather-driven fuel moisture is generally considered a major limiting factor of burned areas (Holden et al 2018, Kelley et al 2019, Gaboriau et al 2020. Besides, direct human influences on fire regimes can be positive and negative (Bowman et al 2011, Knorr et al 2016. On the one hand, humans act as an essential agent of ignition, both intentionally (e.g. slashand-burn agriculture) and accidentally, significantly affecting global burned areas (Pechony andShindell 2010, Kelley et al 2019). But on the other hand, humans play a significant role in fire suppression, mostly since the second half of the 20th century with land-use driven decrease in fuel abundance and continuity, as well as modern wildfire-fighting methods (Knorr et al 2014, Andela et al 2017. A better knowledge of how those mechanisms interacted to drive historical changes in burned areas would help understand and predict current and future fire regimes (Rogers et al 2020).
This study presents a meta-analysis of the effects of climate fluctuations and Euro-Canadian settlement on historical burned areas from 1850 to 1990 in a zone of more than 100 000 km 2 in northern temperate and boreal forests of eastern Canada (figure 1). In this study area, the Little Ice Age (approximatively 1250 AD-1850 AD; Gennaretti et al 2014, Delwaide et al 2021 represented a cooler and dryer fire-prone climatic period, and burn rates experienced a major decrease around 1850 when entering the modern warmer and moister period (Drobyshev et al 2017). Notwithstanding such climate-induced long-term shifts in fire regimes, large variations in burned areas were recorded since 1850 that seem to be strongly linked to decadal-scale climate fluctuations (Bergeron et al 2004, Drever et al 2009, Chaste et al 2018.
The study area corresponds to the ancestral lands of formerly semi-nomadic Indigenous populations that have been-and still are-engaging in cultural activities in the boreal forest for at least 5000 years (Wright 1979, Côté 1996. Over the last 150 years,  Girardin et al (2006) and Girardin and Mudelsee (2008) for the eastern and western parts of the study area. Thick lines show annual raw data (Y-axis left) and vertical bars show the decadal scale anomaly (Y-axis right). The grey background shows the transition from tree-ring data to historical meteorological data. (C) Decadal increase in population (vertical bars, Y-axis left) for the five landscapes that Euro-Canadian populations have settled since 1850 (the five juxtaposed bars for each decade follow the following order of landscapes: AWS, GAS, NOS, TEN, TES). Suppression index values are shown through background blue ribbons (Y-axis right). See the section 2 for more details about these data. the study area experienced massive Euro-Canadian settlement with recorded increases of >10 000 Euro-Canadian persons per decade during the 1st half of the 20th century in some areas (figure 2). Given the generalized use of fire by Euro-Canadian settlers for agricultural clearing during these early periods of colonization (Boucher et al 2014, Terrail et al 2020, the increase in population could have had an impact similar to that of climate upon burned areas (Lefort et al 2003, Bergeron et al 2004. Later on, during the 2nd half of the 20th century, active fire suppression systems were implemented to protect Euro-Canadian communities, infrastructures and forest stands targeted for logging (Martell 2001, Blanchet 2003, Scott et al 2013, Cardil et al 2019. The relative and potentially interactive effects of climate fluctuations and direct human control (i.e. ignition and suppression) on burned areas remain poorly understood. We reanalyzed data from previous studies that reconstructed the historical burned areas per decade over the last few centuries in 12 large landscapes (figure 1) using dendrochronology. These studies (table 1) have documented regional fire regimes, as well as their climatic and environmental determinants. They have shown that burn rates are higher in conifer-dominated boreal forests than in mixed temperate forests, and higher in the more continental western part of our study area than in the eastern part . While the effects of anthropogenic actions are difficult, if not impossible, to separate from other determinants when analyzing landscapes individually (Drobyshev et al 2017), pooling the available data through a metaanalysis and using better adapted statistical methods (Cyr et al 2016) will increase statistical power and is more likely to reveal the influence of those potentially important determinants of past and future fire regimes. Through this approach, our main objective was to assess the relative and potentially interactive effects of (a) climate-driven fuel moisture fluctuations, (b) increased Euro-Canadian population, and (c) active fire suppression following the Little Ice Age (i.e. after 1850).

Material and methods
We based our meta-analysis on previous studies that reconstructed historical burned areas in 12 large landscapes located in eastern Canada (table 1). These landscapes vary in size from ≈1800 to ≈16 000 km 2 , representing a total of ≈100 000 km 2 , and are disseminated in a vast region (>500 000 km 2 ) at the transition between northern mixed temperate forests and coniferous-dominated closed canopy boreal forests (figure 1). From south to north, the study area spans four bioclimatic domains according to Québec's classification system (table 1). Mean annual temperature ranges from 4 • C to −3 • C from south to north, and mean total annual precipitations ranges from 800 mm to 1300 mm along the west-east gradient (1970-2000 climate normals; Fick and Hijmans 2017).

Burned area data
For the 12 landscapes, historical burned areas were estimated using stand initiation maps constructed from aerial photographs, historical archives, and dendrochronological data. The analysis of stand initiation maps and representative samples of the forest age structure at the landscape level, which are assumed to be the result of periodic standreplacing disturbances, provides a quantitative proxy of long-term historical disturbance rates (Wagner 1978, Johnson and Gutsell 1994, Cyr et al 2016. This method is particularly well suited to the northern forests considered in this study, where crown fires are the primary natural standreplacing disturbance, generally leaving few surviving trees with fire scars compared to forests characterized by surface fire regimes (van Horne and Fulé 2006). Stand initiation patterns are measured with dendrochronological dating of initial post-fire tree cohort establishment dates, considered as time since the last stand-replacing fire. Time since fire area per decade was then quantified for the last two or three centuries, depending on the landscapes (table 1). These data were aggregated at the decadal scale since dendrochronological methods only approximate initial postfire tree cohort establishment dates. The downscaling of these data (e.g. one or 5 years scale) and its associated lack of accuracy would seriously undermine our ability to detect correlations between historical burned areas and their potential predictors. In some cases, no trace of past fire events could be detected or precisely dated, and thus it was only possible to estimate a minimum time since fire (i.e. age of the oldest trees sampled), which are referred to as censored data in the analyses.

Climate, settlement, and suppression data
Historical drought data (figure 2) were taken from previous work that reconstructed changes in Drought Code for two large sections of our study area (eastern and western; table 1) over the past 300 years using tree-ring analysis (Girardin et al 2006, Girardin andMudelsee 2008). As part of the Canadian Fire Weather Index system, the Drought Code was developed to act as a daily index of water stored in the soil and moisture variations in deep organic layers of the forest floor (about 20 cm deep on average). The Drought Code has a very slow drying rate, with a timelag of 52 d. Drought Code reconstructions thus take account for both temperature and precipitation and provide meaningful information for inferring seasonal fire danger (Girardin et al 2013). For the period post 1913, we used Drought Code reconstructions based on historical records of meteorological stations. For older periods, we used Drought Code reconstructions based on a network of wellreplicated and accurately-dated tree-ring width chronologies (Girardin and Mudelsee 2008). We used the decadal Drought Code anomaly, calculated as the mean Drought Code per decade minus the mean Drought Code for the whole period analyzed, as a predictor variable. As tree ring width contains a mixture of climatic and non-climatic signals (e.g. biological growth trend, competition), the latter was removed through a detrending procedure before carrying the reconstructions (Girardin et al 2006, Girardin andMudelsee 2008). This procedure has the caveat of removing long-term trends associated with secular climate change (Cook et al 1995, Girardin et al 2006. In the context of our study, detrending of the treering data might have removed part of the long-term climatic trends associated with the transition from the fire-prone Little Ice Age (i.e. before 1850) to the modern, less fire-prone period, while retaining the information about decadal variability in the post-Little Ice Age period needed for our analysis (Girardin et al 2013).
Changes in Euro-Canadian population (figure 2) were reconstructed with the Canadian census data compiled at the subdivision level since 1831. Data from censuses before 1911 have been made available by the Centre interuniversitaire d'études québécoises (CIEQ; Project GÉORIA, version 2003) and data from 1911 to 1990 are freely available on Canada's Century Research Infrastructure website (https:// ccri.library.ualberta.ca). For each of the 12 landscapes, we computed the population density changes per decade, which were used in the analyses. We did not consider any effect of the baseline population because the Euro-Canadian population was almost null before 1850 in these remote areas. We thus assumed that an increase in population would mechanically lead to more land clearing and thus potentially more fire.
Finally, we collected data on fire suppression following Lefort et al (2003), who defined suppression as an ordinal index representing the three main phases of technological and efficiency development (figure 2). There was no active fire suppression effort before the 20th century (index = 0). In the early phase of fire suppression , efforts were limited to the use of shovels, axes and manual saws, and thus efficiency was minimal (index = 1). Active suppression became more efficient in the 1950s with motorized chainsaws, vehicles, water pumps, and with the development of a network of roads and fire observation towers (index = 2). Fire suppression reached its maximum efficiency (index = 3) in the 1970s with the introduction of water bombers. Moreover, while suppression mainly aimed to protect exploited forests, the area in which suppression was carried has changed through time with the northward expansion of forestry activities (supplementary figure 1 (available online at stacks.iop.org/ERL/ 16/114007/mmedia)). For a given decade, the suppression efficiency index was thus attributed only to the part of a landscape within the limits of the fire suppression zone, and an index of zero was attributed to landscapes outside the suppression zone.

Statistical analyses
We limited our analysis to the time since fire data reconstructed for the period between 1850 and 1990 for two reasons. First, the 1850s represent the end of the fire-prone Little Ice Age and the beginning of the modern warming period (Gennaretti et al 2014, Drobyshev et al 2017, Delwaide et al 2021. Second, before 1850, time since fire estimation is less accurate due to the difficulty to precisely date trees aged >150 years with dendrochronological methods (Bergeron et al 2004). A lack of accuracy in the time since fire estimation could seriously undermine our ability to detect some correlations between historical burned areas and their potential predictors (i.e. climate, population, suppression). Moreover, reliable population data was only available since the 2nd half of the 19th century. We also limited the analysis to 1990 because several of the datasets we used ended in this decade (table 1).
We tested for statistical effects of climate, population, and suppression variables on burned area since 1850 with Cox regression models (Cox 1972). This type of survival analysis is well suited as it was developed to examine time-to-event data (e.g. time since fire) and its relationship with one or more covariates, whilst solving two main limitations of our data (Bélisle et al 2011, Cyr et al 2016, Portier et al 2016. First, in Cox regressions no assumptions about the shape of the baseline hazard function is necessary. For example, that allows to implicitly account for long-term variations caused by broad bioclimatic temporal variations, which can be hard or impossible to measure. That also allowed us to focus on measurable effects of predictors, assumed to have been constant over time and additive at a given scale (Cox 1972). Second, Cox regression models (and other types of survival analyses) allow for the inclusion of censored observations (i.e. minimum time since fire estimation as described in the fire data section above), which add some incomplete but meaningful information that cannot be omitted without introducing considerable biases (Moritz et al 2009, Cyr et al 2016. Cox regressions are semi-parametric regressions consisting of two distinct parts. The 1st part is the non-parametric baseline hazard function mentioned above, which corresponds to a non specified null model for hypothetical average conditions. In other words, this baseline hazard curve corresponds to the probability of the landscape to have been affected by fire at each time step (i.e. estimated proportion of the landscape burned per decade). The 2nd part of the model is a parametric estimate of relationships between the tested predictor covariates (i.e. independent variables) and the burning hazard curve (i.e. dependent variable). Regressions were stratified by bioclimatic domains (spruce-moss, balsam fir-white birch, and balsam fir-maple-yellow birch; table 1), that are known to have different fire regimes , Gauthier et al 2015. This stratification procedure allows the Cox model to fit a separate baseline hazard curve (i.e. estimated proportion of the landscapes burned per decade) for each bioclimatic domain. Individual time-since-fire observations were weighted by their estimated burned area in km 2 . Models were fitted with the coxph function of the survival R package (Therneau 2020). We tested ten candidate models, including different combinations of potential predictors and interactions (table 2). We hypothesized that Drought Code anomaly and change in population density would positively affect burned areas, while active suppression would have a negative effect on burned areas. We then selected the best model conform with those predictions, and using the Akaike Information Criterion (AIC). We finally used the simPH package (Gandrud 2015) to compute Cox model post-estimation simulations (i.e. coxsimLinear function with 1000 simulations) and their associated hazard ratio averages and confidence intervals. The full aggregated dataset and R code used to fit Cox regressions models can be accessed at https: //doi.org/10.6084/m9.figshare.14417471.v2.

Results
The AIC values pointed to the most complex model, including Drought Code anomaly, population  (table 2). However, in all the models that included fire suppression, this variable showed an extremely high positive effect on fire risk, inverse to predictions (supplementary table 1). We considered this unexpected positive effect as aberrant considering previous studies on fire regimes in this study region (see section 4) and therefore we focussed our analysis and interpretation on the model which included Drought Code anomaly, population density, and their interaction as the best model with the lowest AIC (table 2). This model produced a supportive goodness-of-fit with a concordance of 0.68 (table 2), which means that for approximatively 68% of all pairs of observations, the modeled pairs of predictions go in the same direction (i.e. higher or lower burned areas) as those observed (Therneau and Atkinson 2020). The Drought Code anomaly, population density, and their interaction all showed significant effects on fire hazard in the selected model ( figure 3). Climate fluctuations had by far the strongest effect (figure 3(A)), with the driest decades (i.e. Drought Code anomaly = 20) showing predicted fire hazard 5-15 times larger (50% confidence intervals of model simulations) than the average decades (i.e. Drought Code anomaly = 0). Conversely, during the moistest decades (i.e. Drought Code anomaly = −10), predicted fire hazards were more than halved compared to average decades. In comparison, the increase in population density had a lower impact on fire hazard ( figure 3(B)). Predicted fire hazards were two to three times larger in decades with the highest increase in population density (i.e. 50% confidence intervals of model simulations with ≈0.7 person km 2 decade −1 ). Moreover, the significant interaction between Drought Code anomaly and increase in population density revealed that once climate exceeded a certain dryness threshold (i.e. Drought Code anomaly ≈8), the population effect became non-significant while burned areas were at their maximum ( figure 3(C)).

Discussion
Our results highlight the overwhelming effect of climate fluctuations on burned forest areas since 1850 in northeastern North America. In our study area, moisture and temperature changes have largely controlled fire activity for the last 10 000 years (Ali et al 2012, Blarquez et al 2015, Remy et al 2017. More recently, the transition from the climatically fire-prone Little Ice Age to the recent warming period (i.e. post-1850) resulted in a decrease in overall burn rates , Drobyshev et al 2017. Our results show that while the study period (1850-1990) was generally not a fire-prone one in our region, the driest decades still engendered massive fire events, culminating in the 1920s. Comparatively, the Euro-Canadian settlement had a much weaker impact on burned areas, even though fire was widely used for agricultural clearing (Boucher et al 2014, Terrail et al 2020. Indeed, our results reveal that the increase in the Euro-Canadian population was significantly associated with an increase in burned areas only during less fire-prone climatic conditions. The most likely interpretation is that settlers used fire cautiously, when weather was not conducive to large fires, which limited their impacts on burned areas. The use of fire for land-clearing by Euro-Canadian settlers was strictly ruled by Quebec's Ministry of Lands and Forests since the beginning of the 20th century, among other things, by delivering a daily burning permit during the fire season (Blanchet 2003). However, when the climate was strongly conductive to fire, burned areas reached maximum values regardless of the increase in Euro-Canadians population density. In boreal forests where lightning is a common cause of fire ignition (Woolford et al 2014, Veraverbeke et al 2017, these results suggest that lightning occurrence is more than enough to trigger very large burned areas during climatic fire-prone periods, independent of the amount of anthropogenic ignitions. This is well supported by an analysis of recent data  in our study area, which shows that despite more frequent , and their interaction on fire hazard ratio (C). The extent of simulations was limited to the 5th and 95th percentiles of each potential predictor's observed values. Dotted lines serve as reference values (i.e. hazard ratio = 1), solid lines show the central tendencies of simulations, and ribbons show the 50% confidence interval (i.e. 25th and 50th percentiles). For the interaction effect in (C), hazard ratios were simulated using a range of values of Drought Code anomaly (5th to 95th percentiles) and two extreme values of increase in population density: 0 (green) and 0.7 (purple), which represent the 5th and 95th percentiles of the observed distribution, respectively. anthropogenic than lightning ignitions, the largest burned areas are associated with lightning ( figure S2).
Some studies have also found the increase in population density to indirectly reduce burned areas through a land-use driven decrease in fuel abundance and continuity (Knorr et al 2014, Andela et al 2017. However, such a phenomenon can be considered negligible in our study area since most settled landscapes have a rather low population density and were predominantly forested through the study period. Our results strongly support our interpretation of no significant negative effects of population density upon burned areas, even when including the suppression in the models (table S1).
We did not retain models with active fire suppression because its effect on burned areas was statistically significant but with a powerful positive effect, which was inverse to that predicted. Active fire suppression can lead to increased fuel abundance (i.e. taller and denser forests) or flammability (e.g. increase in fire-prone late-successional conifers), in turn increasing fire activity in the long-term (e.g. Parisien et al 2020). However, it seems unlikely that such longterm processes have resulted in an increase in burned areas in our study area insofar as effective fire suppression only began in the last decades we took into account in our analyzes ; see section 2). The unexpected positive effect of suppression could be largely linked to the fact that our semi-quantitative suppression variable did not adequately capture the short-term fluctuations of suppression efforts (e.g. Stocks and Martell 2016). The positive effect of active suppression on burned areas may also be linked to its positive correlation with population density (Spearman ρ = 0.23). Moreover, when suppression is included in the models, it removes the significant effect of increased Euro-Canadian population on burned areas (table S1). In our study area, the Euro-Canadian settlement has been identified by several studies as an important driver that has potentially increased burned areas (Lefort et al 2003, Bergeron et al 2004, Boucher et al 2014, Terrail et al 2020, as opposed to suppression or human-driven decreased fuel abundance and continuity. A positive effect of increased Euro-Canadian population on burned areas is thus far more credible than any other direct or indirect effect of suppression. Nevertheless, our results strongly suggest that the progress of active suppression, which culminated in the 1970s with the use of water bombers, either did not reduce burned areas or merely offset the effects of other covariates. That interpretation is supported by the fact that the effectiveness of water bombers-the only available tool to control severe fires-is highly limited during extreme fire weather, especially in remote areas where prompt initial attack is operationally nearly impossible (Cumming 2005, Coogan et al 2019.
Our analysis did not take into account the effect of Indigenous populations on pre-settlement fire regimes. In northeastern North-America, fire use by Indigenous people for diverse land-management purposes has been well documented, both in northern temperate forests (Williams 2003, Munoz andGajewski 2010) and in boreal forests (Lewis and Ferguson 1988, Miller et al 2010, Johnson and Miyanishi 2012. In southern Québec, the significant increase in Indigenous population after 1500 AD has been found to shift local fire regimes from frequent low-intensity fires to less frequent but more severe crown fires (Blarquez et al 2018). However, it is very likely that Indigenous people in our study area used fire cautiously and when the weather was not conducive to fire spread (Miller et al 2010). Recent studies suggest consequently that the impacts of Indigenous land-use practices (including fire use) were mostly at the local scale (1-10 km 2 ) rather than at the regional or continental scale (10 2 to 10 3 km 2 ; Munoz et al 2014, Oswald et al 2020. Hence, large burned areas in pre-settlement times were likely triggered by climateinduced moisture fluctuations and lighting, irrespective of Indigenous ignitions, similar to the situation following Euro-Canadian settlement.

Conclusion
Our analysis shows that climatic drivers mostly control burned areas in northern temperate and boreal forests of eastern Canada, triggering extreme fire events regardless of anthropogenic ignitions. This echoes recent studies concluding that fuel moisture may be the most critical driver of burned areas and extreme fire events at the global scale (Bowman et al 2017, Kelley et al 2019. Such dominance of climate drivers has been shown to be particularly evident in temperate and boreal forests of western North America (Holden et al 2018, Gaboriau et al 2020. Our results suggest that the overall projected increase in climate moisture deficit and drought (Cook et al 2014(Cook et al , 2018, which should be associated with more frequent and intense short-term fire-prone climatic episodes (Price et al 2013), may trigger unprecedented burned areas and extreme fire events regardless of anthropogenic ignition or suppression (Podur and Wotton 2010). However, this should not be taken as evidence that fire mitigation efforts through control of land-use activities and active suppression are worthless. For example, McGee et al (2015) suggest that effective fire suppression in Canada explains why only two communities experienced major losses because of fire since 1938. Thus, fire control efforts should be concentrated on protecting communities and important infrastructure. In addition, land management planning such as forestry or infrastructure development should integrate an inevitable fire risk during increasing fire-prone climatic periods in the future.

Data availability statement
The data that support the findings of this study are openly available at the following URL/DOI: https:// doi.org/10.6084/m9.figshare.14417471.v3.

Acknowledgments
We thank Daniel Lesieur, Annie Claude Bélisle, Héloise Le Goff, Victor Kafka, Patrick Lefort, Eve Lauzon, Daniel J Grenier and C Ronnie Drever for kindly sharing their fire reconstruction data. We also thank Byron Moldofsky, Laurent Richard, and Marc St-Hilaire for their help with historical population data. The first author V D was financially supported through a MITACS-sponsored industrial partnership with Rayonier Advanced Materials forest products company (La Sarre, QC).