The response of wildfire regimes to Last Glacial Maximum carbon dioxide and climate

Climate and fuel availability jointly control the incidence of wildfires. The effects of atmospheric CO₂ on plant growth influence fuel availability independently of climate, but the relative importance of each in driving large-scale changes in wildfire regimes cannot easily be quantified from observations alone. Here, we use previously developed empirical models to simulate the global spatial pattern of burnt area, fire size, and fire intensity for modern and Last Glacial Maximum (LGM; ∼ 21 000 ka) conditions using both realistic changes in climate and CO₂ and sensitivity experiments to separate their effects. Three different LGM scenarios are used to represent the range of modelled LGM climates. We show large, modelled reductions in burnt area at the LGM compared to the recent period, consistent with the sedimentary charcoal record. This reduction was predominantly driven by the effect of low CO₂ on vegetation productivity. The amplitude of the reduction under low-CO₂ conditions was similar regardless of the LGM climate scenario and was not observed in any LGM scenario when only climate effects were considered, with one LGM climate scenario showing increased burning under these conditions. Fire intensity showed a similar sensitivity to CO₂ across different climates but was also sensitive to changes in vapour pressure deficit (VPD). Modelled fire size was reduced under LGM CO₂ in many regions but increased under LGM climates because of changes in wind strength, dry days (DDs), and diurnal temperature range (DTR). This increase was offset under the coldest LGM climate in the northern latitudes because of a large reduction in VPD. These results emphasize the fact that the relative magnitudes of changes in different climate variables influence the wildfire regime and that different aspects of climate change can have opposing effects. The importance of CO₂ effects imply that future projections of wildfire must take rising CO₂ into account.

1 Introduction

Climate influences the occurrence of wildfires both through fire weather, which affects the probability of wildfire start and spread, and the long-term establishment of vegetation which is strongly controlled by temperature and precipitation (Bradstock, 2010; Pausas and Ribeiro, 2013). It has been suggested that current climate change, driven by increasing atmospheric CO₂ levels, will increase wildfire risk in many regions through increased fuel dryness whilst potentially reducing wildfire risk in some regions due to decreasing fuel availability (e.g. Abatzoglou et al., 2019; Bowman et al., 2020; Harrison et al., 2021; Rogers et al., 2020). However, atmospheric CO₂ levels also affect fuel loads independently of climate through physiological effects on photosynthesis, which cascade into plant growth rates (Bond et al., 2003; Bond and Midgley, 2012; Kgope et al., 2010). Much emphasis has been placed on recent and future changes in fire weather (see e.g. Abatzoglou et al., 2019; Betts et al., 2015; Flannigan et al., 2013; Jolly et al., 2015). However, increases in atmospheric CO₂ concentrations promote vegetation productivity, thus altering fuel availability and loads, as well as affecting fuel types through e.g. woody thickening (Buitenwerf et al., 2012; Donohue et al., 2013; Knorr et al., 2016; Martin Calvo et al., 2014; Martin Calvo and Prentice, 2015; Pausas, 2014). Fuel properties have different effects on different aspects of the fire regime, with fire size being strongly constrained by fuel continuity and with fire intensity being limited by fuel loads (Archibald et al., 2013; Haas et al., 2022). Thus, CO₂-induced changes in vegetation properties will most likely affect these aspects of wildfire regimes differently.

One reason the impact of CO₂ on wildfires is poorly constrained is the difficulty of isolating it based on observations alone. Satellite records only span ∼ 25 years, a relatively short period to monitor the effect of changing CO₂ levels on the vegetation properties that influence wildfires. Furthermore, changes in atmospheric CO₂ levels and climate are temporally correlated, and since both affect vegetation, it is difficult to attribute changes in observations to one or the other. An alternative approach is to use process-based fire-enabled vegetation models which explicitly account for the physiological effects of CO₂ and can be used to examine the temporal and spatial patterns of wildfires under different conditions. Process-based models have been used to examine the impact of climate and atmospheric CO₂ changes on both vegetation and wildfire at the Last Glacial Maximum (LGM; 21 000 years ago) (Martin Calvo et al., 2014; Martin Calvo and Prentice, 2015). The LGM is a useful out-of-sample experiment since the climate forcing is of similar magnitude to the change expected by the end of the century in high-end scenarios, though it is of opposite sign (Kageyama et al., 2021). The LGM had a generally colder and drier climate than today, with CO₂ levels of ∼ 185 ppm. Palaeorecords show reduced vegetation productivity and forest cover (Harrison and Prentice, 2003; Kaplan et al., 2016; Moreno et al., 2018), and ice core and sedimentary charcoal records indicate reduced biomass burning globally (Albani et al., 2018; Harrison et al., 2022; Marlon et al., 2016; Rubino et al., 2015). Although this reduction could reflect the colder and drier conditions, model experiments suggests that low CO₂ also played a crucial role. Experiments using the coupled biogeography and biogeochemistry model BIOME4 (Kaplan et al., 2003) showed that it was necessary to include the direct effect of CO₂ to simulate observed global and regional reductions in forest cover during the Glacial (Bragg et al., 2013; Harrison and Prentice, 2003). Similarly, Martin Calvo et al. (2014) showed that low CO₂ was necessary to simulate the observed reduction of biomass burning in LGM experiments using the LPX (Land Surface Processes and eXchanges) fire-enabled vegetation model.

