Reconstructing burnt area during the Holocene: an Iberian case study

Charcoal accumulated in lake, bog or other anoxic sediments through time has been used to document the geographical patterns in changes in fire regimes. Such reconstructions are useful to explore the impact of climate and vegetation changes on fire during periods when human influence was less prevalent than today. However, charcoal records only provide semi-quantitative estimates of change in biomass burning. Here we derive quantitative estimates of burnt area from vegetation data in two stages. First, we relate the modern charcoal abundance to burnt area using a conversion factor derived from a generalised linear model of burnt area probability based on eight environmental predictors. Then, we establish the relationship between fossil pollen assemblages and burnt area using tolerance-weighted weighted averaging partial least-squares regression with a sampling frequency correction (fxTWA-PLS). We test this approach using the Iberian Peninsula as a case study because it is a fire-prone region with abundant pollen and charcoal records covering the Holocene. We derive the vegetation–burnt area relationship using the 31 records that have both modern and fossil charcoal and pollen data and then reconstruct palaeoburnt area for the 113 records with Holocene pollen records. The pollen data predict charcoal-derived burnt area relatively well (R² = 0.44), and the changes in reconstructed burnt area are synchronous with known climate changes through the Holocene. This new method opens up the possibility of reconstructing changes in fire regimes quantitatively from pollen records, after regional calibration of the vegetation–burnt area relationship, in regions where pollen records are more abundant than charcoal records.

1 Introduction

Fire is an important element in many ecosystems and in the Earth system (Bowman et al., 2009; Resco de Dios, 2020). It impacts vegetation dynamics, ecosystem functioning and biodiversity (Harrison et al., 2010; Ward et al., 2012; Keywood et al., 2013). It also affects climate through vegetation changes and the release of trace gases and aerosols. Fire directly impacts socio-economic assets (e.g. Stephenson et al., 2013; Thomas et al., 2017) and has deleterious effects on human health though the release of smoke and particulates into the atmosphere (e.g. Johnston et al., 2012; Yu et al., 2020). These impacts make it important to understand what controls the incidence and severity of fires.

Analyses of fire regimes during the satellite era have shown that multiple factors play a role in determining the occurrence of fire, including climate and fire weather, vegetation properties and human activities (e.g. Harrison et al., 2010; Brotons et al., 2013; Bistinas et al., 2014; Knorr et al., 2014; Andela et al., 2017; Forkel et al., 2019a, b; Kuhn-Régnier et al., 2020). However, the satellite record only covers a short time period (ca. 20 years), and the impact of anthropogenic changes to land cover on suppressing fire during this interval is strong (Andela et al., 2017). Reconstructing changing fire regimes during the pre-industrial Holocene (12 000 BP to ca. CE 1850) provides an opportunity to investigate the controls on fire over timescales when human influences on the landscape, including fire regimes, were more localised and less profound than they have become during the industrial era.

Sedimentary charcoal, preserved in lakes, peatbogs and other anoxic environments, has been widely used as an indicator of past changes in fire regimes (Marlon et al., 2008, 2016; Power et al., 2008; Daniau et al., 2012; Vannière et al., 2016; Connor et al., 2019). Evaluations that combine charcoal-inferred palaeofire reconstructions with past hydrological, vegetation and archaeological data support the idea that there are strong relationships among climate, fire, vegetation and human activities (Carrión et al., 2007; Marlon et al., 2008; Gil-Romera et al., 2010; Turner et al., 2010; Vannière et al., 2011; López-Sáez et al., 2018; Morales-Molino et al., 2018). However, charcoal records only provide a semi-quantitative index of fire activity rather than quantitative estimates of burnt area or biomass loss. Attempts to calibrate the charcoal record to provide quantitative estimates of proximity or area burnt are either site-specific (Duffin et al., 2008; Hennebelle et al., 2020) or rely on modelling (Higuera et al., 2007). Furthermore, although the number of charcoal records is increasing, there are still comparatively few sites compared to other types of palaeoenvironmental data, and this can make it difficult to make regional reconstructions of changing fire regimes.

