Interactive comment on “ Comparing modelled fire dynamics with charcoal records for the Holocene ”

General comments:  1. This manuscript stands to benefit from a thorough language edit. However, these issues could easily be resolved.  The revised manuscript will include language edits (see also response to reviewer #2).   2. I agree with Reviewer #2 that the recently published paper by Molinari and coauthors (2013) should be drawn into the manuscript. Further, a recently published fire data-­‐model paper by Feurdean and co-­‐authors (2013) should also be compared against the results presented here.  Feurdean, A., Liakka, J., Vannière, B., Marinova, E., Hutchinson, S. M., Mosburgger, V., & Hickler, T. (2013). 12,000-­‐Years of fire regime drivers in the lowlands of Transylvania (Central-­‐Eastern Europe): a data-­‐model approach. Quaternary Science Reviews, 81,48-­‐61.  Thanks for pointing us to that paper. In the revised manuscript we include it in the introduction to inform about local studies dealing with data-­‐model comparison. We did, however, not include a comparison of our results against the one published in Feurdean et al. (2013). Our study compares large area means of fire activity over the Holocene. Feurdean et al. focuses on a small area, which is out of the scope of our study.     3. Aside from a few minor grammatical issues, I appreciate that the authors state their research questions in the Introductory of the manuscript. I encourage the authors to revisit (and answer) their questions in the Conclusions as well.


Introduction
The current interglacial period, the Holocene, that started about 11 650 cal yr BP (Masson-Delmotte et al., 2013), is characterized by a relatively stable climate.The main climatic and environmental changes resulted primarily from slow changes in orbital forcings.The maximum of incoming summer solar irradiance in the Northern Hemisphere was about 11 000 cal yr BP, when it was higher by ca.50 W m −2 (65 • N: 528.45 W m −2 at 11 000 cal yr BP; Berger, 1978;Berger et al., 1998;Tzedakis et al., 2012).However, until about 9000 cal yr BP, the climate in the northern temperate and high latitude regions was strongly affected by the remains of the northern hemispheric ice sheets.Based on terrestrial proxy records, a time slice of 6000 cal yr BP was chosen as a data-model comparison reference period for the mid-Holocene by the Paleoclimate Modelling Intercomparison Project (PMIP, Braconnot et al., 2007).
Because ample records of climate and environmental changes are available for the Holocene (e.g.Wanner et al., 2008), this period is a well-suited one to test climate and biospheric models (e.g.Fang et al., 2013) and compare these results with the syntheses of geological archives (e.g.BIOME6000).Simulated and reconstructed changes in climate and vegetation cover (e.g.Claussen, 1997;Kutzbach and Liu, 1997) have often been compared to the 6000 cal yr BP time slice.More recently, the palaeomodelling research has focused on simulating transient changes, for example, in the sea surface temperature (Lorenz et al., 2006), sea ice (Fischer and Jungclaus, 2010), land surface climate (Renssen et al., 2004), and comparison of T. Brücher et al.: Comparing modelled fire dynamics with charcoal records for the Holocene available pollen records with the tree cover changes (Kleinen et al., 2011).Progress in the synthesis of the Holocene land cover and land use changes (e.g.Gaillard et al., 2010) and climate proxy records (Marcott et al., 2013) now provides a basis for detailed model-data comparison throughout the Holocene period.
Fire is an important process that affects climate through changes in CO 2 emissions, albedo, and aerosols (Ward et al., 2012).In addition, it disturbs vegetation distribution (Sitch et al., 2003).Fire-history reconstructions from charcoal accumulations in sediment have indicated that biomass burning has increased since the last glacial maximum (Power et al., 2008;Marlon et al., 2013).Recent comparisons with transient climate model output suggest that this increase in global fire activity is linked primarily to variations in temperature and secondarily to variations in precipitation (Daniau et al., 2012).A new aspect of the recent generation of Earth System Models (ESM) is the implementation of fire models (Arora and Boer, 2005;Kloster et al., 2010;Thonicke et al., 2010;Pfeiffer et al., 2013) that allow testing hypotheses generated through reconstructions of palaeofire data on the controls of fire across a range of spatial and temporal scales.Fire model outputs have included simulated burned areas and CO 2 emissions that can be evaluated against presentday observations from remote-sensing products.For example, Kloster et al. (2010) demonstrated that the Arora and Boer (2005) fire model implemented in the Community Land Model (CLM-CN) reproduces reasonable patterns and annual cycles of burned area and carbon emissions within the range of satellite-based observations.Simulated mean burned area (327 Mha yr −1 ) was at the lower band of satellite observations (329 to 401 Mha yr −1 ), and modelled carbon emissions varied between 2.0 and 2.4 Pg C yr −1 , also close to the low end of satellite-based products (2.3 and 2.7 Pg C yr −1 ).
Recent progress in the analysis and syntheses of charcoal records led to a Global Charcoal Database with hundreds of records from around the world (Power et al., 2008;Daniau et al., 2012).To analyse the temporal changes in the charcoal database, individual charcoal records are transformed and standardized to allow comparisons across the multiple types of records (e.g.Power et al., 2008;Marlon et al., 2009;Daniau et al., 2012;Marlon et al., 2013).The standardization allows analysis of trends in different regions using charcoal records obtained from a wide range of depositional environments and quantified with different laboratory techniques.
Some recent studies compare fire model results on a local (Feurdean et al., 2013) or regional scale (Molinari et al., 2013) to a composite or single site measurements of charcoal.In this study we use the advantages offered from a mechanistic global model and ask "What factors are causing the variations in fire?" and analyse how well fire models reproduce reconstructed trends in fire activity during the last 8000 yr.Second, we ask "How well can global fire models reproduce reconstructed fire activity?" and from a methodological perspective, "What is the best way to compare global fire model output with charcoal records?".The last question is important because (i) the fire model provides quantitative information about burned area and fire-related emissions of CO 2 , but charcoal-based palaeofire data only provide information about relative changes in biomass burning for specific regions and time periods; and (ii) since the charcoal records are interpreted via a non-linear power transformation, the model can be used to relate trends in charcoal to trends in quantitative changes in burned area or fire carbon emissions; for example, does a 2 × increase in the standardized charcoal reconstruction reflects a 2 ×, 200 ×, or some other increase in area burned?The essential steps in answering these questions are to test how well the model can reproduce the reconstructions, develop metrics for characterizing the comparisons, and to understand how the data transformation and standardization affects the model output.
In the following, Sect. 2 describes the data and model used in this study in detail, including the transformation and standardization of model output to Z-scores.Section 3 discusses the drivers of variation, comparison of model results and charcoal data on (i) large (continental) and (ii) regional (sub-continental) spatial scales followed by a summary of key findings in Sect. 4.

Charcoal reconstructions
The Global Charcoal Database (GCD version 2.5) is used to determine regional and global palaeofire trends from 218 sedimentary charcoal records covering part or all of the last 8000 yr.To retrieve regional and global composites of changes in fire activity over the Holocene, charcoal accumulation in sediments is compiled and transformed (Power et al., 2008;Marlon et al., 2009) as a standardized measure frequently used by the palaeofire community to compare aggregated values of past fire activity.The transformation aims to homogenize the variance within individual time series using a Box-Cox transformation, and then rescales the values and calculating anomalies to identify time periods of lower-and higher-than-modern fire activity.Modern is usually defined as a time window prior to the industrial period.The Box-Cox transformation is necessary, as reconstructions of charcoal records are based on a wide range of processing techniques and various types of sites and differ by several orders of magnitude.The Z-score transformation is based on three major steps: (1) rescaling values using a minimax transformation; (2) homogenization of variance using the Box-Cox transformation; and (3) rescaling values once more to Z-scores.The Z-score transformation is applied to each time series (ci s ) of charcoal influxes (ci) for each site (s) separately.First, the linear minimax transformation ci s (Eq. 1) is used to scale the values ci s between 0 and 1, by subtracting the minimum value min(ci s ) and divide by the amplitude (max(ci s ) − min(ci s )) of the time series.ci s = ci s − min (ci s ) max (ci s ) − min (ci s ) (1) The non-linear Box-Cox transformation ci * s (Eq.2) reduces high numbers and removes outliers to achieve a more Gaussian-like distribution (e.g.Fig. 3 in Power et al., 2008): (2) Here, the parameter λ s is estimated by a maximumlikelihood method (Venables and Ripley, 1994) for each time series at each charcoal site (s) separately.To avoid a division by zero, a small constant α is added (here: α = 0.01).Afterwards, the minimax and Box-Cox transformed time series ci * s (Eq.2) is normalized by subtracting the mean value of a predefined base period ci * s and dividing the anomalies by the standard deviation S * ci s of the minimax and Box-Cox transformed time series ci * s (Eq.2): To retrieve regional, aggregated values Z region out of the site information ci Z s (Eq.3), all time series ci Z s are linearly averaged (Eq.4) by deviding with the number of sites (N sites ): Even though the Z-score transformation is not linear it is still rank conserving.In order to calculate area composites, each record was sampled (without interpolation) at 20 yr intervals and afterwards smoothed (LOcally-WEighted regression Scatter plot Smoother -LOWESS) by running a 250 yr moving window.Furthermore, a bootstrap analysis (sampling by site) was used to obtain the 95 % confidence intervals for the Z-score of a region.

CLIMBA -a fast global carbon cycle model
The computational efficiency of ESMs is rather low, therefore it is a challenge to do interactive simulations over the Holocene, as, for example, done by Fischer and Jungclaus (2010) et al., 2000;Ganopolski et al., 2001) and JSBACH (Reick et al., 2013;Schneck et al., 2013;Brovkin et al., 2009;Raddatz et al., 2007), which is the land component of the Max Planck Institute Earth System Model (MPI-ESM; Giorgetta et al., 2013).While CLIMBER-2 simulates the atmosphere and land processes at roughly 51 Because the base climate is chosen from the pre-industrial simulation, the yearto-year variability is given by the variability of the base climate alone.The climate forcing includes values for temperature, precipitation, radiation balance, and atmospheric CO 2 concentration.Given this climate and orbital forcing, JS-BACH runs offline for one model-year.The simulated change in the total land carbon feeds back to CLIMBER-2 as a carbon flux to the atmosphere (negative or positive) and closes the carbon cycle.This coupling scheme does not allow us to account for biogeophysical feedbacks linked to changes in burned area, as the CLIMBER and JSBACH models are only coupled via the carbon cycle.The biogeophysical effects from changes in vegetation distributions are solely simulated by the land surface model of CLIMBER-2 (Brovkin et al., 2002).JSBACH includes a dynamical vegetation scheme and calculates disturbance of vegetation by natural fire occurrence and windthrow.While the default JSBACH version (Reick et al., 2013) uses a simple fire module (Brovkin et al., 2009), the fire algorithm used in this study is based on Arora and Boer (2005), which simulates fire occurrence as a function of fuel availability, fuel moisture, and ignition source (Kloster et al., 2010;Krause et al., 2014).The fire model calculates a potential burned area, which is simulated as a function of fuel moisture and wind speed.Fuel moisture is not explicitly simulated in the model.Soil moisture is taken as a proxy for fuel moisture.The associated carbon emissions are calculated as a function of area burned and available fuel load.There is no distinction as to whether carbon is emitted as CH 4 or CO 2 .Kloster et al. (2010) shows, that the Arora and Boer (2005) algorithm implemented in the Community Land Model (CLM-CN) reproduces reasonable patterns and annual cycles of burned area and carbon emissions within the range of satellite-based observations.
In this study we keep lightning constant while recognizing it is a weather-driven process and therefore variable within T. Brücher et al.: Comparing modelled fire dynamics with charcoal records for the Holocene the changing climate of the Holocene.We also do not account for human ignitions and anthropogenic land use or anthropogenic land cover changes and acknowledge these factors can impact the carbon cycle.Therefore, the research focus is designed to explore potential vegetation and natural fire occurrence.

