Will climate change drive 21st century burn rates in Canadian boreal forest outside of its natural variability: collating global climate model experiments with sedimentary charcoal data

Natural ecosystems have evolved within ranges of conditions that can serve as references for setting conservation targets or assessing the current ecological integrity of managed ecosystems (Landres et al. 1999). Disturbance regimes are key processes in many types of ecosystems and they contribute to a large extent to the creation of a variety of ecological conditions that exist through both space and time (Reynolds 2002). Disturbance regime characteristics, including frequency, spatial extent and severity, are particularly important in generating this natural variability at various spatial and temporal scales.

In this paper, we focus on the influence of disturbance frequency, which can be expressed as the mean fire interval (MFI) or its opposite, that is, the percentage annual burn rate (1/MFI; % year^–1). When stand‐replacing fires are predominant, which is the case in the continental regions of the boreal forest, the MFI is largely responsible for the creation of a complex landscape mosaic consisting of stands varying in age, composition and structure, within which other disturbances and processes interact (Wein and MacLean 1983; Johnson 1992; Payette 1992; Niklasson and Granström 2000).

Reconstruction of fire history in western Quebec boreal forests for the last 7000 years is well documented (Bergeron et al. 2004; Cyr et al. 2009). The high variability observed during the Holocene, despite slight changes in vegetation composition (Carcaillet et al. 2001, 2010; Ali et al. 2008), highlights the fact that the natural landscape mosaic has natural resilience to change in fire frequency.

Because of their climate determinism, forest fires are likely to have effects that could exacerbate biophysical and socioeconomical vulnerabilities in the context of climate change (Le Goff et al. 2005; Flannigan et al. 2009). Climate change may be defined as a change in the state of the climate that can be identified by changes in the mean or variability of its properties and that persists for an extended period of time, typically decades or longer. However, as in the United Nations Framework Convention on Climate Change, we herein more specifically define climate change as ‘a change of climate which is attributed directly or indirectly to human activity that alters the composition of the global atmosphere and which is, in addition to natural climate variability (e.g. solar, orbital forcing), observed over comparable time periods’ (United Nations 1992). Already, climate change has had a significant influence on area burned in Canada. Forest fire activity has increased steadily over the second half of the 20th century (Podur et al. 2002; Stocks et al. 2003), and part of this rise has been attributed to climate change (Gillett et al. 2004). Confounding influences may also have added up, notably from changes in the atmospheric circulation patterns governing regional fire activity (Macias Fauria and Johnson 2006; Skinner et al. 2006) or changes in land use (Podur et al. 2002).

Climate change has direct and indirect effects on forest ecosystems. Direct effects include the alteration of species growth, reproduction and migration whereas indirect effects correspond to modifications of disturbance regimes such as forest fires, insect outbreaks and diseases (Dale et al. 2001). Indirect effects, such as changes in fire regimes, may have more dramatic results than direct effects, which take a long time to materialise (Overpeck et al. 1990; Weber and Flannigan 1997).

The objective of this paper was to evaluate future trends and rates of change in fire activity under climate change in eastern Canadian boreal forests and determine whether these changes were included in the observed variability during the Holocene. First, we used regression analyses to model the historical (1959–99) relationship between climate and area burned by large forest fires (size >200 ha). These models were then used as transfer functions for predictions of future area burned. This was done by substituting historical climate conditions by future ones simulated from an ensemble of seven global climate model (GCM) experiments driven by four scenarios of radiative forcing. These predictions were used to generate changes in burn rate between three reference periods: 1959–99, 2046–65 and 2081–2100. Next, we compared these predicted changes with Holocene MFI as reconstructed by sedimentary charcoal analyses from three lakes (Cyr et al. 2009). Finally, we discuss implications of predicted future fire activity in the context of the natural range of variability observed during the Holocene.

The study area (Fig. 1) is located in the boreal zone of western Quebec and eastern Ontario (Canada). This part of the Precambrian Shield is covered with glaciolacustrine clays that have been reworked by glacial surges, resulting in a relatively compact deposit composed of clay and gravel that is called the Cochrane Till (Veillette 1994). The mostly flat region (mean altitude ∼250 m above sea level) shows three major soil types: Luvisols, Gleysols and Organic Soils (Soil Classification Working Group 1998), which reflect the variable thickness of the organic layer (∼10 to >200 cm). The dominant forest types are black spruce (Picea mariana (Mill.) B.S.P.)–feathermoss and black spruce–Sphagnum, with an understorey dominated by ericaceous shrubs (Rhododendron groenlandicum, Kalmia angustifolia, Vaccinium spp). Jack pine (Pinus banksiana Lamb.) and trembling aspen (Populus tremuloides Michx.) also occur in pure or mixed stands, and secondary tree species include balsam fir (Abies balsamea (L.) Mill.), paper birch (Betula papyrifera Marsh.), and tamarack (Larix laricina (DuRoi) K. Koch) (Gauthier et al. 2000). The region's mean annual temperature is –0.7°C (mean temperatures for January and July are –20.0° and 16.1°C respectively), and annual precipitation is 905 mm, 35% of which falls during the growing season (Matagami weather station (49°46′N, 77°49′W), Environment Canada 2006).

