Variations in biomass burning across the boreal biome (a transcontinental circumpolar band located between 50 and 70^oN latitude) were analyzed based on 88 sedimentary charcoal records from lakes and peat bogs characterized by robust dating control. Even though some geographic areas are less well represented than others, a broad view of boreal forests’ fire activity during the investigated period is provided (Figure 1). Charcoal data were selected from the latest version of the Global Charcoal Database (GCD v4.0.6) compiled by the Global Palaeofire Working Group (GPWG¹) or from a previous published synthesis (Molinari et al., 2020). Metadata including detailed information about each record are presented in Table 1.

To reconstruct regional and continental composite curves of changes in fire activity over the last century, all charcoal data were first converted to charcoal accumulation rates (CHARs, i.e., the number, area or weight of particles cm-2 year-1). A standardization procedure was then applied in order to allow comparison within and between charcoal records obtained from various depositional environments and quantified with different laboratory techniques (Power et al., 2008). This technique includes: (1) a min-max rescaling of CHAR values, (2) a Box-Cox transformation to homogenize the within-record variance, and (3) a Z-score conversion using a base period between 1900 and 2000 CE. Palaeofire analyses generally omitted this time interval from the base period because of the large impacts of anthropogenic activities on ignitions and suppressions during the 19^th and the 20^th century (Power et al., 2008), which, instead, is what we want to capture in this study. Transformed charcoal data were then bootstrap re-sampled 1000 times with a moving window procedure using non-overlapping bins of 1-yr. Re-sampled time series were then smoothed using a LOcally WEighted regression Scatterplot Smoother (LOWESS) with a window half-width of 10 years. Charcoal composite series mean and 97.5% confidence intervals (hereafter CI) were calculated by averaging the smoothed and bootstrapped data series (Daniau et al., 2012). For each sub-area taken into consideration, composite curves were finally produced with the method implemented in the R (R Core Team, 2016) package “paleofire” version 1.2.4 (Blarquez et al., 2014).

We used the LPJ-GUESS-SIMFIRE-BLAZE model to perform a transient simulation over the last century of vegetation and fire dynamics as influenced by changes in climate, land use and population density. Here, the empirical fire model SIMFIRE (SIMple FIREmodel, Knorr et al., 2014) produces an annual burned area for each grid-cell (0.5^o x 0.5^o) containing a charcoal record based on fire weather, fuel continuity, fire-biome and human population density. Using a monthly fire climatology, daily values of burned area (as fraction of grid-cell) are computed and used as a fire probability for the combustion model BLAZE (Blaze-Induced Land–Atmosphere Flux Estimator, Rabin et al., 2017). In LPJ-GUESS (Smith et al., 2001, 2014) natural vegetation is represented by a multitude of patches with different “time since last disturbance” to account for cohorts of different age and height. Each patch is tested against the aforementioned fire probability on a daily basis. In case of a successful ignition, the whole patch is affected by fire and BLAZE simulates fire driven fluxes within the biosphere (live to litter) and from the biosphere to the atmosphere. The combustion model computes survival-probabilities for stands of trees based on vegetation properties (like height or diameter-at-breast-height, DBH) and potential fire-line intensity [based on a revised set of FullCAM parameters following Surawski et al. (2012)] derived from available fuels and local meteorology. In the event of a fire, the stochastic survival of each tree is tested individually. Total carbon fluxes (the total amount of C emitted by fire to the atmosphere) are then computed from the amount of tree-deaths and as fractions of above-ground litter combusted based on the fire-line intensity. SIMFIRE-BLAZE incorporates the HYDE 3.1 population density dataset (Klein Goldewijk et al., 2011) within LPJ-GUESS.

Patterns and annual cycles of burned area and carbon emissions within the range of satellite-based observations have been previously demonstrated by Kloster et al. (2004) to be reasonably simulated by the Arora and Boer (2005) fire model implemented in the Community Land Model (CLM-CM).

For each grid-cell containing a charcoal record, the probability to burn per year (fraction year^–1, hereafter fire probability Fp) and the total carbon flux emitted to the atmosphere per year (kg C m^–2 year^–1, hereafter carbon emitted Ce) simulated by LPJ-GUESS-SIMFIRE-BLAZE during the period 1900–2000 CE were treated analogously to sedimentary charcoal data in order to estimate regional and continental fire activity. To reduce the high year-to-year variability, a 10-yr running mean filter was applied before the Z-scores were derived. This same data processing minimized inconsistences in the statistical analysis and maximize comparability between the different datasets. In this way, discrepancies in the data-model comparison can only be explained by model deficiencies.