Experimental set-up
A transient simulation spanning 8000 yr has been performed starting at the mid-Holocene at 8000 cal yr BP (8 ka = 8000 cal yr BP; 0 ka = 1950 AD).The initial conditions are taken from a control simulation with CLIMBA at 8000 cal yr BP with an atmospheric CO 2 concentration of 260 ppm, in accordance with the reconstruction from ice core (Elsig et al., 2009).The transient simulation includes an orbital forcing after Berger (1978), fixed greenhouse gases and aerosol concentration, and ignores changes in sea level and land ice.At the end of the transient simulation, atmospheric CO 2 concentration is simulated as 272 ppm, which is lower than observed by ca. 8 ppm (Monnin et al., 2004).As the focus of this study is on natural vegetation and natural fire occurrence, anthropogenic land use emissions (e.g.Pongratz et al., 2009;Ruddiman, 2003) are neglected.If such emissions were included, they would lead to higher atmospheric CO 2 concentrations at pre-industrial (PI) times, for example, about 18 ppm, if a land use emission scenario based on HYDE (Goldewijk, 2001) would be used.
Reconstructions and model data are restricted to between 8000 and 200 cal yr BP to exclude the onset of industrialization and the large human impacts on fire activity during the subsequent centuries.Furthermore, within this period the charcoal database has the greatest number and sample density of palaeofire records and therefore the highest data quality.