In this analysis, we use three empirical models (Haas et al., 2022) to explore the relative importance of climate and of CO₂ to the global spatial patterns of burnt area, fire size, and fire intensity. We performed two experiments under realistic modern CO₂ and climate conditions (MOD climate/MOD CO₂ and LGM climate/LGM CO₂). We also performed two counterfactual sensitivity experiments to quantify the sensitivity of each wildfire property to climate and CO₂ independently (MOD climate/LGM CO₂ and LGM climate/MOD CO₂). Comparisons to LGM charcoal records from the Reading Palaeofire Database (RPD) (Harrison et al., 2022) were used to examine which experiments provided the most realistic spatial patterns.

2 Methods

Haas et al. (2022) developed empirical models of the global spatial patterns of burnt area (BA), fire size (FS), and fire intensity (FI) using generalized linear modelling (GLM) of modern observations. Here we use these models to simulate the global spatial patterns of burnt area (BA), fire size (FS), and fire intensity (FI) under four climate and CO₂ scenarios (Fig. 1). We used two realistic scenarios: (a) MOD climate and CO₂ conditions and (b) LGM climate and CO₂. We ran two sensitivity experiments (a) combining MOD climate/LGM CO₂ and (b) combining LGM climate/MOD CO₂ levels. The empirical models use climate, vegetation, topography, lightning ignitions, land cover, road density, and human population density as predictors to represent the environmental controls on each of the wildfire properties.

Figure 1 Flowchart of the method to obtain each of the four scenarios: MOD climate/MOD CO₂, LGM climate/LGM CO₂, MOD climate/LGM CO₂, and LGM climate/LGM CO₂.

Modern (MOD) climate data (daily temperature (T), daily precipitation (P), photosynthetic photon flux density (PPFD), monthly wind speeds (wind), vapour pressure deficit (VPD), monthly specific humidity (huss), cloud cover (cld), monthly pressure (Pa)) were obtained from the WFDE5 bias-adjusted ERA5 database (Cucchi et al., 2020) for 2010 to 2015. The number of monthly dry days (DDs) (days with ≤ 1 mm of precipitation); the monthly diurnal temperature range (DTR) (daily maximum temperature – daily minimum temperature); and the monthly vapour pressure deficit (VPD), which is a function of specific humidity, temperature, and pressure, were all calculated following the methodology in Haas et al. (2022). Seasonal climatologies were derived for all variables eliminating inter-annual variability. For each grid cell, values from the month with (on average) the maximum number of DDs, the largest DTR, and the highest VPD were selected. Wind speed value was taken from the hottest month of the year (determined from the WFDE5 2 m air temperature; Cucchi et al., 2020). For lightning, the mean value over the seasonal climatology was selected. A seasonality predictor to account for wet vs. dry seasons was constructed by dividing the range of monthly values from the seasonal DD climatology by the mean value of all 12 months. Expanded ice sheets in North America, Fennoscandia, Greenland, and Antarctica resulted in global sea levels ∼ 120 m lower than today at the LGM. The modern climate data were extrapolated out onto the exposed shelves using the ICE-6G_C (Peltier et al., 2015) boundary conditions and a nearest-neighbour approach from the GeoInterpolation package in R.

Figure 2 Latitudinal distribution of the LGM-MOD climate anomalies for MPI.ESM1.2 (orange), AWI-ESM1.2 (pink), and CESM1.2 (brown) for (a) the maximum number of dry days, (b) the seasonality of dry days, (c) maximum monthly VPD, (d) maximum monthly DTR, (e) maximum monthly mean wind speeds, (f) mean monthly temperature, and (g) mean monthly total precipitation. The zero-intercept line represents no change between LGM and MOD climate, with negative values representing lower values at the LGM and positive values representing higher values at the LGM.