In order to evaluate the agreement between reconstructed and modeled fire dynamics in boreal forests since the beginning of the 20^th century, we calculated the Spearman’s rank correlation coefficient (Spearman, 1904) between Z-scores of transformed charcoal influxes (hereafter tCHAR Z-scores) and Z-scores of transformed model output of fire probability (hereafter tFp Z-scores) or carbon emitted (hereafter tCe Z-scores) for the boreal forests as well as for each region and continent. The Spearman correlation measures the strength of the monotonic relation between two variables, is robust to the presence of outliers, especially when there is sufficient sample size, and is applicable for any continuous or ordinal scale (Sheskin, 2007). The statistical analyses were run for the whole time period taken into consideration (1900–2000 CE), as well as for the two different sub-periods (1900–1950 CE and 1950–2000 CE).

Reconstructions of fire activity across the boreal forests based on sedimentary charcoal records during the period 1900–2000 CE show different trajectories (Figure 2 and Table 2).

In Alaska (Figure 2A), despite a downturn between 1910 and 1920 CE, a general increasing trend is recorded until 1980 CE. Then tCHAR Z-scores progressively decrease, with only a slight increase between 1998 and 2000 CE. In Canada (Figure 2B), charcoal values show a decrease during the period 1900–1955 CE, while two peaks (followed by downturns) are registered around 1965 and 1985 CE. An increase is recorded between 1998 and 2000 CE. The charcoal composite curves for North America (Figure 2F) and for the whole boreal forests (Figure 2H) record values around the long-term mean between 1900 and 1990 CE. A decreasing trend characterize the following years, with a little increase between 1998 and 2000 CE.

Despite large uncertainties, tCHAR Z-scores for Norway (Figure 2C) show a decreasing trend during the first half of the 20^th century, increasing values between 1950 and 1985 CE and a decrease at the end of the century. Sweden (Figure 2D) and Fennoscandia (Figure 2G) document decreasing charcoal trends between 1900 and 1930 CE. Increasing values are recorded during the period 1930–1960 CE, followed by a decrease with a minimum at 2000 CE. After increasing values at the beginning of the century, Finland/Russia Karelia (Figure 2E) registers a general decreasing trend between 1905 and ca 1930 CE. Increasing charcoal values are recorded until 1990 CE, followed by a decrease.

Changes in tFp and in tCe Z-scores simulated by LPJ-GUESS-SIMFIRE-BLAZE during the period 1900–2000 CE show a general progressive rise toward the end of the 20^th century for all the areas taken into consideration (Figure 2 and Table 2). Large oscillations occur as a direct consequence of climatic drivers.

Observed and modeled boreal forest fire activity since the beginning of the 20^th century are compared by investigating trends in tCHAR, tFp and tCe Z-scores transformed time series (Table 2). For the whole time period taken into consideration (1900–2000 CE) as well as for the two different sub-periods (1900–1950 CE, 1950–2000 CE), the Spearman rank correlation coefficient (Spearman, 1904) is also calculated.

In Alaska, based on 9 sites, an increasing trend in tCHAR (not-significant), tFp and tCe Z-scores is recorded when the whole 20^th century is considered and, only in tCHAR and tFp, during the time interval 1900-1950 CE. These results are confirmed by the statistical analyses, indicating positive correlations (only significant among tCHAR and tFp, Table 3) for these time periods. Instead, during the second half of the 20^th century the correlation between observations (decreasing) and simulations (increasing) is negative (Tables 2, 3).

For what concerns North America, on a total of 49 sites, an agreement between the different trends (all increasing but significant only for tFp) is recorded during the time interval 1900–1950 CE (Table 2). The statistical analyses confirm these results, showing a positive correlation (only significant among tCHAR and tFp) for charcoal data and model outputs during this period (Table 3). On the contrary, observations decrease and simulations increase between 1900 and 2000 CE and for the second half of the 20^th century (Table 2), with a negative correlation between the variables (Table 3).

In Canada (40 sites), Sweden (22 sites), Fennoscandia (39 sites) and for the whole boreal forests (88 sites) observations (increasing, but non-significantly for Canada between 1950 and 2000 CE and for the whole boreal forests between 1900 and 1950 CE) and simulations (decreasing) follow opposite trends during the 20^th century as well as for the two distinct periods (Table 2). Despite that, in Canada a positive non-significant correlation exists among tCHAR and tFp and tCe during the first half of the 20^th century (Table 3), while in the whole boreal forests a positive non-significant correlation is registered between tCHAR and tFp for the period 1900–1950 CE.

In Norway, based on 9 sites, tCHAR, tFp, and tCe record an increasing trend during the second half of the 20^th century. These results are confirmed by the statistical results, showing a positive relationship between the variables (only significant between tCHAR and tFp, Table 3). Instead, during the two different time intervals the correlations among observations (decreasing) and simulations (increasing) are negative (Tables 2, 3).