Processing of model output for the comparison with charcoal data
For the comparison with reconstructions of palaeofire, two output variables of the fire model at each grid box (g) have been used (1) the fraction of grid box area burned per year f g [yr −1 ] and (2) the total carbon flux to the atmosphere To compare aggregated model results of burned area (F ) and carbon emissions (C) with regional estimates of fire activity out of the charcoal database reported as Z-scores (Z), the model output is processed using two different approaches: (i) at time t the grid box values f g (t) and c g (t) related to the region under investigation are weighted by its area a g [m 2 ] and summed up to get accumulated, regional numbers (Eqs.5 and 6): (ii) Furthermore, Z-score transformed (Eq. 3) time series f Z g and c Z g are derived from f g and c g of each grid box.Then they are linearly averaged to achieve regional time series of burned area F Z region and carbon emissions C Z region (Eqs.7 and 8).Without using an area weighting function the local information of vegetated area get lost by just dividing with the number of grid boxes per region N g : N g (7) These two different approaches of absolute values (Eqs.5 and 6) and regional averaged Z-scores (Eqs.7 and 8) of model output are used in the following.
To reduce the high year-to-year variability, a 250 yr running mean filter is applied before the Z-scores are derived.For reconstructions and model data, the used Box-Cox transformation and the normalization afterwards are based on the full period (7800 yr), in particular the same period is used to calculate the mean and standard deviation used in Eq. ( 3).Hence, charcoal influxes and model output are treated in the same way to minimize inconsistencies in the statistical analysis and maximize comparability.Therefore, regional disparities in the model-data comparison cannot be explained by differences in data processing.More specifically, only model deficiencies can explain data-model differences.