LGM climate data were obtained from three models participating in the Palaeoclimate Modelling Intercomparison Project (PMIP) contribution to the sixth phase of the Coupled Model Intercomparison Project (CMIP6), namely AWI-ESM-1-1-LR (short name: AWIESM1) (Lohmann et al., 2020; Sidorenko et al., 2015), MPI_ESM1.2 (Mauritsen et al., 2019), and CESM1.2 (Li et al., 2013; Tierney et al., 2020), to represent a range of LGM climates (Fig. 2). A seasonal climatology was derived for each climate variable from the PMIP picontrol experiment (pre-industrial conditions, PI) and the PMIP lgm experiment of the PMIP4 and CMIP6 simulations. The difference between the PI and LGM values (LGM-PI climate anomalies) were calculated and added to the MOD climatology (LGM-MOD climate anomalies) (see Fig. 1). We use the term climate anomalies to refer to the difference between the MOD climatology for each climate variable and the computed bias-adjusted LGM climatology for the same variable, consistently with the PMIP4 protocol (Kageyama et al., 2021). The use of anomalies is designed to minimize the impact of systematic model biases on the derived climate. This approach provided three LGM climate scenarios, resulting in 12 experiments for BA, FS, and FI.

We obtained MOD and LGM vegetation and gross primary production (GPP) using the coupled biogeography and biogeochemistry model BIOME4 (Kaplan et al., 2003) and a simple optimality based model of GPP, the P Model (Wang et al., 2017; Stocker et al., 2020). BIOME4 was used to simulate biome distribution with modern-day climate data (T, P, cld), setting CO₂ levels to 395 ppm (the annual mean from 2010–2015) and 185 ppm in turn. LGM biome distributions were simulated using the three different LGM scenarios, again setting CO₂ levels to 395 and 185 ppm. We derived mean fractional tree, shrub, and grass cover for each of these 12 experiments using the mean values for each biome from ESA CCI (Climate Change Initiative) Land Cover (Li et al., 2018). We also calculated fAPAR (fraction of absorbed photosynthetically active radiation) for each experiment from the leaf area index (LAI) computed by BIOME4 and obtained the fractional cover of C₄ plants (see Sect. S1 in the Supplement). We computed global monthly C₃ and C₄ photosynthesis using the P model using appropriate combinations of climate (T, VPD, ppfd, and Pa), BIOME4-derived fAPAR, and CO₂ concentration for the MOD and LGM scenarios (see Fig. 1). Total GPP was calculated as follows:

$$\begin{array}{}\text{(1)}& {\mathrm{GPP}}_{\mathrm{monthly}}={\mathrm{GPP}}_{{\mathrm{C}}_{\mathrm{3}}}\left(\mathrm{1}-{\mathrm{C}}_{{\mathrm{4}}_{\mathrm{fraction}}}\right)+{\mathrm{GPP}}_{{\mathrm{C}}_{\mathrm{4}}}{\mathrm{C}}_{{\mathrm{4}}_{\mathrm{fraction}}},\end{array}$$

with ${\mathrm{GPP}}_{{\mathrm{C}}_{\mathrm{3}}}$ and ${\mathrm{GPP}}_{{\mathrm{C}}_{\mathrm{4}}}$ representing monthly C₃ and C₄ GPP values from the P model and ${\mathrm{C}}_{{\mathrm{4}}_{\mathrm{fraction}}}$ representing the fractional C₄ cover from BIOME4 (see Table 1).

Table 1 Total annual gross primary production (GPP) (in PgC) estimates for each scenario.

This approach led to estimates of total BA, median FS, and median FI under modern conditions of similar magnitudes to the original GLM models and other global estimates (Andela et al., 2019; Humber et al., 2019) (Table 2).

Topographic and lightning variables were assumed not to change between the LGM and the present day. We used modern values, extrapolated out onto the exposed shelves, for the LGM experiments. The GLMs (Haas et al., 2022) include predictors associated with human activity, specifically human population density, road density, and cropland cover. Population density is used as a measure of potential human ignitions, and road density and cropland cover are used as measures of landscape fragmentation. Including these anthropogenic predictors in the GLM models was found to be essential to capture the global drivers of the observed spatial patterns of wildfires (Haas et al., 2022). This is because modern fire regimes are influenced by human activity at a global scale (e.g. Marlon et al., 2008; Bowman et al., 2020; Harrison et al., 2021). However, although the practice of foraging for plants by some hunter-gatherer communities at the LGM has been shown (Liu et al., 2013), we presume that there was no large-scale agriculture (or road networks) at the LGM. Additionally, information about pre-agricultural population sizes is limited and highly uncertain (see e.g. A. N. Williams et al., 2013; Gautney and Holliday, 2015), and though some regional models of human population do exist (Tallavaara et al., 2015), a reliable global product is not yet available. To avoid confounding effects due to the high uncertainty of human impacts on global wildfire regimes, we decided to exclude these anthropogenic predictors in all the experiments by setting them to 0. This ensured that differences between the experiments were driven solely by climate and CO₂. We performed sensitivity analysis to examine the impact of setting human predictors to 0 under modern and LGM conditions (see Sect. S2). Whilst BA and FS increase in the modern sensitivity analyses (especially in areas with high road density and cropland density such as Europe and India), the effect was negligible for FI, highlighting the sensitivity of BA and FS to human activity. Under LGM conditions, the effect of including human population was negligible for all three fire properties. This reflects the slight and localized human impact on the natural landscape at the LGM (Black et al., 2007; Fuller et al., 2014; Portenga et al., 2016).