For Finland/Russia Karelia, based on 8 sites, observed and modeled data show an analogous (increasing, non-significantly for tCHAR between 1900 and 2000 CE) trend when the whole 20^th century and during the second half of the 20^th century (Table 2). This is in agreement with the Spearman rank correlation analyses indicating positive correlations (non-significant for tCHAR-tCe between 1900 and 2000 CE) among the variables during these periods (Table 3).

Granted all simplifications, generalizations and uncertainties in LPJ-GUESS-SIMFIRE-BLAZE structure, processes and input data (Hamilton et al., 2018), the use of sedimentary charcoal data to evaluate fire model simulations requires a careful consideration of the uncertainties connected to palaeofire data. In fact, even if time series of sedimentary charcoal abundance have become a widely used tool for reconstructing long-term trends in biomass burning (Whitlock and Bartlein, 2004), interpretation rests on understanding the source area of charcoal particles (i.e., the area burned and the distance to the site) and the ability to accurately infer other properties of the fire regime, such as frequency, intensity, and fuel type (Lynch et al., 2004). Constraining the charcoal source area is not a trivial task because charcoal abundance is influenced by vegetation, fire climate, and fire weather as well as site-specific processes related to charcoal deposition and burial (Clark, 1998). Despite advances in model-data integration (e.g. Higuera et al., 2007, 2011; Kelly et al., 2013) have helped validate charcoal-based fire reconstructions and identify uncertainties associated with natural (e.g., charcoal production, transport, and deposition) and analytical processes, reconstructing the area burned thus remains a challenge (Hennebelle et al., 2020). Furthermore, variability in the charcoal record comes from a variety of sources, including temporal variations in primary charcoal production, transport, and deposition; bioturbation and addition of secondary (remobilized) charcoal; differences in sampling, chemical treatment, and charcoal identification and quantification that can obscure the relationship between charcoal accumulation rates and fire (Iglesias et al., 2015). More specifically for our study, the uppermost sediments of the charcoal records might be less tightly constrained due to: (1) the difficulty to fully recover sediment/water interface during coring, (2) the high water content which tend to lower charcoal concentrations, and (3) the potential impact of human activities on catchment processes and sedimentation (Gale, 2009). Despite the recent development of Bayesian age-models (Blaauw et al., 2020), uncertainties in charcoal records can also be associated with a poor chronological control for recent times (where most of the dating are based on radiocarbon results) and to the different processes of charcoal production, transportation and deposition in the catchments, partially solved by field experiments that enable a better understanding of charcoal accumulation through time (Higuera et al., 2007; Hennebelle et al., 2020). Although charcoal is a proxy of carbon released by biomass burning, sedimentary-based reconstructions are currently insufficient for providing quantitative estimates of carbon released by fires because the links between stratigraphic charcoal and gaseous carbon fluxes are difficult to assess (Carcaillet et al., 2002). In boreal forests carbon release comes primarily from the ground layer (woody debris, humus, peat), which is more susceptible to becoming completely consumed than is the aboveground biomass (Kasischke et al., 2000). Finally, our analyses are based on charcoal records primarily collected from Canada (40 sites) and Sweden (22 sites), and this may have partly biased the results because some trends might be less robust than others.

The decline in fire activity over the 20^th century shown by our charcoal-based reconstructions for most of the boreal regions taken into consideration reiterate the results of previous palaeoecological studies at the global scale (Marlon et al., 2009a), in Australia (Mooney et al., 2011), in Europe (Molinari et al., 2013) and in the Americas (Power et al., 2013). A number of causes have been suggested for this downturn, including: (1) land-use change, resulting in landscape fragmentation and a generally less flammable fuel (Pyne and World, 1995), (2) active fire suppression management policies (Granström and Niklasson, 2008; Marlon et al., 2009a), and (3) an economic and cultural transition from traditional fire practices (such as slash-and-burn cultivation) to modern agriculture and forestry activities, especially in Fennoscandia (Wallenius, 2011). On the contrary, the increase in fire probability and, thus, in carbon emitted simulated by LPJ-GUESS-SIMFIRE-BLAZE are the result of the factors which determine the probability of a grid-cell to burn per year, such as: (1) the shift of many grid-cells taken into consideration from tundra to boreal fire-biome with a different flammability due to a change in CO₂-fertilization and climate (rise in both total annual precipitation and temperature), leading to an increase in the maximum annual fAPAR (fraction of Absorbed Photosynthetically Active Radiation, a proxy for the amount of vegetation growing in a grid-cell), and (2) a rise in the maximum annual Nesterov index, which is a proxy for fire-weather conditions and the curing of fuels. Both these factors lead to higher amounts of fuels and more fire-promoting weather conditions with a consequent increase of fire occurrence. So, when looking only into climatic drivers, an increase in boreal forest fires is expected under present and future climate change (Flannigan et al., 2005, 2009a; Ali et al., 2009). This is probably the reason of the agreement between trends in sedimentary charcoal data and model outputs during the 20^th century in Alaska, a region where fire activity was mainly driven by warm temperatures and vegetation flammability (Hoecker and Higuera, 2019). The good correlation between observations and simulations detected in Finland/Russia Karelia is instead probably due to reworking of charcoal particles in the uppermost levels of the sequences due to ditching and clearcutting in the forests (Pitkänen and Huttunen, 1999).