Changes in fire activity at 8000 cal yr BP
The natural variability and trends in fire occurrence simulated by the model are driven by changes in climate and climate induced changes in vegetation cover.We discuss in the following the simulated changes in fire occurrence in conjunction with changes in precipitation and temperature as the dominant drivers for vegetation and fire activity.Furthermore, we compare modelled and observed palaeofires and discuss the advantages of transforming model results into aggregated and Z-score-transformed time series.
During the mid-Holocene, the Northern Hemisphere received more solar irradiation during the summer season relative to pre-industrial (Berger, 1978).This led to substantial summer warming of up to 2-3 • C which was most pronounced in high northern latitudes (Renssen et al., 2009).Northern subtropics, including North Africa, were substantially wetter, presumably due to an intensification of the monsoon circulation.This led to the significant increase of vegetation cover in subtropical drylands and in the Sahel/Sahara region (Prentice et al., 1992;Claussen, 1997).While the Northern Hemisphere was warmer over the Holocene, the temperature anomalies in southern extra-tropics (30-60 • S) were small (Wanner et al., 2008).These general features of the mid-Holocene climate changes are well reproduced by the CLIMBER-2 model as also seen by previous studies (e.g.Claussen, 1997).Temperature anomalies simulated by the CLIMBER-2 model on a yearly mean basis are within the range of −0.5 to 0.5 • C except the area of western African and Indian monsoon with a strong dipole of a cooler, wetter region and a warm area (not shown).
Our simulated total burned area (Figs.1a and 3a) for the mid-Holocene is at 514 Mha yr −1 and increases slightly by 14 Mha yr −1 (ca.2.5 %) to 528 Mha yr −1 .The hotspots of burned area are located in tropical Africa, central North America, central South America, Australia and partly the Mediterranean region plus South Asia.The increase is mainly due to more fires in the Amazon, North America, South Asia, and the east coast of Australia (Fig. 1b).The main regions showing a decrease in fire activity (burned area fraction) are tropical West Africa, Central Australia, and Europe (Fig. 1b).Looking at carbon emissions (Fig. 1c), the pattern is similar to the pattern of burned area (Fig. 1a) capturing similar hot spots of palaeofire activity (note, the scale of the colour bars differ by a factor 1000).An increase of almost 5 % (0.29 Gt C yr −1 ) is calculated by the transient simulation.This trend in carbon emissions appears stronger than trends in burned areas.Although, there are some regions (e.g.central South Africa) showing a decrease in burned area along with an increase in carbon emissions.Compensating effects of declining burned area, but a higher vegetation fraction (Fig. 2b) along with more carbon in living biomass due to the atmospheric CO 2 fertilization could partly explain these opposing trends.Although the fertilization effect in JSBACH is high (Arora et al., 2013), it is probably insufficient to explain these trends.As a result, changes in the underlying vegetation may also play a significant role in driving palaeofire trends.Regions with an increasing higher burned area but decreasing carbon emissions -both linked to natural fire -are not simulated by JSBACH.
To understand the temporal evolution of regional fire occurrence Hovmöller diagrams of anomalies (8000 cal yr BP minus 200 cal yr BP) are shown (Fig. 2) for precipitation (Fig. 2a), desert fraction (Fig. 2b), carbon in living biomass (Fig. 2c), and burned area fraction (Fig. 2d; note, this is land data only).In general, the patterns of anomalies do not identify strong gradients with time, even before applying the 250 yr smoothing window.So, there are no abrupt changes with time simulated, however, the latitudinal gradient does identify some sharp boundaries, which are likely an effect of the coarse resolution of the 10 • latitudinal bands of CLIMBER-2 anomalies driving the land model JSBACH.
There are two dominant patterns of anomalies apparent: within the northern tropics the intensified and northward shifted monsoon system (Fig. 3a) during the Holocene leads to widespread greening (Fig. 3b) with an increased biomass (Fig. 3c) between 8000 and 4000 cal yr BP as previously observed by, for example, Prentice et al. (1992), Claussen (1997) and Brovkin et al. (2002).Because of the humid climate in the northern tropics, decreases in fire activity become  amplified with time.Therefore, we hypothesize the trend in fire activity is moisture driven and not determined by fuel.
The second prominent pattern lies around 20-30 • S, where zonal means of yearly precipitation point to drier conditions at 8000-6000 cal yr BP.For South Africa, all but the southern tip was drier during 8000 cal yr BP, as the monsoon system shifted northward (Fig. S1a in the Supplement).For South America a dipole of drier Amazonia and a wetter region south of 20 • S is simulated (Fig. S1b in the Supplement).Therefore, the small increase (up to 100 mm yr −1 ) in annual precipitation in Australia almost counterbalances on latitudinal means the decrease of precipitation in South America.Between 20 and 30 • S, the vegetated area increases over the whole period, while there is no prominent shift in the zonal averaged numbers of green biomass.The increase in vegetated area parallels the increase in fire activity.This trend is apparent because of the higher fire occurrence in South America and Africa, while a modelled dipole of changes in burned area over Australia amplifies the trend at ca. 20-30 • S and lessens the increase south of 30 • S.
For the northern extra tropics, the patterns are noisier (Figs. 1 and 2) and an overall prominent decline of the boreal forest is not as strong (Figs.2a and S1g in the Supplement) as expected after, for example, Claussen (1997) and Kleinen et al. (2010).The increased vegetated area at ca. 60 • N can mainly be linked to some greening spots in Asia and Alaska without impacting the fire activity (Figs.S1g and h in the Supplement).