When modelled GPP values were 0, BA, FS, and FI were automatically set to 0. Modelled BA values smaller than 0.001 were assumed to imply no burning; thus, under these conditions, FS and FI were also assumed to be 0 since both GLM models were trained on data of existing fires (see Sect. S3).

The resulting BA, FS, and FI anomalies refer to the difference between the MOD climate/MOD CO₂ experiment and the three other experiments since each experiment is considered to represent the long-term average spatial pattern for each fire property under the set experimental conditions. We used the sensitivity experiments to quantify the separate effects of CO₂ and climate on BA, FS, and FI independently. We then used the realistic experiments to identify which predictors were driving the largest change between MOD and the three LGM scenarios by excluding one predictor at a time from the GLM models, re-running the LGM experiments, and identifying which excluded variable caused the greatest change in the BA, FS, and FI MOD-LGM anomalies in each grid cell. Comparing these results to the BA, FS, and FI MOD-LGM anomalies of the full GLM models allowed us to determine whether the predictor was responsible for an increase or a decrease in BA, FS, and FI.

We also compared the spatial patterns of BA, FS, and FI with sedimentary charcoal data from the Reading Palaeofire Database (RPD; Harrison et al., 2022). Sedimentary charcoal records provide a record of fire activity but may reflect changes in both burnt area or completeness of combustion (Power et al., 2008), so this comparison allowed us firstly to establish which of the fire regime properties was most closely reflected in these records and secondly which of the scenarios produced the most realistic patterns of burning. Model outputs and the charcoal records were re-gridded to the coarsest resolution of the three climate models (2.5^∘ × 1.875^∘ resolution). We calculated the number of correctly predicted BA, FS, or FI anomalies (same sign within a given grid cell), separating positive and negative BA, FS, or FI anomalies to assess the rate of false positives and false negatives for each scenario and each LGM climate scenario.

3 Results

Global BA was substantially reduced compared to the realistic MOD scenario under all three realistic LGM scenarios, decreasing by 72 % for the coldest CESM1.2 LGM scenario, 62 % for the MPI-ESM1.2 LGM scenario, and 41 % for the warmest AWIESM1 LGM scenario (see Fig. 3). The largest decreases were observed in sub-Saharan Africa (excluding the tropical regions), as well as in northern Australia and the Indian subcontinent (MPI-ESM1.2 and CESM1.2 LGM scenarios). Some increases in BA were observed in Alaska (MPI-ESM1.2 and AWIESM1 LGM scenarios), as well as in southeast Asia, Indonesia, Papua New Guinea, and the northern tip of Australia. Increases in Somalia and Central America were also observed (MPI-ESM1.2 and AWIESM1 LGM scenarios). The number of grid cells (excluding ice-covered cells) in which no burning occurred was 3 times higher in the MPI-ESM1.2 and AWIESM1 LGM scenarios and 4 times higher in the CESM1.2 LGM scenario compared to the realistic MOD scenario. This was driven by the expansion of desert and tundra biomes at the LGM. The Arabian plate, the Middle East, inland China, and Australia, as well as the tips of South America and Africa, saw burning reduced to 0. Nearly all burning above 60^∘ N was excluded, except for that in Alaska under the MPI-ESM1.2 and AWIESM1 LGM scenarios, with the exclusion extending down to 50^∘ N for the CESM1.2 LGM scenario (see Sect. S3).

Figure 3 Experiments for BA, FS, and FI for MPI-ESM1.2 LGM scenario (MOD 395 ppm and LGM 185 ppm represent the realistic modern-day simulation and LGM simulation, whilst MOD 185 ppm and LGM 395 ppm represent the CO₂ and climate sensitivity experiments respectively. The ice is shown in grey). The other experiments can be found in Sect. S3.