The influence of anthropogenic factors (e.g., fire suppression, forest management and ignition control) on boreal fires dynamics is likely to be the main cause of the data-model disagreement. In Canada, the decrease in charcoal-based reconstructions during the first half of the 20^th century is mainly the result of forest management (Bergeron et al., 2001). For the same period, a gradual abandonment of slash-and-burn cultivations in Norway and Finland/Russia Karelia (Lehtonen and Kolström, 2000; Rolstad et al., 2017), and of an early use of fire management practices in Sweden (Zackrisson, 1977; Granström and Niklasson, 2008) resulted in a similar decrease in these regions. The good correspondence between fire dynamics observations and simulations in Alaska and in North America is instead connected to a rise in wildfire extent due to an earlier impact of climate warming in these northern latitude areas (Chapin et al., 2008; Kelly et al., 2016). Regarding the second half of the 20^th century, the negative correlation between sedimentary charcoal trends and model projections recorded in Alaska and Canada is primarily related to the high efficiency of active fire suppression policies thanks to increased landscape fragmentation (Weir et al., 2000; Kasischke et al., 2010). For Canada this is in contrast with the data included in the Large Fire Database (LFDB, which includes information on fires larger than 200 ha in area) indicating, as suggested by LPJ-GUESS-SIMFIRE-BLAZE outputs, a general increase in area burned during the period 1959–1997 CE (Stocks et al., 2003). Discrepancies between charcoal-based reconstructions and simulations are also registered in Sweden, where the decline in forest fires is linked to both active fire suppression (Niklasson and Granstrom, 2000) and climate variations [particularly water availability as suggested by Drobyshev et al. (2012)]. Norway and Finland/Russia Karelia are the only sub-areas where both observations and model projections register an increase probably connected to a comeback of prescribed burning practices for biodiversity goals (Bleken et al., 1997; Lindberg et al., 2000) or to sediment remixing at the top of the sequences due to forest management (Pitkänen and Huttunen, 1999).

As in most fire models (Hantson et al., 2016), the only anthropogenic driver mediating fire probability (and thus carbon emitted) in LPJ-GUESS-SIMFIRE-BLAZE is human population density which, in the high latitudes investigated by the present study, is at very low levels and its general mild increase toward the end of the century is outcompeted by the climatic drivers which promote a rise in fire probability. However, even under these low population densities, people can have large-scaled impacts on fire occurrence as indicated in the previous subsection. Detailed representation of forest management (and the associated stand structure) and landscape fragmentation (and their absence in the past) should be taken into account in high-latitudes fire-simulations for a better representation of the human impact on fire regime. Furthermore, a more detailed simulation of fires in the boreal vegetation, especially peatland and shrubs, with a special focus on soils, is needed to more accurately represent fires in these regions. The same goes for the representation of fire and forest management (e.g., fire prevention, prescription). While, in fact, the impact of such factors on the global scale is probably small, it is clear that in practise these can override the predominant climate signal currently simulated by LPJ-GUESS-SIMFIRE-BLAZE on a regional scale.

This study further indicates the need for additional sampling and high-resolution analyses of charcoal records from areas inadequately represented by the present dataset and covering the more recent years (i.e., the beginning of the 21^st century). Due to the low number of charcoal records available for the last 20 years, it was in fact impossible to perform the statistical analyses beyond the end of the 20^th century. Despite that, the rise in both charcoal-based reconstructions and model outputs at the end of the study period in Alaska, Canada, North America and for the whole boreal forests (Figures 2A,B,F,H) make us assuming increasing trends until the present time as already suggested by Turetsky et al. (2011), Sanford et al. (2015), and Tymstra et al. (2020). On the contrary, the decline in charcoal-based reconstructions register in Fennoscandia (Figures 2C–E,G), probably connected to the effectiveness of fire suppression policies in the modern forests to support a growing timber industry that makes large fire events rare in these regions (Pinto et al., 2020), is in contrast with an increasing fire risk in these areas raised by the outbreak of exceptionally large fires in Sweden in recent years, and especially during summers 2014 and 2018 (San-Miguel-Ayanz et al., 2019). A better understanding of the spatio-temporal variability of fire activity is thus needed in order to maximize the possibility of data-model comparison and better evaluate the ability of ecosystem models to more realistically simulate future biomass burning.