Comparison on hemispheric scale
To compare the modelled and reconstructed numbers of aggregated Z-score values on a hemispheric scale, we investigate five spatial domains separately: global [90 • S-90 • N] (Fig. 3a), northern extra tropics [90-30 • N] (Fig. 3b), northern tropics [30-0 • N] (Fig. 3c), southern tropics [0-30 • S] (Fig. 3d), and southern extra tropics [30-90 • S] (Fig. 3e).These domains are chosen in analogue to Daniau et al. (2012).Shown are Z-score transformed time series of reconstructed (charcoal influx, Z) and modelled (burned area F Z and carbon emissions C Z ) fire activity.There is also the untransformed model output of burned area (F ) given, as well as the number of sites the reconstructions are based on.The untransformed figures for carbon emissions by fire (C) are not shown, as on regional to hemispheric scale F and C are almost linearly correlated (not shown here).Charcoal data are given as the median and the ±95th percentile of a bootstrap analysis (see Sect. 2.1).A large spread between median and the ±95th percentile highlights a low number of records in the database or a heterogeneous domain (e.g.Fig. 3c).
For three different time periods (8000 cal yr BP-PI, 8000-4000 cal yr BP, 4000 cal yr BP-PI) the rank correlation after Spearman (1908) is used as an objective measure to quantify the agreement between modelled and reconstructed values.The numbers are given for the correlations between (i) Z and F (ρ(Z, F )), (ii) Z and F Z (ρ(Z, F Z )), and (iii) Z and C Z (ρ(Z, C Z )) as a bar chart next to the global or large-scale zonal averages.
On global scale (Fig. 3a) an increase within the burned area is observed (red line = burned area).With a base of 218 sites, the spread between the ±95th percentiles and the median is rather low, suggesting that the global signal out of the reconstructions is rather robust.Modelled data do agree in the trend and magnitude (14 Mha, ca.2.5 %).If the same fire model is used within an earth system model, findings by Krause et al. (2014) suggest values in the same order of magnitude (ca.580 Mha burned area per year) by counting for present-day climate, pre-industrial land areas and no human ignition.However, if the fire model is used in a more realistic set-up for present-day climate (Kloster et al., 2010), a mean simulated burned area of 327 Mha yr −1 stays at the lower band of satellite observations from 329 Mha yr −1 (GFED) to 401 Mha yr −1 (L4JRC).Model results of carbon emissions vary between 2.0 to 2.4 Pg C yr −1 and are close to the numbers of satellite products (2.3 to 2.7 Pg C yr −1 ).While the rank correlation between the untransformed model data is small (ρ(Z, F ) = 0.18), the numbers are much higher (ρ(Z, F Z ) = 0.73 and ρ(Z, C Z ) = 0.77) if the Z-scores of charcoal values are compared to Z-scores of modelled burned area.If the full period of 7800 yr is divided into two subsets (8000 to 4000 cal yr BP and 4000 cal yr BP to PI), the correlation is reduced, especially for the latter one, which suggests the overall Holocene trend is reproduced, but not the sub-millennial variability.The decrease in fire can be related to drier conditions on a global scale.While changes in temperature are relatively small (ca.0.1 K, see Fig. S2a in the Supplement), decreases in yearly precipitation by ca.40 mm yr −1 and the small increase in biomass (i.e.fuel availability) dominate the effect on driving fire activity.The continuous increase in the carbon stocks is also supported by CO 2 fertilization (12 ppm increase; Fig. S2a in the Supplement; Keenan et al., 2013).
For the northern extra tropics, the charcoal data show a small increase, while the modelled fire activity stays almost constant (Fig. 3b) at ca. 164 Mha.A rank correlation gives a negative correlation coefficient for the Z-scores of burned area and carbon emissions with the charcoal reconstructions, while the correlation of burned area and the charcoal reconstructions is not significant (Fig. 3b).On the large area mean, the temperature decreases by 0.2 K, the climate gets drier (40 mm yr −1 ), and the biomass decreases.The shift toward reduced fuel availability and drier conditions seems to be compensated by drier conditions, which lead to almost no change in modelled burned area.
For the tropics (Fig. 3c and d) a strong increase in burned area is reported after 7000 cal yr BP in the charcoal database.The model results reflect this increase over the full period and all time series are positively correlated, ranging from ρ(Z, F Z ) = 0.42 to ρ(Z, C Z ) = 0.48 (for further details see Table S1 in the Supplement).While the numbers of Z-score transformed data suggest a strong change, the modelled numbers vary only by 6 Mha from 8000 cal yr BP to PI (northern and southern tropics), which suggests a change by roughly 4 %.In the case of untransformed modelled data, the correlation shrinks by factor of three, which supports the standardization of model output to improve comparability with charcoal influx values.In terms of large-scale averages, climate is not changing significantly over the 7800 yr in the southern tropics, while precipitation decreases in the northern tropics by ca. 10 %. (80 mm yr −1 ).As precipitation is the controlling parameter for tropical vegetation, biomass decreases in the northern tropics and increases slightly in the southern tropics.As both areas point toward an increase in fire activity, it seems that on a large-scale fire in the tropics is primarily determined by fuel abundance and moisture availability.
For the southern extra tropics (Fig. 3e), the level of reconstructed natural fire activity stays almost constant with a small decrease around 4000 cal yr BP.However, the model results show a strong increase in burned area and carbon emissions over the entire period of 7800 yr simulated.The rank correlation shows, that less than 15 % of variability is explained by the Z-score transformed data (ρ(Z, F Z) = 0.24; ρ(Z, C Z ) = 0.22) and the absolute, not transformed values of burned area do not correlate significantly (p < 0.05) with the reconstructions.An explanation for the disagreement could be the small land area in general and a large fraction of coastal area within this region (southern part of Africa, Patagonia, and partly Australia), which are in general difficult to represent in global climate models.Alternatively, the domain of southern extra tropics is the area with the highest simulated burned area (nearly 70 %).However, the simulated drop in temperature (by 0.2 K) is rather small, but the decrease in precipitation, leading to drier conditions, appears more significant.In the model the increase in Z-score values of burned area is likely linked to changes in precipitation.The modelled precipitation change does, however, not necessarily reflect observed changes, which require further investigations.
In general, the trend of simulated burned area and their carbon emissions point to a small increase in fire occurrence, which is reflected in the observed charcoal reconstructions as well.The regional correlation analysis explains a maximum of 25 % of the variance for the different areas.Running the correlation analysis for the two time segments (8000 to 4000 cal yr BP and 4000 cal BP to 200 cal yr BP) separately, the correlation decreases or even becomes negative, which suggests that the general trend of increasing burned area over the entire period (8000 cal yr BP until 200 cal yr BP) is partly reproduced and responsible for the correlation coefficient over the 7800 yr.
To test the robustness of these modelled trends, we created a four-member ensemble based on identical initial conditions.Due to the randomness inherent in the base climate (see Sect. 2.1) all simulations develop independently.The analysis of the ensemble shows that the simulated trends in fire occurrence are all similar, indicating that our results are robust (Fig. S2 in the Supplement).