Globally, there was a large decrease in global median FS and FI when considering all grid cells (not covered in ice) because of the overall global reduction in burning (see Fig. 3). Under all three LGM scenarios, global median FS and FI were reduced to 0 compared to ∼ 5 km² for FS and 40 W km⁻¹ for FI. However, when excluding grid cells in which no burning occurred, global median FS increased compared to the realistic MOD scenario (by ∼ 16 % under the two less conservative scenarios (MPI-ESM1.2 and AWIESM1) and by 12 % under the CESM1.2 LGM scenario). The main increases in FS occurred in Central America, Amazonia, tropical Africa, the Indian subcontinent, Europe, and Asia between 30 and 60^∘ N (except for CESM1.2, which had very few positive FS anomalies). The largest reductions were observed in North America, southern Australia, the Middle East, and the rest of Eurasia. Global median FI also increased in regions that were burning under two of the LGM scenarios – by 11 % under the CESM1.2 LGM scenario and by 4 % under the MPI.ESM1.2 LGM scenario. Under the AWIESM1 LGM scenario, global median FI decreased by 2 % even when excluding grid cells that were not burning. Despite this, changes in FI were spatially consistent across all three LGM scenarios, with increases in FI occurring primarily across the American and African continents, as well as in the Mediterranean Basin and Europe, and with decreases occurring in Asia and inland Australia.

Under low CO₂ levels with MOD climate (MOD climate/LGM CO₂), global BA decreased by ∼ 70 % under all three LGM scenarios (72 % for CESM1.2 and AWIESM1, 73 % for MPI-ESM1.2). Despite larger global decreased BA compared to the realistic LGM scenarios, the number of grid cells in which no burning occurred was only 1.7 times higher for the MPI-ESM1.2 and AWIESM1 LGM scenarios and 1.5 times higher for the CESM1.2 LGM scenario compared to the realistic MOD scenario. The spatial pattern was consistent across all three LGM scenarios, with very few grid points showing a positive BA anomaly relative to the MOD experiment. Though FS increased slightly under this sensitivity experiment when burning did occur, this increase was concentrated in the tropical regions of South America and Africa (mainly Amazonia), except in the case of AWIESM1, where increases were observed across Eurasia. In burning grid cells, global median FI increased by ∼ 15 %–18 % in this sensitivity experiment (18 % for MPI-ESM1.2 and CESM1.2 and 15 % for AWIESM1). This spatial pattern was also consistent with BA, with very few negative FI anomalies, except for regions ∼ 20–30^∘ N and ∼ 20–30^∘ S.

Under MOD CO₂ and LGM climate, BA decreased by 41 % compared to the MOD experiment for the CESM1.2 LGM scenario and by 4 % for the MPI-ESM1.2 LGM scenario, but it increased by 48 % for the AWIESM1 LGM scenario, showing a strong sensitivity to climate. The number of grid cells in which no burning occurred was of similar amplitude to the previous sensitivity experiment for the MPI-ESM1.2 and AWIESM1 LGM scenarios (∼ 1.8 times higher compared to the realistic MOD scenario) but was much higher for the CESM1.2 LGM scenario (∼ 3.5 increase). When burning occurred, the global median FS increased under all LGM scenarios by 17 % for CESM1.2, 25 % for MPI-ESM1.2, and 23 % for AWIESM1. These increases were concentrated in tropical Africa, central America, and Russia, with decreases shown in North America and South Africa. Global median FI also increased under this sensitivity experiment by 2 %–3 % for AWIESM1 and MPI-ESM1.2 but decreased by 5 % for the CESM1.2 LGM scenario, with decreases concentrated in Eurasia and North America.

Figure 4 Boxplots showing relative importance of each predictor (GPP – gross primary production; GPP.s. – GPP seasonality; tree – tree cover; shrub – shrub cover; grass – grass cover; DDs – dry days; DD.s. – dry-day seasonality; VPD – vapour pressure deficit; DTR – diurnal temperature range; wind – wind speed) in driving the BA, FS, or FI anomaly between the MOD 395 ppm and LGM 185 ppm experiment. For each grid cell common to both experiments (on modern-day continental shelves and masking the LGM ice sheets), the predictor which caused the largest change in the anomaly between the two experiments when it was excluded from the GLM model was retained; it is the change in anomaly that is shown here. This was taken as an indicator of the relative importance of that predictor in driving the observed change for (d–f) the AWIESM1 LGM scenario, (a–c) the MPI-ESM-1.2 LGM scenario, and (g–i) the CESM1.2 LGM scenario. A positive anomaly indicates that the variable caused an increase in BA, FS, or FI at the LGM, and a negative anomaly indicates that the variable caused a decrease in BA, FS, or FI at the LGM.