Fig. 1.  Location of the three lakes for which sedimentary charcoal analyses were conducted. The study area (shaded) delineates forest fire data (circles) from the Large Fire Database (LFDB, Stocks et al. 2003) that were used in the present study to develop the annual area burned model.

Forest fire data from the Large Fire Database (LFDB, Stocks et al. 2003) are used in the current study to develop the predictive model for annual area burned. These large fires (size >200 ha) represent only a very small percentage of fires but account for ∼97% of the area burned in Canada. The LFDB contains information on start location, estimated ignition date, cause and size of each fire. Fires that occurred within the territory defined as 83 to 75°W and 47 to 53°N were compiled (Fig. 1) and a time series of annual area burned covering the period 1959–99 was created. Fire is highly variable among years and the use of a large number of samples provided some statistical smoothing that improved the final calibration results and minimised the prediction error. This territory of 32.9 Mha was therefore found to offer the best compromise between selecting a region with a disturbance regime representative enough of the surroundings of the sampled lakes and maximising the signal‐to‐noise ratio in the fire dataset. Although they don't apply to the exact same periods, the recent burn rates estimated for the 32.9‐Mha territory from the LFDB during the 1959–99 period was 0.24% per year whereas it was estimated to 0.27% during the 1920–98 period for a smaller area located in the surroundings of the sampled lakes from a dendroecological study (Bergeron et al. 2004).

The predictors of area burned used in this study included monthly means of daily mean temperatures and the monthly Drought Code (hereafter called MDC), an index representing the net effect of changes in evapotranspiration and precipitation on cumulative moisture depletion in deep, compact organic layers (Girardin and Wotton 2009). The MDC was previously shown to be a robust predictor of area burned across circumboreal forests (Girardin et al. 2009a). Monthly mean temperatures also proved to be a robust predictor of area burned across North American boreal forests (Balshi et al. 2009), and therefore it is used in our model calibration.

The MDC is a generalised monthly version of the daily Drought Code widely used across Canada by forest fire management agencies in their monitoring of wildfire risk. The MDC was developed by Girardin and Wotton (2009) to be used in seasonal drought characterisation analyses when the daily weather data necessary for computation of the daily Drought Code are not available.

The Drought Code represents the moisture content of organic matter that is on average ∼18 cm thick and 25 kg m^–2 dry weight, for a bulk density of 138.9 kg m^–3. The equation linking the Drought Code (DC) to its moisture equivalent (Q; unitless) is:

The 400 constant in Eqn 1 represents the maximum theoretical moisture content of the fuel represented by the Drought Code, which roughly corresponds to the water‐holding capacity of the soil, i.e. 100 mm (Van Wagner 1987). In its daily version, the Drought Code has a response time of 62 days at 15°C and 44 days at 30°C. This long response time served as the basis for the development of a monthly approximation model.

The MDC formulation may be summarised as follows. First, potential evapotranspiration E (unitless quantity) over month m follows the method of Thornthwaite and Mather (1955) and is given by:

where T_mx is the monthly mean of daily maximum temperatures (°C), L_f is the standard day‐length adjustment factor (Van Wagner 1987), and n is the number of days in the month. Next, computation of the MDC is carried out twice in a month to reduce bias that may arise when the forest floor becomes saturated in the spring. Drying taking place over the first half of the month (DC_HALF) is calculated as:

where MDC₀ is the MDC from the end of the previous month. Total monthly rainfall (r_m) is simulated to occur in the middle of the month, and the moisture equivalent in the layer after rain Q_mr (unitless) is calculated as:

Drying taking place over the middle of the month (DC_mr) is calculated as:

In these equations, r_m is reduced to an effective rainfall (RM_EFF; mm) after canopy and surface fuel interception using RM_EFF =0.83r_m. An estimate of the MDC value at the end of the month (MDC_m; unitless), for which total rainfall and mean temperature apply, is calculated as:

Finally, MDC₀ and MDC_m are averaged to find a mean drought value for the month:

When calculating the next month's MDC, the value of MDC_m from the previous month then becomes the new MDC₀.