Comparison on regional scale
To investigate the model performance on more homogeneous regions, the methodology used in the previous Sect.3.2.1 is applied on regional scale, following domains defined and discussed in Marlon et al. (2013) (Figs. 4 and 5).Most of the charcoal Z-score time series point to an increase in fire activity, except for the Central America tropics (Fig. 4c), where a u-shaped pattern of Z-score influxes comes across with a high centennial variability including a high uncertainty possibly caused by a low number of charcoal records available (n = 9).The Asia monsoon region (Fig. 4f) also has a limited number of charcoal records (n = 10).
On a first glance, JSBACH simulates an increase in burned area and carbon emissions in all regions except Europe (Fig. 4b), where burned area decreases between 6000 and 4000 cal yr BP by 25 % and remains afterwards at ca. 18 Mha (Fig. S4b in the Supplement).In addition, for all regions, the modelled timing of local minima or maxima does not match the reconstructions.The rank correlation on these time series indicates overall agreement that there was an increase in fire activity, with positive correlation coefficients between ρ = 0.32 and ρ = 0.66 (p < 0.05).The highest correlation is found for reconstructed charcoal influxes and modelled carbon emissions by wildfire in North America (Fig. 4a), which is also the region with the most charcoal data available.The charcoal reconstructions for Africa (Fig. 4d) fit the model results rather well (ρ = 0.46; p < 0.001), even though this domain is limited by a low number of charcoal sites.Repeating this correlation exercise on split time frames (8000 to 4000 cal yr BP and 4000 cal yr BP to PI), results in a decrease in explained variance, and partly negative correlations (Fig. 5).Only the African region (Figs.4d and 5) differs from this general pattern; as the correlation coefficient between the Z-scores of modelled carbon emissions and charcoal influxes increases over time; it has a lower value for the entire period (ρ 8k-PI = 0.46; p < 0.001) a higher value during the first half (ρ 8K-4K = 0.61; p < 0.001), but is negative for the second time interval.Another interesting aspect of the African palaeofire simulation is the slight decrease in burned area (and carbon emissions) with time over the entire period, although an increase is notable for the Z-score transformed values of carbon emissions and burned area.As a Z-score transformation removes the absolute value of a variable, averaging the Zscore time series of several grid cells can change the sign of the trend.As an example, two neighbouring grid boxes with one showing an increase by 50 % and the other showing a decrease by 50 % of fire activity can result in a Z-score of zero, as both trends cancel each other.If the boxes, however, differ in their burned area magnitude (with respect to square metre burned area), the absolute change would be different from 0.
Another unique result for Africa (Fig. 4d) is the interplay between burned area and carbon emissions.The correlation coefficient between Z-scores of modelled burned area and modelled carbon emissions is ρ = 0.98 (p < 0.001; Table S1 in the Supplement) at a minimum for all but one region (Figs. 3 and 4).Only the analysis of the African domain stands out with a lower value of ρ = 0.84 (p < 0.001).
By considering significant and positive rank correlations, we find for all regions -except for the Central America tropics (Fig. 4c) -a higher rank correlation between Zscore transformed data and reconstructions (ρ(Z, F Z ) and ρ(Z, C Z )), than for untransformed model data and reconstructions (ρ(Z, F )).Given the current study design and model configuration, the main driver of changes in Holocene climate (last 8000 yr) is orbital forcing; changes due to increasing CO 2 are of secondary order.Human impact in the experimental set-up is neglected.Furthermore, this study focuses on large, heterogeneous regions and differences in the results based on (Z, C Z ) vs. (Z, F Z ) are therefore minor.Looking at the ensemble results (Fig. S4 in the Supplement), it is apparent, that the trends simulated by each ensemble member are rather close to each other and differences between the model results and reconstructions are systematic.