Reductions in BA between the MOD and LGM scenarios were driven primarily by changes in GPP; grass cover; VPD; and, to a lesser extent, dryness (dry days (DDs) and dry-day seasonality (DD.s)). Changes in FI were driven by changes in GPP and VPD, with changes in GPP seasonality also leading to increased FI in inland regions, reflecting both changes in climate and CO₂ levels for BA and FI. Increased FS was largely driven by increased wind speeds; DDs and diurnal temperature range (DTR), reflecting a strong climate effect; and GPP seasonality. Decreases in FS were driven by changes in GPP and grass cover, as well as by VPD under the CESM1.2 LGM scenario and DTR under the AWIESM1 LGM scenario (Fig. 4). Changes in GPP and grass cover were responsible for the largest reductions in burning, with these vegetation effects being concentrated across Africa and much of Eurasia (see Fig. 5). In Amazonia, changes in DDs were the most important factor, reducing BA and FS (except for MPI-ESM1.2, which saw increased FS driven by DDs). Increased BA in western Alaska was driven by GPP in the MPI-ESM-1.2 and AWIESM1 LGM scenarios. Increased BA in tropical regions was driven by grass cover, GPP, and DD changes. Changes in VPD across the northern latitudes, especially of North America and Europe, led to decreased BA in the most conservative CESM1.2 LGM scenario. FS decreased across the Americas and Eurasia in the CESM1.2 LGM scenario because of low VPD values, which reduced the occurrence of burning and offset the increases caused by wind speed and DTR in the other two LGM scenarios. Low values of VPD drove increases in FI across eastern North America, South America, western Africa, and southeast Asia.

Figure 5 Map showing selection of four variables (GPP in green, DDs in red, VPD in orange, and wind in grey) responsible for some of the most important grid cell drivers in reducing BA and increasing FS and FI for (d–f) the AWIESM1 LGM scenario, (a–c) the MPI-ESM1.2 LGM scenario, and (g–i) the CESM1.2 LGM scenario. Maps of most important grid cell drivers for all variables and all experiments can be found in Sect. S3.

Comparing the spatial patterns of the simulated BA anomalies with charcoal-based reconstructions of the sign of changes in biomass burning (RPD; Harrison et al., 2022) showed that the best overall match occurred when both the climate and CO₂ effect were considered, with a success rate of ∼ 39 %–45 % depending on the climate scenario. The MPI-ESM1.2 and AWIESM1 LGM scenarios produced the best overall matches. None of the MOD climate/LGM CO₂ experiments identified any of the positive BA anomalies shown by the charcoal records. The LGM climate/MOD CO₂ experiments identified around half (∼ 10 %–17 %) of the negative BA anomalies identified by the realistic experiment (17 %–20 %) and the MOD climate/LGM CO₂ sensitivity experiment and only performed marginally better than the realistic experiment in identifying the positive BA anomalies (Table 3). Thus, although this sensitivity experiment produced a similar overall agreement with the reconstructions as the LGM climate/LGM CO₂ simulations, only the realistic scenarios produced similar success rates for both the negative and positive BA anomalies. Climate change alone produced too few negative anomalies matches; CO₂ changes alone resulted in no positive anomaly matches.

Figure 6 Comparison of BA anomalies between the experiment outputs from the MPI-ESM1.2 LGM scenario with charcoal records from the Reading Palaeofire Database (RPD) for (a) the realistic LGM experiment (b) the LGM climate/MOD CO₂ sensitivity experiment and (c) the MOD climate/LGM CO₂ sensitivity experiment The modelled positive LGM-MOD anomalies are shown in red and LGM-MOD negative anomalies in blue. Dotted red (positive anomaly) and blue (negative anomaly) points show the location of the RPD records for the LGM. The LGM ice sheets are shown in dark blue.

Table 2 Comparison of sign in BA anomalies (between the MOD climate/MOD CO₂ experiment and other three experiments respectively) at the location of each RDP (Reading Paleofire Database) charcoal-based reconstruction record. A positive anomaly represents increased biomass burning, and a negative anomaly represents decreased biomass burning. A successful identification means that the sign of the experiment anomaly and the sign of the RPD charcoal-based reconstructions are the same.

The sign of the charcoal records could reflect changes in FS or FI, as well as in BA. However, the success rates in predicting the sign of the charcoal anomalies (both positive and negative) were not as good for FS (27 %–31%) and FI (24 %–30 %) as those obtained for BA for the realistic LGM experiment. Furthermore, both FS and FI did not perform any better than BA under any experiment, with the sensitivity experiments matching the charcoal anomalies slightly better for FS and FI than under the realistic LGM experiment (see Sect. S4).

4 Discussion