Girardin and Wotton (2009) found that the MDC generally followed trends in the monthly means of the daily Drought Code closely (r² ranging from 0.87 to 0.95 for n = 612 sample months). There are no absolute guidelines as to the meaning of MDC values but, generally speaking, Drought Code values below 200 are considered low whereas values around 300 may be considered moderate in most parts of the country. A Drought Code rating of 400 or more indicates that fire could involve burning of deep subsurface and heavy fuels (e.g. de Groot et al. 2009). The index generally peaks in mid‐ to late August, beyond which it either declines or maintains the same value; during extreme years with late‐season fires, this does not always hold true (McAlpine 1990; Girardin et al. 2004). The reversal in August is only attributed to a change in day length, and is not a function of seasonal precipitation.

Monthly means of daily mean and maximum temperatures and monthly precipitation totals were obtained for the 1959–99 period using BioSIM (Régnière and Bolstad 1994). Data were obtained for 100 locations distributed across 83 to 75°W and 47 to 53°N by interpolating data from the four closest weather stations to each location, and adjusting for differences in latitude, longitude and elevation between the data sources and each of the locations. The MDC and mean monthly temperatures were computed at each location and data were then averaged to create regional mean monthly data. In order to give greatest weight to climate anomalies in fire‐prone regions (Gillett et al. 2004; Girardin and Wotton 2009), the distribution of locations was made such that regions having higher fire frequencies also had higher location replicates (see Fig. 1).

Prediction of future annual area burned was made using simulated monthly temperature and precipitation data collected from an ensemble of seven GCMs (Table 1). Model selection was made accordingly with the availability of monthly means of daily maximum temperature outputs necessary for simulation of MDC. GCMs are time‐dependent numerical representations of the atmosphere and its phenomena over the entire planet, using the equations of motion and including radiation, photochemistry, and the transfer of heat, water vapour, and momentum. Future climate scenarios are built based on the effects of various concentrations of greenhouse gases and other pollutants within the atmosphere on the earth–atmosphere system. Monthly temperature and precipitation data were collected from four to six cells (depending on model resolution) located near the region, and averaged. Data were collected for horizons 1961–99, 2046–65 and 2081–2100.

Table 1.  General circulation models from the Intergovernmental Panel on Climate Change Fourth Assessment Report (AR4; Meehl et al. 2007) and their greenhouse gas (e.g. radiative) forcing scenarios

GCMs can have large biases when it comes to reproducing the regional features of climate. To account for differences between the actual climate data derived from BioSIM and the GCM predictions, we adjusted the monthly simulations relative to the absolute difference from the 1961–99 monthly means of actual data (e.g. Balshi et al. 2009). A correction was also applied to the interannual variability by changing the width of the distributions so that mean monthly GCM predictions and data derived from BioSIM had equal standard deviations over their common period 1961–99 (Girardin and Mudelsee 2008). The following algorithm was used in doing these corrections:

where GCM_adj is the corrected monthly value output of the GCM predictions, BioSIMµ is the mean monthly value (across all years) for the period 1961–99 derived from BioSIM, GCM_month is the monthly value output by GCM predictions, GCMµ is the mean monthly value (across all years) for the period 1961–99 derived from the GCM monthly predictions, and δ is the ratio of the standard deviation in monthly values (across all years) for the period 1961–99 derived from BioSIM compared with that derived from GCM monthly predictions (δ = BioSIMσ/GCMσ). This ratio inflates or deflates the simulated interannual variations.