Conclusions
Running a fast global carbon cycle model over the Holocene until pre-industrial times (here defined as 1750 AD), a transient climate change is simulated that closely tracks reconstructed patterns from proxy records.By developing the capability to simulate reconstructed trends in natural fire activity, we were able to compare simulations of paleofire activity at continental and regional scales using Z-score transformed charcoal influx data.
Over the past 8000 yr we simulate an increase of ca. 14 Mha (from 512 to 526 Mha) in burned area.The absolute numbers of burned area are high, however, considering that the human dimensions in terms of fire ignitions and fire suppression are neglected.In addition, the model accounts for fire activity given only potential natural vegetation (i.e. with no land use and land use cover change included).The increase itself and the variability on millennial timescales varies between and among regions.When the modelled time series are transformed to Z-scores and compared to reconstructed charcoal Z-scores, model and data agree for most of the regions.Thus, to address the question "How well can fire models reproduce fire activity?",we can state that our model simulates most of the trends in the reconstructions on millennial scales.Further investigations would be necessary to test whether the model performs well for the right reasons.
One limiting component of our model set-up is the choice of an EMIC as the climate driver.Especially on regional scales, our fire model results are limited by the quality of the climate forcing, which can (partly) explain our limited ability to reproduce centennial or millennial variations.Another regional model study by Molinari et al. (2013) compared charcoal records in Europe with an earth system model that simulated dynamic vegetation using LPJ-GUESS and two independent land-cover scenarios.The climate simulations and analysis of land cover change suggests biomass burning across Europe was primarily explained by vegetation, precipitation and temperature-related parameters during the early Holocene.Charcoal-based observations of increased fire activity during the mid-late Holocene were primarily driven by changes in anthropogenic land-cover, and secondarily by changes in vegetation and temperature (Molinari et al., 2013).As our study does not count for anthropogenic changes, we can conclude that we can conclude that the most parsimonious explanation suggest inaccuracies in our simulated climate is the most likely source of discrepancies between model and data.
For most of the investigated regions the model simulates an increase in burned area and carbon emissions.The trends in the carbon emissions were higher than trends detected in burned area.We propose several reasons for this observation: (i) CO 2 fertilization by the increasing concentration of atmospheric CO 2 increases the emissions per square metre burned area due to a higher carbon stock in the vegetation.(ii) The carbon stock of the fuel can increase with shifts in vegetation types due to changes in the climate, or fire occurrence may be altered because of climate changes.A rank correlation analysis points to the overall agreement between simulated and observed trends in fire activity over the whole period of study, while a rank correlation on time segments shows that the model does not match the centennialor millennial-scale variability.Agreement on variability on these timescales is not expected, as regional climate affects local fire activity, and there is no reason why the timing of the modelled climate variability should coincide with the climate of the past.The differences and similarities between reconstructions and model results are stable.The analysis of a fourmember ensemble shows rather small differences between the individual simulations and a consistent trend given by the ensemble members (Figs. S4 and S5 in the Supplement).
From a modelling perspective, this study helps to validate the capability of a model to simulate past fire activity.In addition, given that the fire model is not tuned by the charcoal data, the overall agreement (on hemispheric and regional scales) illustrates the ability to reconstruct broad changes in palaeofire activity from syntheses of palaeofire records in the Global Charcoal Database.Even regions which are sparsely covered by reconstructions correlate positively with the model results, pointing to the importance and increased confidence of using both data and models together to provide more complete spatial coverage of past fire activity.
Z-score transformed data do not provide information about quantitative changes in burned area, because the transformation is rank conserving but not linear.So, a given difference in Z-score values [ ] does not imply an increase or decrease by the same magnitude in burned area [Mha] among different Z-scores.This suggests regional averages of transformed and untransformed data may not necessarily result in the same trends.Thus, with respect to the research question "What is the best way to compare fire model output with charcoal records?", we conclude that it is more meaningful to convert the time series of modelled burned area or carbon emissions to Z-score for comparing modelled and observed palaeofire variability.While we do see some general agreement between model results and reconstructions, it is still unclear whether the absolute values of simulated burned area are capturing the right magnitude for past fire activity.While high fluctuations could suggest huge changes (e.g.Fig. 4e), the absolute change is rather small (2 Mha, ca.7 %).Future studies should consider methods of transforming model output variables and palaeo-proxy data consistently to increase the comparability of simulated and observed data.In this study the Z-score transformation helps to validate modelled natural fire occurrence and compare it to reconstructed values of charcoal influxes reported as Z-scores.