Our simulations show a global reduction in burning at the LGM but increased median fire size and intensity when burning did occur. BA, FS, and FI were all sensitive to changes in vegetation driven directly by CO₂ levels alone. BA and FI were most sensitive to this effect, with the climate effect dampening the effect of CO₂ alone when both are included. The largest reductions in burning occurred when only the CO₂ effect was considered, although this experiment had fewer regions in which burning was excluded completely. This suggests that the reduction in burning was more spatially consistent and widespread under these conditions than when both effects were accounted for. The sensitivity of BA to CO₂ is explained by the reduction in fuel availability under low-CO₂ conditions, a strong constraint on burnt area. For FI, including a CO₂ effect also amplified the overall global signal. This CO₂ effect is most likely driven by the negative relationship between GPP and FI fitted by the empirical model. Whilst this relationship might seem counter-intuitive, it has a sound basis. The most intense fires occur in regions with a seasonal variation in productivity rather than in the most productive environments such as tropical forests (Archibald et al., 2013). High productivity can (under some climate conditions) increase the frequency of burning, which also reduces fuel loads (Rodrigues et al., 2019). Under appropriate climate conditions, there can be long-term fuel build-up in areas of low productivity that is not offset by frequent burning. All these factors help to explain why FI is not reduced at the LGM when burning occurs even though BA is. Low CO₂ decreased FS, except in tropical regions, and reduced the impact of climate in the realistic scenarios. We hypothesize that this is because of decreased productivity leading to patchier vegetation and hence reduced fuel continuity, which is a factor limiting wildfire spread (Dial et al., 2022; Schertzer et al., 2015).

Changes in climate alone also affected all three modelled wildfire properties. The climate effect was larger than the CO₂ effect across all models for FS, with increases in wind, DD, and DTR driving the change. BA was particularly sensitive to the amplitude of climate change: climate change alone greatly reduced BA under the coldest LGM scenario (CESM1.2), had a limited effect in the intermediate LGM scenario (MPI-ESM1.2), and increased BA in the warmest LGM scenario (AWIESM1). The amplitude of change in VPD, a measure of atmospheric moisture, relative to other climate variables was especially important in influencing overall trends. In the case of BA, large decreases in VPD under the CESM1.2 climate scenario led to much more substantial reductions, most likely due to an increase in fuel moisture. Additionally, though stronger winds and increased DTR were the main drivers of larger wildfires at the LGM, low VPD values in CESM1.2 severely limited FS and FI in the northern latitudes. VPD has been shown to influence wildfire ignition and wildfire spread (Sedano and Randerson, 2014), and our results suggest that high atmospheric moisture can inhibit fire spread. When vegetation was sufficiently abundant, however, low VPD values were key in driving intensity. Although vegetation productivity was lower at the LGM, decreased VPD may have contributed to larger fuel build-ups, thus increasing fuel loads. This highlights the sensitivity of the fire regime not just to overall climate change but to the relative amplitude of change in individual climate variables.

Our model results reproduce the global reduction of biomass burning at the LGM observed from ice cores and sedimentary charcoal records (Daniau et al., 2012; Harrison et al., 2022; Power et al., 2008; Rubino et al., 2015). Some studies have indicated the occurrence of high-intensity wildfires on the Palaeo-Agulhas Plain of South Africa and in tropical regions, northern Australia, and central China at the LGM (Kraaij et al., 2020; Power et al., 2008; Rowe et al., 2021; Ruan et al., 2020; Song et al., 2023). Our results are consistent with the trends in these regions. The LGM simulations of BA that account for both climate and CO₂ appear to fit the charcoal records best. The spatial patterns of BA at the LGM were more consistent with the patterns shown by sedimentary charcoal records than FS and FI, consistent with the assumption that charcoal abundance can be used as a measure of biomass burning. The FS and FI anomaly patterns were less consistent than that of BA, suggesting that a regime of less burning but larger and more intense wildfires at the LGM could be consistent with the charcoal records. Whilst FI has been reconstructed from charcoal (e.g. Duffin, 2008; Snitker, 2018), there are currently no comparable measures that record FS or FI changes globally. Charcoal records are not available from some regions, further limiting our ability to evaluate the models, particularly in Eurasia and inland South America, where low CO₂ leads to large reductions in BA that are not observed when only climate is considered.

Our results are based on simple empirical models for BA, FS, and FI. However, the inferred changes in BA are like those of Martin Calvo et al. (2014), who used the Land surface Processes and eXchanges (LPX) dynamic global vegetation model. Empirical models have been shown to perform as well as more complex process-based models in simulating burned area under modern-day conditions (Hantson et al., 2020). Thus, our conclusions about the relative impact of climate and CO₂ changes on fire properties are unlikely to be adversely affected by the relative simplicity of the models used. Their simplicity facilitates running multiple scenarios and diagnosis of the factors influencing changes in wildfire properties.