Although the occurrence of fire is influenced by multiple factors, analyses of present-day fire relationships globally using satellite-derived data have shown that vegetation properties determining fuel availability are the strongest determinants of fire occurrence (Bistinas et al., 2014; Forkel et al., 2019a, b; Kuhn-Régnier et al., 2020). This suggests that palaeovegetation data could provide a way of reconstructing burnt area in the past, particularly at times when human influences on land cover were less important. This would also allow us to capitalise on the more extensive site networks for palaeovegetation.

In this study, we present a new method to reconstruct quantitative changes in fire regimes over the Holocene. We relate the relative scale of modern charcoal abundance to absolute burnt area using a conversion factor derived from a generalised linear model (GLM) of fire probability based on burnt area data. We then derive quantitative relationships between pollen assemblages and inferred burnt area using tolerance-weighted weighted averaging partial least-squares regression with a sampling frequency correction (fxTWA-PLS; Liu et al., 2020). The vegetation–burnt area relationship is then used to reconstruct changes in burnt area through time from pollen assemblages, including at sites with no charcoal record. We use the Iberian Peninsula as a test case. The Iberian Peninsula is the most fire-affected region in southern Europe (Jesus et al., 2019; Molina-Terrén et al., 2019). Although the modern fire regime is partly driven by human activities, the patterns also reflect the strong climate gradients across the region. Although much of the Iberian Peninsula has a typical Mediterranean climate, parts of the region are influenced by proximity to the Atlantic Ocean or the Mediterranean Sea and by the mountainous topography, giving rise to complex weather and climate patterns and large gradients in vegetation diversity (Loidi, 2017). We reconstruct fire regimes across the Iberian Peninsula through the Holocene and discuss the implications of the reconstructed changes.

Figure 1 Flow chart of the methodology of burnt area reconstructions.

2 Methods

The central premise of our approach is that fire frequency is one of the factors that influences vegetation assemblages (see the Supplement) and therefore that specific aspects of differences in vegetation assemblages – identified by a numerical technique that can isolate the effects of any one controlling factor on taxon composition – can be used to reconstruct fire. The vegetation–fire relationship can be derived by comparing changes in pollen assemblages and charcoal records through time. However, since the charcoal records from different sites consist of different size fractions, and the records must be normalised to facilitate comparisons, it is necessary to derive site-specific conversion factors between modern charcoal abundance and present-day burnt area fraction. This calibration is then applied to the charcoal record in order to derive an estimate of the palaeoburnt area for each pollen sample.

2.1 Iberian pollen and charcoal data