Fig. 1 .
Fig. 1.Yearly burned fraction of grid cell area [m 2 m −2 ] of natural fire activity (a) and carbon emissions [g C m −2 yr −1 ] (c) for the mid-Holocene (8 ka = 8000 cal yr BP) and their anomalies (b, d) to burned fraction with pre-industrial (PI = 200 cal yr BP) climate.

Fig. 3 .
Fig. 3. Averaged values for reconstructed and modelled biomass burning during the present interglacial as global values (a), for extra tropics (b, e), and tropics (c, d) separately.Reconstructions are shown by Z-scores of charcoal influxes (Z, pink), and Z-score transformed values of modelled burned area (F Z , black) and carbon emissions by fire (C Z , blue).Untransformed model output of burned area (F , red) and the number of sites used in the reconstructions (green) are also given.For all time series a running mean of 250 yr is applied.Please note the varying, relative scale of modelled burned area.The scale for the modelled Z-scores of burned area is determined by the maximum amplitude of Z-score transformed charcoal influxes (reconstructions).On the right side, the corresponding rank correlations ρ (after Spearman) are shown.Significant, positive values are given by filled bars for three different time windows: 8 ka-PI, 8-4 ka, and 4 ka-PI.

Fig. 4 .
Fig. 4. Same as Fig. 3, but for continental scale regions North America (a), Europe (b), Central America tropics (c), Sub Saharan Africa (d), Australian monsoon region (e), and Asia monsoon region (f).The definition of the domains is taken from Marlon et al. (2013).For the rank c or relations see Fig. 5.