The effect of human activity was not considered in this analysis, and, as such, no conclusions can be drawn on how human activity may affect these trends. Although this is a limitation, we believe it is unlikely that human activity would substantially impact the response of wildfire regimes to the changes in climate and CO₂ observed here. Pre-agricultural hunter-gatherer populations used fire for land management, for example to facilitate hunting and to promote the local abundance of food plants (Bowman, 1998; Gott, 2005), although recent work indicates that the burning regimes they practised tended to reduce fire overall compared to the natural state (see e.g. Constantine et al., 2023). However, the areas suitable for hunter-gatherer populations was much reduced at the LGM by generally colder and drier climates, and hunter-gatherer populations were confined to climatically suitable refugia (see e.g. A. N. Williams et al., 2013; Blinkhorn et al., 2022). Furthermore, although the estimates of population density are highly uncertain, the LGM population of Australia was less than 5 % of the modern population, and the reduction in Africa was even larger (Gautney and Holliday, 2015). Palaeoecological evidence from Australia suggests that the use of fire by pre-agricultural hunter-gatherers had a low impact on the environment before the late Holocene (e.g. Black et al., 2007; Fuller et al., 2014; Portenga et al., 2016). Thus, it is unlikely that human activities during the LGM would have substantially increased fire or offset the impact of the changes in climate and CO₂ on fire regimes. Previous studies show a weak influence of population and land use change in driving global wildfire trends prior to the 18th century (e.g. Pechony and Shindell, 2010; Bowman et al., 2020) and a sharp human-driven decline in wildfire activity since the mid-nineteenth century (e.g. Marlon et al., 2008; Wang et al., 2010). This recent reduction in global biomass burning was most likely driven by population growth and land use change leading to increased landscape fragmentation, which tends to suppress fire spread (e.g. Knorr et al., 2014; Andela et al., 2017; Harrison et al., 2021).

These results add to a growing body of literature highlighting the importance of considering not only changes in wildfire weather but also vegetation properties in projections of future wildfire regimes (e.g. Harrison et al., 2021; Kuhn-Régnier et al., 2021; Pausas and Keeley, 2021). The impact of rising CO₂ levels will most likely enhance vegetation growth and litter accumulation, which are important controls on fuel availability, continuity, and load. However, climate and specifically VPD may have opposing effects to that of rising CO₂ levels. Since VPD controls plant growth, increasing VPD can limit ecosystem productivity and tree growth, in turn reducing fuel loads (A. P. Williams et al., 2013). Nevertheless, VPD has also been shown to increase litter fall, thus increasing available dead fuel (De Dios et al., 2021; De Faria et al., 2017). As such, it is important to consider how temporal and spatial scales affect the response of vegetation to changing VPD (Grossiord et al., 2020). Although the trade-offs between future increases in CO₂ and reductions in productivity due to higher temperatures and atmospheric dryness are not fully understood, this work highlights the importance of considering both. These effects will most likely not be evenly distributed across the globe (Gonsamo et al., 2021; Piao et al., 2020; van der Sleen et al., 2015), and CO₂ effects may be more important in some regions than others. In fuel-limited ecosystems, CO₂ fertilization could increase fuel loads and fuel continuity, increasing overall burnt area but also the potential for larger and more intense wildfires. This is particularly worrying in regions with anticipated decreases in atmospheric moisture, especially since evidence suggests rising VPD may only counteract a small proportion of CO₂-induced plant growth (Song et al., 2022). Increased woody thickening, for example in tropical southern Asia (Kumar et al., 2021; Scheiter et al., 2020), may also alter fuel loads in regions that are likely to be vulnerable to ignition under a drier and warmer atmosphere (Clarke et al., 2022). Whilst climate variables such as DDs and DTR have also been shown to be strong controls of global wildfire regimes (e.g. Bistinas et al., 2014; Forkel et al., 2019; Kuhn-Régnier et al., 2021), this study highlights the importance of VPD relative to other climate variables in driving spatial patterns of BA, FS, and FI. This is in line with previous studies that have highlighted the important role of VPD in promoting fuel loads and fire spread (e.g. Diffenbaugh et al., 2021; Grillakis et al., 2022; Duane et al., 2021; Balch et al., 2022). Correctly projecting changes in fuels in the next century will require considering both the effect of VPD and the effects of CO₂ on plant growth and fuel loads.

Our results stress the importance of accounting for the effects of CO₂ on vegetation when considering how future fire regimes may evolve. Different aspects of the fire regime respond differently to changes in fuel properties. Without accounting for this crucial effect, our understanding of future risks will remain limited.