Four scenarios of projected changes in greenhouse gas emissions (Nakicenovic et al. 2000) are used in the present study (from ‘worst’‐ to ‘best’‐case scenarios): A2, A1B, B2 and B1. Each scenario reflects a specific storyline of future development such as global population growth, economic development, and technological change. The storylines describe the relationships between the forces driving greenhouse gas and aerosol emissions and their trajectories over the 21st century (http://sedac.ciesin.org/ddc/sres/index.html, accessed 22 May 2009):

• A2 storyline (intense forcing): a very heterogeneous world with continuously increasing global population and regionally oriented economic growth that is more fragmented and slower than in other storylines.

• A1B storyline (intense forcing): a future world of very rapid economic growth, global population that peaks mid‐century and declines thereafter, and rapid introduction of new and more efficient technologies with low carbon emissions.

• B2 storyline (intermediate forcing): a world in which the emphasis is on local solutions to economic, social and environmental sustainability, with continuously increasing population (lower than A2) and intermediate economic development.

• B1 storyline (intermediate forcing): a convergent world with the same global population as in the A1B storyline but with rapid changes in economic structures towards a service and information economy, with reductions in material intensity, and the introduction of clean and resource‐efficient technologies.

These storylines are intended to cover a wide spectrum of alternative futures to reflect relevant uncertainties and knowledge gaps associated with climate change issues (Nakicenovic et al. 2000).

Development of a predictive model for area burned was carried out using Multivariate Adaptive Regression Splines (MARS) (Friedman 1991). MARS is a technique in which non‐linear relationships between a predictand (i.e. variable to predict) and a predictor are described by a series of linear segments of differing slopes, each of which is fitted using a basis function. Breaks between segments were defined by an inflection point in a model that initially overfitted the data, and was then simplified using a backward and forward stepwise cross‐validation procedure to identify terms to be retained. At each step, the model selected the inflection point and its corresponding pair of basis functions that gave the greatest decrease in the residual sum of squares. Selection was proceeded until some maximum model size was reached, after which a backward‐pruning procedure was applied in which those basis functions that contributed least to model fit were progressively removed. The sequence of models generated from this process was then evaluated using generalised cross‐validation, and the model with the best predictive fit was selected. For the present study, annual area burned was regressed against the MDC and mean temperatures from April to October over the 1959–99 period (total of 14 potential predictor variables). Once the model is developed, it may then be used as a transfer function for prediction of future area burned. This is done by substituting historical climate conditions by future ones generated from the GCM experiments. For other applications of MARS, see Leathwick et al. (2005), Balshi et al. (2009), and Girardin et al. (2009b). The R package ‘earth’, specifically designed for MARS using the techniques described by Friedman (1991), was used for calibration of model parameters (R Development Core Team 2007). The Generalised Cross Validation (GCV) penalty per knot was set to two. Other parameters were kept as in the ‘earth’ default settings.

A long‐term fire history was reconstructed using sedimentary charcoal from three lakes and dated using ¹⁴C and ²¹⁰Pb isotopes (Carcaillet et al. 2001, 2010). These three lakes were considered as representative of the study area. Only charcoal fragments larger than 150 µm were considered as particles of this size generally do not travel more than a few hundreds of metres from a fire (Higuera et al. 2007). Charcoal accumulation peaks were then isolated to build fire event chronologies beginning as far back as 7600 calendar years before present (BP) for one of the lakes, although only results spanning the last 6800 years were available for all three lakes (see method in Carcaillet et al. 2001). The fire intervals from the three lakes were pooled to estimate the MFI across the landscape, which we define as the average number of years between two successive fire events at a given point (cf. Merrill and Alexander 1987). We assumed that the ratio between the charcoal source‐area and the typical size of fires in this area is negligible and therefore, we considered these lakes as a point‐based representation of fire activity. Approximately 80% of the total area burned between 1959 and 1999 was indeed affected by individual fires larger than 8500 ha and charcoal source‐area represented by these lakes was most likely between several tens of and a few hundred hectares (Higuera et al. 2007). Li (2002) demonstrated that the derivation of point‐based estimates of fire activity such as the MFI to obtain an area‐based one is possible when the sampled sites are independent of each other and if the observations cover a period of several times the MFI. Both conditions are met as our paleoecological fire history reconstruction covers 6800 years and as the nearest lakes are located 56 km from one another. A similar paleoecological reconstruction conducted nearby from sedimentary charcoal also shows an asynchrony of recorded fire events in lakes that are actually closer to one another (Ali et al. 2009). Our estimations of MFI were then reversed to obtain estimations of average burn rates (burn rate = 1/MFI; cf. Li 2002).

To assess the long‐term natural variability in MFI (and burn rate), we first estimated MFI within periods of relatively constant regimes using the two‐parameter Weibull probability density distribution fitted over the observed fire intervals (not shown; see Carcaillet et al. 2001 or Cyr et al. 2009). The Weibull‐modelled MFI during these periods were used to define a first range of long‐term natural variability that will be referred to as the conservative range, which is suggested as a management guideline by Cyr et al. (2009), whereas the corresponding 95% confidence intervals were used to define an extended range of natural variability. The second approach was a smoothing method, which we used to fit the Weibull distribution within a moving window of 13 observations and report the estimated MFI along our time series. This moving window roughly corresponded to a little more than 1000 years during periods with low fire frequency and ∼300 years during periods with high fire frequency. The smoothing method shows shorter‐term variations, which are reported for comparative purposes.

Simulations from the seven GCMs showed agreement as to the direction of future temperature changes in our study area, regardless of radiative forcing scenarios (Fig. 2). Warming is expected to occur throughout all months of the year, with winter months showing the greatest rates of change. By the end of the 21st century, daily maximum temperatures during summer months are predicted to be warmer by ∼4°C when compared with 1961–99.

Fig. 2.  Predicted changes in monthly averages of daily mean and daily maximum temperatures, and total monthly precipitation between (a) 2046–65 and 1961–99 and (b) 2081–2100 and 1961–99. Simulations were obtained using an ensemble‐mean of seven global climate models forced by various scenarios of greenhouse gas emissions (Table 1), for a total of 19 simulations. The boxplot shows the distribution of the 19 simulation results for January (J) to December (D). The boundaries of the box denote the 25th and 75th percentiles, and the line within the box marks the median of the simulated mean monthly anomalies. Error bars above and below the box indicate the 90th and 10th percentiles.

The rate of change in precipitation predicted from the ensemble of GCMs is not as unequivocal as the one in temperature. Approximately half of the simulations predicted a decrease in summer precipitation (down by 18% when comparing 1961–99 with 2081–2100 horizons), whereas an increase in precipitation was predicted by the other half (up by 11%). An increase in precipitation is predicted for fall (autumn), winter and spring months by all models (from 8 to 41% in spring and from –1 to 20% in fall).

The influence of climate on forest fire involves a trade‐off between the amount of precipitation and temperature. Our estimates of the net effects of changes in evapotranspiration and precipitation on cumulative moisture depletion in soils suggest that the uncertainty in modelled precipitation changes will be exacerbated by the large increases in temperature. In 17 out of 19 of our simulations (i.e. 90% of our simulations), MDC was predicted to increase by the end of the 21st century, resulting in an increase median ensemble MDC in all months from April to October (Fig. 3). In the few cases where the MDC was predicted to decrease (particularly during the mid‐century), the emission scenario was forced by a B1 storyline of low population growth with rapid technological changes (intermediate forcing). It is noteworthy to mention that in spite of this general agreement with regard to the direction of summer moisture changes, uncertainty remains important, as suggested by the 10th and 90th percentiles, up to 80 MDC units for August (Fig. 3).

Fig. 3.  Predicted changes in Monthly Drought Code (MDC) between (a) 2046–65 and 1961–99 and (b) 2081–2100 and 1961–99. See Fig. 2 for definition.

We regressed the annual area burned (AAB) against the MDC and mean temperatures over 1959–99 using MARS. The regression model explained 42% of the deviation between the variable to predict and the predictors (P < 0.001; GCV R² = 0.36) and took the following form:

where MDC_August is the August MDC. The calibration error in AAB estimates was 0.061 Mha per year. In the MARS model, the inflection point for the MDC variable takes on a value of 192.36, and the basis function (BF₁) takes on a value of MDC – 192.36 when MDC is >192.36, but otherwise takes on a value of zero (Fig. 4). Coefficients applied to the basis function define the slopes of the non-zero sections. Dividing the calibration interval (see Girardin and Mudelsee 2008) led to partial changes in the relationship. Specifically, July MDC showed up as a potential AAB predictor depending on period divisions. However, the July MDC was discarded as it is highly collinear with the August MDC (Pearson r = 0.89, P < 0.001). Monthly means of daily mean temperature did not show a significant relationship with AAB.

Fig. 4.  Relationships between the August Monthly Drought Code (MDC) and annual area burned (AAB) modelled using Multivariate Adaptive Regression Splines.

Predicted changes in burn rates from the ensemble of seven GCMs forced by various scenarios of greenhouse gas emissions are shown in Fig. 5. A wide array of model and scenario outcomes suggests a fairly large amount of uncertainty in predictions of future burn rates. Simulations from the Canadian and the National Institute for Environmental Studies GCMs (i.e. CGCM3T63 and MIROC3.2 medres) yielded the greatest rates of change. For instance, the CGCM3 A2 simulation suggested a change in burn rate from a current value of 0.24% year^–1 to one of 0.69% year^–1 for the end of the 21st century. Other models, like the GISSAOM, suggested moderate changes in burn rate in the course of the 21st century. All models suggested that the future burn rate will remain below 1.00% year^–1.

Fig. 5.  Predicted changes in burn rate over the 21st century simulated from an ensemble of seven global climate models forced by various scenarios of greenhouse gas emissions (Table 1). A 95% bootstrap confidence interval for the ensemble‐mean burn rate is shown (grey‐shaded). The latter was obtained after randomisation of the annual area burned (AAB) simulations and recomputation of the burn rate for each horizon using the upper and lower AAB percentiles. The yellow shading denotes the three horizons (1961–99, 2046–65 and 2081–2100) used in the computation of the future burn rate. See Table 2 for individual ensemble‐mean scenarios. The burn rate estimated from the Large Fire Database (LFDB) is indicated by a red square.

Modelling uncertainty was dealt with by randomising of the AAB simulations using a bootstrap method and recomputing the burn rate for each horizon using the 2.5, 50 and 97.5% AAB percentiles. This yielded an ‘ensemble‐mean’ burn rate and a 95% confidence interval (CI) for each horizon (Fig. 5). Not only does this ‘uncertainty band’ take into account climate modelling and greenhouse gas emissions uncertainties, but it may also partly take into account errors owed to various assumptions in our data and the calibration error (a calibration error of 0.061 Mha per year corresponds to an error in the estimated burn rate of ±0.18% year^–1). Burn rates for ensemble‐mean scenarios A2, A1B and B1 were also computed (Table 2). However, one should note that ensemble‐mean scenarios were obtained from limited GCM experiments (six GCM experiments per scenario), which may explain the wide uncertainty bands. The ensemble mean of all 19 simulations suggests an increase in the burn rate of our study area from a current value of 0.20% year^–1 (95% CI (0.12, 0.28)) to 0.36% year^–1 (95% CI (0.25, 0.48)) by the mid‐21st century, and 0.45% year^–1 (95% CI (0.32, 0.59)) by the end of the 21st century. Upper and lower bounds of this ensemble‐mean CI (Fig. 5) closely approximate burn rates calculated from ensemble‐mean scenarios A2 and B2 respectively (Table 2). Will these values be within the historical variability?

Table 2.  Predicted burn rates (% year^–1) from ensemble‐mean scenarios for horizons 1961–99, 2046–65 and 2081–2100
Each ensemble simulation shown here consists of an average of six general circulation model experiments undertaken with identical forcing scenarios (A2, A1B and B1). A 95% bootstrap confidence interval (CI) for the ensemble‐mean burn rate is indicated in parentheses. General circulation models (GCMs) used in each climate change scenario

The period between 6800 and 3200 years BP was characterised by a relatively low burn rate, whereas the period after 3200 years BP shows a considerable increase with a subsequent decrease during the last millennium (Fig. 6). Climatic drivers are currently the most plausible explanation for these changes, as pollen records show no relationships between vegetation, based on pollen composition and diversity, and fire activity during this period (Carcaillet et al. 2001), both in terms of trends and timing of paleofire and vegetation modification (Carcaillet et al. 2010). Using these periods of relatively constant regimes, the conservative range of variability in burn rate varies between 0.37% year^–1 (MFI = ∼267 years) and 0.90% year^–1 (MFI =∼111 years), whereas the extended range varies between 0.24% (MFI = ∼419 years) and 1.22% year^–1 (MFI = ∼82 years) (Cyr et al. 2009). The smoothing method shows shorter‐term variations, which are generally encompassed by the extended range of natural variability.

Fig. 6.  Long‐term variability in burn rate compared with the burn rate for 2081–2100. Average burn rates during periods of relatively constant regime (straight lines) define the conservative range of natural variability, whereas the 95% confidence intervals (shaded area) define the extended range of natural variability. Shorter‐term variations are illustrated for comparison purposes. The extreme experiment of MIROC3.2 medres A2 is indicated by an asterisk.

The increase in burn rate that is predicted by the end of the 21st century, +0.45% year^–1 (95% CI (0.32, 0.59)), falls within the long‐term variability that was derived from the paleoecological reconstruction. This predicted burn rate, which corresponds to an MFI of 222 years, is relatively low when compared with periods of higher burn rates, i.e. between 3300 and 900 years BP, when the burn rate was oscillating at ∼0.9%. Moreover, the 95% CI associated with the predicted burn rate (corresponding to a range of MFI varying between 169 and 313 years) is also well contained and quite smaller than the conservative and extended ranges of natural variability.

This study used an ensemble‐mean of 19 GCM experiments for prediction of future burn rates in eastern Canadian boreal forests. These experiments investigate the effects of different, but equally plausible, initial atmospheric and oceanic conditions and CO₂ forcing scenarios on climate and fire behaviour (Intergovernmental Panel on Climate Change 2007). Individual GCM experiments provide indications of the magnitude of the natural variability in the system. This variability is the main cause for the large uncertainty in GCM experiments under a given CO₂ scenario (see confidence intervals in Table 2). However, the ensemble mean allows one to distinguish the approximate climate change signal from the natural variability. The ensemble mean allows the identification of the effect response of burn rate to the climate change signal. The only other study of future fire activity in boreal forests that we are aware of in which an ensemble‐mean of GCM experiments was used is that of Scholze et al. (2006).

The main result of our study is that, whatever the CO₂ emission scenarios or the GCM used, the future burn rate is predicted to remain within the natural range of variability for this region of the boreal forest. Burn rates predicted by the ensemble mean of all 19 GCM experiments (Fig. 5) and by individual ensemble‐mean scenarios (Table 2) are lower than assessed by sedimentary charcoal data during the Holocene period. Even under the most extreme experiment of MIROC3.2 medres A2, a 0.8% year^–1 burn rate (equivalent to an MFI of 125 years) was still found to be slightly below what was observed between 3300 and 900 years BP (Fig. 6; see also Cyr et al. 2009).

Our assessment is based on long‐term variability in burn rate reconstructed from sedimentary charcoal data. The plausibility of the obtained results rests on the validity of two principal assumptions. One is that forest fire distribution is relatively homogeneous in the study area and across time, such that locations of the selected lakes are representative for the whole study area. Second, MFI is estimated from peaks in charcoal sediments originating from a local source area around the lakes. However, present forest fire distribution is not homogeneous in the study area (see Fig. 1), and this for a number of factors including variability in land‐use and in fire detection and suppression (Lefort et al. 2003). Also, a charcoal source area may sometimes be undefined; there is always a probability with charcoal sediment records that fire events originating from large and distant fires will be detected as small local fires (Higuera et al. 2007). Such uncertainty may be taken into account by construction of CI (Fig. 6) that may encompass some range of sampling errors.

To increase confidence, results may also be supplemented by other reconstruction methods of past fire activity. Past fire activity in this part of the boreal forest has been well documented with time‐since‐fire distributions reconstructed from dendrochronological dating of forest stands and fire scars (Bergeron et al. 2004, 2006). Consistent with our results from the charcoal sediments record, reconstruction of burn rates for the last 300 years for a territory covering 15 000 km² (Bergeron et al. 2004) reports variations in burn rate that encompass the entire array of predicted burn rates for 2100 AD (Fig. 5; Table 2). The 0.8% year^–1 burn rate predicted for 2100 AD associated with the most extreme experiment (MIROC3.2 medres A2) is indeed lower than what was observed during the Little Ice Age (before 1850 AD), where the burn rate was the highest recorded in recent history, at ∼1.09% (equivalent to an MFI of 92 years). Girardin and Mudelsee (2008) in their analysis of late Holocene fire risk variability inferred from tree‐ring records collated to GCM modelling also came to a similar conclusion. Their results suggest probabilities of having future levels of extreme fire risk conditions within the historical range.

The conclusions drawn from this work should also be limited to the effect of late spring and summer climate variability on forest fire activity. Factors that are not directly taken into account in the seasonal drought severity component could modulate the projected increase in burn rates over the 21st century and could distort the predicted trend. Effects of forest composition (coniferous v. hardwood) and age structure on fuel availability and moisture regimes, which are important determinants of fire activity under a given climate (Hély et al. 2001), are not taken into account in the current simulations. These limitations also apply to the work done by Flannigan et al. (2005), Bergeron et al. (2006), Girardin and Mudelsee (2008), and Balshi et al. (2009). Changes in fire regimes could lead to shifts in vegetation composition and structure that could provide feedback to fire activity. However, such a feedback doesn't appear to be important in the region as no significant or delayed effect on the pollen‐inferred vegetation was detected in relation to the variation in burn rate that was observed during the last 7000 years (Carcaillet et al. 2010; see section below for additional details). Furthermore the dominance of black spruce and a thick organic layer could very likely limit rapid changes in vegetation composition (Lecomte et al. 2006). Other factors not considered by the predictive model that may be of importance include changes in the frequency of small precipitation events and their effects on fine fuels, and changes in wind velocity and their effects on fire behaviour (Li et al. 2000). That being said, these variables are currently not important predictors of area burned in the study area (Balshi et al. 2009; their grid points 80.0°W, 50.0°N and 80.0°W, 52.5°N presented in their table A1). Other factors that were not taken into account in the present study include changes in ignition agents (lightning frequency and human‐caused ignition; Price and Rind 1994; Wotton et al. 2003), changes in land use (e.g. fragmentation of landscapes), interactions with other natural disturbance agents such as insect outbreaks and diseases, feedback to the climate system through increases in trace gas emissions (Gillett et al. 2004), and alteration of surface energy exchanges (Chambers and Chapin 2003). Simulations of future fire conditions using dynamic climate–vegetation models (de Groot et al. 2003; Keane et al. 2004) would be relevant when attempting to account for changes in lightning probability, vegetation and fuel types under a changing climate regime.

Despite major changes in burn rate, pollen‐inferred vegetation reconstructions were unable to report clear concurrent shifts in the regional vegetation linked to fire during the last seven millennia (Carcaillet et al. 2001, 2010). These studies suggest that vegetation is rather well adapted to fluctuations in fire activity, although some species might locally have behaved individually (Ali et al. 2008). The stability of available pollen records is particularly strong in black spruce‐dominated areas (see also eastward, Garralla and Gajewski 1992), a species that can cope well with the range of fire‐free intervals occurring under natural conditions. Its seeds stored in serotinous cones allow for quick regeneration following fire, whereas its shade tolerance and layering capacity allow for its persistence in the landscape without fire (Lecomte and Bergeron 2005). Although less abundant, both early successional species, such as jack pine and aspen, and late successional species, such as balsam fir or white spruce, are currently present in the forest mosaic (Gauthier et al. 2000). This emphasises the possibility that the burn rate was never sufficiently high or low to totally exclude them from the landscape, although locally significant effects are observed on balsam fir (Abies balsamea) close to its northern range limit (Ali et al. 2008). Minor changes in species’ relative abundances or in the forest age structure could therefore occur with expected changes in burn rates, but available pollen records suggest that these changes are unlikely to cause a major shift towards an alternative state.

Moreover, the considerable variations in burn rate that were observed during the last 300 years (Bergeron et al. 2004; Girardin and Mudelsee 2008) show that the rate of change can be quite high under natural conditions. This observation suggests that systems should be able to cope with the predicted change in the burn rate of +0.2% during the next century (ensemble mean of all 19 GCM experiments, Fig. 5).

Fire activity is already an important issue for public safety and timber protection (Martell 1994). The expected doubling of the burn rate calculated for 2100 will have very significant economical and social consequences that will require important adaptive measures. Adaptive measures suggested previously include increased fire suppression efforts, fuel management, salvage logging, regeneration enhancement, functional zoning and risk assessments that are considered a priori in general strategic planning and annual allowable cut (AAC) calculations (Le Goff et al. 2005; and references therein). Independently of the implementation and success of all these measures, the maintenance of ecosystem resilience is a prerequisite to the sustainable management of the boreal forest, as it prevents the system from shifting towards an undesirable alternative state and thus reduces uncertainty.

Forest ecosystem management, an approach that is gaining in popularity and that aims to decrease the difference between managed and natural forests (Gauthier et al. 2009) by emulating natural disturbances through harvesting, may provide additional solutions (see also Attiwill 1994; Angelstam 1998; Harvey et al. 2002; Kuuluvainen 2002). The rationale for such an approach is that fire is perhaps the most important determinant of landscape compositions and stand structures and that maintaining these through forest ecosystem management practices should in turn allow us to maintain biological diversity, essential ecological functions and ecosystem resilience. As burn rates get closer to the upper boundaries of the natural range of variability, it decreases the necessary ‘room for manoeuvre’ that allows the substitution of fire for harvesting as a high‐severity disturbance, hence limiting the extent to which we can emulate fire (Bergeron et al. 2006). To maintain a comparable AAC while limiting the cumulative effects of these two stand‐initiating disturbances (fire and clear‐cutting), which together are driving the system out of its natural range of variability (Cyr et al. 2009), we suggest increasing the relative importance of uneven‐aged management, which can be designed to emulate finer‐scale disturbances or individual mortality within stands. First, uneven‐aged management would contribute to maintaining old forest attributes in a landscape where they were naturally and historically dominant. Second, it would contribute to preventing the additional costs of regeneration failures that occur more frequently when high‐severity disturbances such as clear‐cutting, fire and insect outbreaks happen within short intervals of time (Jasinski and Payette 2005), a phenomenon that is likely to increase as the pressure from low‐retention harvesting and fire does.

Our results support precedent studies reporting that climate change will cause an increase in burn rate over this part of the boreal forest (Flannigan et al. 2005; Girardin and Mudelsee 2008). The increase in burn rate that we predict, from the present value of 0.20 to 0.45% year^–1 with 95% (0.32, 0.59%) by the end of the 21st century, might appear important at first but it is relatively modest in comparison with the natural range of variability. Although our results suggest that the predicted increases in burn rate per se will not move this ecosystem to new conditions never before encountered in the past, the cumulative effects of fire and clear‐cutting or other low‐retention types of harvesting, which still prevail in this region, remain preoccupying. It has already been shown that clear‐cutting has considerably altered this system in terms of age‐class representation at the landscape level by diminishing the amount of stands exceeding in age the length of the typical harvest rotation (Bergeron et al. 2006; Cyr et al. 2009). An excessive use of even‐aged management, therefore, contributes to eroding the ecological resilience by reducing ecosystem variability in time and space (Drever et al. 2006), a process that will be exacerbated by the predicted increase in burn rate.