Pollen data were obtained from the European Pollen Database (EPD; http://www.europeanpollendatabase.net, (last access: 5 April 2021) or provided by the authors (Table S2 in the Supplement). Non-pollen palynomorphs (e.g. fungi and algae), introduced species and fire-insensitive plants (e.g. obligate aquatics) were removed from the assemblages before analysis. Some pollen taxa are not identified consistently by palynologists or occur at very few sites, so some pollen types were amalgamated to higher taxonomic groups (mostly genera for trees and families for herbaceous taxa) for consistency across the records (Table S3). Charcoal data were obtained from the Global Charcoal Database (Power et al., 2010; Marlon et al., 2016) or provided by the authors (Table S4). The original age models for both the pollen and the charcoal records were constructed using different methods and different calibrations of radiometric to calendar ages. We created new age models for all the records using the IntCal20 calibration curve (Reimer et al., 2020) and the BACON Bayesian age-modelling tool in the rbacon package (2.5.0) in CRAN (Blaauw and Christeny, 2011). Charcoal concentration data were converted to charcoal accumulation rate (influx: $\mathrm{particles}{\mathrm{cm}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$) before analysis by multiplying concentration by the background sedimentation rate.

2.2 Development of the generalised linear model

We obtained modern burnt area for the Iberian Peninsula from the fourth version of the Global Fire Emissions Database (GEFD4) (Randerson et al., 2017). The GLM was initially developed using 13 environmental variables covering climate, vegetation and human activities (Table S5). Some environmental data sets were only available at 0.5^∘ × 0.5^∘ resolution, so all data sets were aggregated to this resolution using bilinear interpolation prior to analysis. Analyses were done for the common period between the data sets (January 2001 to December 2016) using annual values of all variables. The GLM was run using the stats package in R (version R.3.6.0) and used the logit link function and assumed a quasi-binomial distribution (R Core Team, 2019). We tested combinations of environmental predictors and selected the most parsimonious model with statistically significant variables and high prediction ability as assessed using pseudo-R² (McFadden, 1973). The GLM-fitted burnt area was disaggregated from 0.5^∘ × 0.5^∘ to 0.0083^∘ × 0.0083^∘ by bilinear interpolation in order to extract present-day burnt area at each of the sites with modern charcoal records.

Figure 2 Map showing the location of the 31 entities (31 sites) with modern charcoal used to derive the fire–vegetation relationship (red triangles) and the 113 entities (111 sites) used for burnt area reconstructions (blue circles) in the Iberian Peninsula.

2.3 Quantitative reconstructions of burnt area

We derived the relationship between the pollen assemblage and burnt area using 31 records with modern pollen and modern charcoal (Fig. 2). Rare pollen taxa with fewer than or equal to five occurrences in the data set were removed because they have been shown to have little predictive power in WA-PLS climate reconstructions (Turner et al., 2021). The charcoal records included some sites with only macroscopic and some with only microscopic charcoal. Since this had little impact on the patterns of change through time (Fig. S4 in the Supplement), we used both types, although we used macroscopic charcoal at sites with both macroscopic and microscopic charcoal. The sampling resolution varies between the individual records. To ensure comparability across records, the charcoal and pollen data were temporally binned prior to analysis: the modern bin covers the post-industrial period (CE 1850 to the present); a 100-year bin width was used for earlier intervals. Pollen counts were summed and converted to percentage of the total count in each bin. To standardise the values for different charcoal measurement units, we used a maximum transformation to convert mean charcoal accumulation rates to a 0–1 range (Eq. 1).

$$\begin{array}{}\text{(1)}& x_{i,j}^{\prime }={\displaystyle \frac{x_{i,j}}{x_{j,max}}},\end{array}$$

where $x_{i,j}^{\prime }$ is the transformed value of the ith sample (x_i) in the jth entity. $x_{j,max}$ is the maximum value of all samples in this entity.

The maximum transformation resulted in a similar scale of variability between entities in fire-prone areas and areas with little fire. We therefore applied a conversion factor to rescale the relative charcoal abundance to absolute burnt area for each of the records:

$$\begin{array}{}\text{(2)}& {\text{conversion factor}}_j={\displaystyle \frac{{\text{present-day burnt area fraction}}_j}{{\text{modern charcoal data}}_j}},\end{array}$$

where modern charcoal data are the core-top binned charcoal data in the jth entity, and the present-day burnt area fraction in the jth entity was obtained from the GLM.

The palaeoburnt area fraction for the ith sample in the jth entity was then derived by multiplying the conversion factor of the jth entity by the charcoal value for this sample (Eq. 3):

$$\begin{array}{}\text{(3)}& \begin{array}{rl}{\text{palaeoburnt area fraction}}_{i,j}=& {\text{conversion\hspace{0.17em} factor}}_j\\ & \times {\text{charcoal\hspace{0.17em} data}}_{i,j}.\end{array}\end{array}$$

We applied Box–Cox transformation (Box and Cox, 1964) with λ = 0.25 to the palaeoburnt area fraction in order to reduce skewness prior to the fxTWA-PLS analyses (see the Supplement). The fire–vegetation relationship was determined using the last significant component in fxTWA-PLS, assessed using the p value, to avoid overfitting.

We applied the fxTWA-PLS-derived relationship between pollen abundance and burnt area to the binned pollen data from the 113 pollen records available from the Iberian Peninsula (Fig. 2) to reconstruct changes in fire regimes through the Holocene. Some of these 113 entities included pollen taxa that were not present in the data used to derive the vegetation–burnt area relationship; these taxa were therefore removed prior to analysis. We used composite plots and maps of specific times to show the spatial and temporal changes of reconstructed palaeofire regimes through the Holocene. We used loess smoothing with a window half-width of 300 years to construct the composite plots, with the uncertainty of the reconstruction estimated by bootstrap resampling of the individual reconstructions 1000 times (Efron, 1979; Efron and Tibshirani, 1993). We tested the robustness of our method by comparing the reconstructed burnt area composite with the trends shown by raw charcoal data for those records with fossil charcoal, where the uncertainty is estimated again by bootstrap resampling of the individual charcoal records 1000 times. Maps were created using the reconstructed burnt area for individual sites in the bin covering the time period of interest.

Table 1 Generalised linear model of the modern burnt area fraction.

Notes: ^⋅ p<0.1. ^* p<0.05. ${}^{**}$ p<0.01. ${}^{***}$ p<0.001.

Figure 3 Mean (over 16 years) of observed (a) and fitted (b) values of burnt area fraction.

3 Results

3.1 The GLM

Several of the environmental predictors of burnt area were highly correlated to one another (Fig. S1). We tested the impact of including/removing highly and moderately correlated variables before selecting the final GLM (see the Supplement). The final model was constructed using eight variables (Tables 1 and S6 in the Supplement) and has a pseudo-R² = 0.20. However, all of these variables show statistically significant relationships (P < 0.1) with burnt area, and most have p values < 0.05. Gross primary production (GPP) shows a very strong positive relationship (t = 10.10) with burnt area fraction (Table 1, Fig. S2). Dry days per month (t = 8.46) and non-tree cover (t = 7.34) also show strong positive relationships with burnt area (Table 1, Fig. S2). These relationships make sense given that much of the Iberian Peninsula is relatively arid: increasing GPP and increasing non-tree cover are indices of increased fuel availability in arid, fuel-limited regions and promote increased burnt area. The number of dry days per month determines fuel dryness, and hence there is a positive relationship between the number of dry days and the burnt area. Predictions of burnt area from the final model show reasonably good agreement with the observed average burnt area (Fig. 3). Hotelling's T-squared test shows that there is no statistically significant difference between observed and fitted values (p value = 1). Both the observations and the model show the highest burnt area in northern Portugal and moderate burnt area in southern Portugal. Both observed and simulated burnt area are low along the northern coast and in the Pyrenees where fire is limited by wet conditions and in the dry interior where fire is limited by fuel availability.

Table 2 Leave-out cross-validation fitness of the fxTWA-PLS method, showing results for all the components. The last significant number of components are shown in bold. RMSEP is the root mean square error of prediction. ΔRMSEP is the percent change of RMSEP using the current number of components (ncomp) than using one component less. b₀, b₁, b₀⋅SE, and b₁⋅SE are the intercept, slope, standard error of the intercept and standard error of the slope of the linear regression using the cross-validation result and burnt area data converted from charcoal abundance.

Figure 4 The fitted plot and residual plot of TWA-PLS method, with fx correction. Panel (a) is the reconstructed burnt area using the last significant number of components, which is four here. The x axis is the burnt area fraction derived from charcoal data, and the y axis is the burnt area fraction reconstructed from pollen data using TWA-PLS with fx correction. The 1:1 line is shown in black, and the linear regression line is shown in red to show the degree of overall compression. Panel (b) shows the residuals of reconstructed burnt area fraction using the last significant number of components. The x axis is the burnt area fraction derived from charcoal data, and the y axis is the residual of burnt area reconstruction using TWA-PLS with fx correction. The zero line is shown in black, and the locally estimated scatterplot smoothing is shown in red to show the degree of local compression. The low compression zone is shown by grey shading.

3.2 The pollen–burnt area relationship

The fxTWA-PLS-derived relationship, based on the last (fourth) significant component, has good predictive power (R² = 0.44) (Table 2). A linear regression of the cross-validation results and the burnt area data has a slope of 0.526 (Table 2, Fig. 4a), which shows that the degree of overall compression towards the centre of the sampled range is relatively low. The degree of local compression, which is assessed by whether the residuals are around zero across the burnt area range in locally estimated scatterplot smoothing, indicates that the low-compression zone where reconstructed values after Box–Cox transformation are most reliable is between −3.25 and −2.5, in other words, between 0.12 % and 1.98 % of the grid cell area (Fig. 4b). Comparison with results using WA-PLS and tolerance-weighted WA-PLS (TWA-PLS) confirms that fxTWA-PLS produces a large reduction in compression in the central part of the burnt-area range and has a higher predictive power (Table S8, Figs. S5 and S6). However, although fxTWA-PLS reduces the compression bias, it does not remove it completely: burnt area is overestimated at the low end and underestimated at the high end of the burnt area (Fig. 4).

Figure 5 Composite plots comparing maximum-transformed charcoal values and the reconstructed burnt area for these entities for the 51 entities with charcoal. Maximum-transformed charcoal is shown in blue; burnt area fraction is shown in red. The loess smoothing is done with a span of 0.04.

Charcoal values are not expected to be directly comparable with the reconstructed burnt area but should show comparable temporal trends. A composite plot of reconstructed burnt area for the 51 entities that have both pollen records used to reconstruct burnt area and charcoal records, and therefore can be compared, shows similar trends to the composite plot derived from the maximum-transformed charcoal (Fig. 5). This suggests there is little distortion of the signal caused by deriving burnt area using the fxTWA-PLS relationship.

Figure 6 Composite curve of reconstructed burnt area using fxTWA-PLS, using the locfit() function with a half-width of 300 and 1000 bootstrap samples. The locally estimated scatterplot smoothing is shown in blue. The upper and lower 95th percentile confidence intervals are shown in grey.

Figure 7 Spatial patterns of reconstructed burnt area fraction at key times in the Holocene.

3.3 Fire history of the Iberian Peninsula through the Holocene

The composite plot based on all 113 pollen records from the Iberian Peninsula (Fig. 6) shows a moderate peak in burnt area around 13 ka followed by a marked increase at 12.5 ka and a subsequent peak in burnt area around 11.5 ka. Although this early part of the record is based on relatively few sites, and so the confidence intervals are large, the pattern corresponds to high fire activity during the Bølling–Allerød (14.6–12.9 ka) warm interval and low fire activity during much of the Younger Dryas (12.9–11.7 ka) cold phase, with an increase in burnt area associated with the rapid warming at the end of the Younger Dryas. Burnt area is relatively low at the beginning of the Holocene. Although there is a gradual increase in burnt area between 9 and 0.6 ka, the burnt area fraction is lower than present until at least 2 ka. The increase in burnt area is quite marked after around 4.5 ka and peaks at 0.6 ka. The burnt area fraction at 0.6 ka is larger than at any time in the record. Burnt area declines after 0.6 ka, although the modern reconstructed value is still higher than the values obtained for most of the Holocene.

The spatial coverage of sites (Fig. 7) for the earlier part of the record is sparse, but coverage is good from 7 ka onwards. The pattern of lower burnt area in eastern than in western Iberia, seen in the modern observations, is generally preserved both in high and low fire intervals. However, some of the records from northern Iberia (e.g. Saldropo, Puerto de Los Tornos) show extremely high burnt area which exceeds the scope of the low-compression zone during the last millennium, and particularly at 0.6 ka. This may reflect the persistent bias at the high end of the fx-TWA-PLS reconstruction range.

4 Discussion

We have shown that it is possible to derive trends in burnt area through time by applying a quantitative relationship between pollen assemblages and charcoal-derived burnt area to palaeovegetation records from the Iberian Peninsula. Our analyses exploit the multivariate nature of vegetation and hence pollen assemblages. Vegetation patterns, and the distribution of individual species, are controlled by many factors including seasonal temperature and precipitation regimes, disturbance (including wildfires) and human activities. Pollen-based palaeoclimate methods have long exploited the multivariate nature of pollen assemblages to reconstruct different aspects of climate (see, for example, the discussion in Bartlein et al., 2011). The canonical correspondence analysis (CCA) shows that there is sufficient information in the pollen assemblages to assess the independent contribution of fire to vegetation assemblages. The overall relationship between pollen and charcoal-derived burnt area is reasonably strong (R² = 0.44), reflecting the importance of vegetation properties (gross primary production and non-tree cover) in driving the occurrence of fire – as seen in the GLM analysis of satellite-derived modern burnt area patterns. The overwhelming importance of vegetation properties in influencing modern fire occurrence is consistent with results from global analyses (e.g. Moritz et al., 2012; Pausas and Ribeiro, 2013; Bistinas et al., 2014; Forkel et al., 2019b). Nevertheless, the GLM analysis shows that climate factors, in particular the occurrence of dry intervals, are important controls on modern fire patterns in Iberia. Again, this is consistent with global analyses of the modern drivers of fire occurrence.

We have used fxTWA-PLS (Liu et al., 2020) to make the burnt area reconstructions because this technique reduces the compression bias characteristic of many other reconstruction techniques by accounting for differences in the tolerance of individual taxa in an assemblage and for the frequency of the reconstructed variable in the training data set. However, although the bias is apparently reduced, there is still an overestimation at the low end and an underestimation at the high end of the burnt area range. This is reflected by the extremely high burnt area values reconstructed for some sites in northern Iberia in the recent millennium which exceed the upper limit of the low-compression zone. This remaining bias may also be explained by the comparatively small sample size (1106 binned samples) compared with the much larger data set of 6458 samples used by Liu et al. (2020) to carry out climate reconstructions for Eurasia. It would be useful to test whether the problem of compression bias in the reconstruction of burnt area could be overcome by expanding the training data set to cover a wider range of vegetation types and fire regimes.

Although palaeoburnt area reconstructions have only been obtained from a limited number of records, they nevertheless show interesting patterns over the past ca. 15 kyr. The high fire intervals at the beginning of the record, between 14–13 ka and between 12–11 ka, correspond to the Bølling–Allerød (14.6–12.9 ka) warming interval and to rapid warming at the end of the subsequent Younger Dryas (12.9–11.7 ka) cold phase (Alley et al., 1993). A similar response to these climate events has been seen in charcoal records from eastern North America (Marlon et al., 2009). Burnt area is less than today through the Early and Middle Holocene (10–5 ka), an interval when pollen, speleothem and lake records suggest the Mediterranean region was wetter than today (Prentice et al., 1996; Magny et al., 2002; Bartlein et al., 2011; Roberts et al., 2011). Reconstructions of fire activity anomalies (FAAs) for the south-eastern part of the Iberia Peninsula also indicate low-level fire activity in the Mid-Holocene between 7.5 and 6 ka (Gil-Romera et al., 2010). Burnt area continuously increases during the Medieval Warm Period (MWP; 1–0.7 ka) and peaks at 0.6 ka, consistent with the warm and dry conditions recorded during this period in the Iberian Peninsula (Moreno et al., 2012). During the Little Ice Age (LIA; 0.6–0.1 ka), the reconstructed fire indicates a sharp downturn, which may be associated with subsequent cold and wet climate (Ramos-Román et al., 2016; Abrantes et al., 2017). Thus, the broadscale patterns of trends in reconstructed burnt area are consistent with known Holocene climate changes in this region.

There is a distinct west–east gradient in burnt area across Iberia today, and this gradient of high fire in the west and less fire in the east is also present during other intervals of the Holocene. This pattern is likely related to the regional gradient in fuel availability and drought (Pausas and Fernández-Muñoz, 2012). However, the west–east gradient in burnt area is less pronounced during the Mid-Holocene, consistent with a less pronounced gradient in precipitation and moisture availability shown by other studies (e.g. González-Sampériz et al., 2017; Liu, 2019). Reconstructed patterns in fire were also more homogenous after 1 ka, and again this is consistent with the fact that temperature and humidity gradients were less pronounced at that time than they are today (Sánchez-López et al., 2016).

Our analyses show that climate, and climate-induced changes in vegetation, have influenced the fire regimes of the Iberian Peninsula during the Holocene. However, many studies have suggested that human activities could also have been important (Blanco-González et al., 2018; Connor et al., 2019; Feurdean et al., 2020). Land clearance during the Neolithic agricultural transition has been associated with increases in fire activity in some sites from the Iberian Peninsula (e.g. García-Ruiz et al., 2016; Carracedo et al., 2018). Although initiation of agriculture was not synchronous everywhere, the regional onset of agriculture is registered around 7.5 ka (Zapata et al., 2004; Fyfe et al., 2019; Harrison et al., 2020) when the burnt area reconstructions do not indicate high fire activity. However, the gradual increase in reconstructed burnt area between 5 and 0.6 ka may be an indication of increasing human activity, since the initial increase is broadly consistent with increased population shown by summed probability distributions (SPDs) of radiocarbon dates (Balsera et al., 2015; Lillios et al., 2016; Harrison et al., 2020). Human activities, such as deforestation and appropriation of land for agriculture, may have been an important driver of fire patterns from the Bronze Age onwards (Morales-Molino et al., 2013; Morales-Molino and García-Antón, 2014; González-Sampériz et al., 2017), while the competing effects of land abandonment and fire suppression may have contributed to the changes in burnt area in recent times (Turco et al., 2016; Silva et al., 2019). Nevertheless, our GLM analysis indicates that the intensity of human influence, as measured by crop or grazing land area or by population density, consistently has a negative effect on burnt area under modern conditions. It seems likely that human influence on Holocene fire regimes may have been complex, with agricultural expansion both promoting and suppressing fire occurrence. More detailed comparisons of the reconstructed burnt area and archaeological data are required to test this.

The limited availability of charcoal records has meant that the analysis of past fire regimes has tended to focus on large-scale zonal or continental-scale patterns (e.g. Marlon et al., 2008; Power et al., 2008; Daniau et al., 2010; Vannière et al., 2011). Our new methodology opens up the possibility of reconstructing changes in fire regimes from pollen data and thus of examining finer-scale patterning that might reflect climate or human influences on fire. Spatially explicit reconstructions of burnt area would also be useful to evaluate the simulated response of fire to changing environmental drivers in the past (Thonicke et al., 2005; Brücher et al., 2014; Martin Calvo et al., 2014; Marlon et al., 2016; Kraaij et al., 2020) since comparisons based on qualitative inferences from charcoal are inconclusive (e.g. Brücher et al., 2014).

5 Conclusion

We have developed a novel method to reconstruct palaeoburnt area quantitatively from vegetation records, based on fire–vegetation relationships derived using fxTWA-PLS and the calibration of modern charcoal against GLM modelling of modern burnt area. We have applied this approach to reconstruct changes in burnt area through the Holocene for the Iberian Peninsula. The good predictive power of the fxTWA-PLS-derived fire–vegetation relationship and the plausibility of the palaeofire reconstructions with respect to known climate changes in the region suggest that this calibration approach could be applied more generally to provide quantitative reconstructions of past fire regimes in other regions where there are limited charcoal data, and pollen data are more abundant.