A model for global biomass burning in preindustrial time : LPJ-LMfire ( v 1 . 0 )

Fire is the primary disturbance factor in many terrestrial ecosystems. Wildfire alters vegetation structure and composition, affects carbon storage and biogeochemical cycling, and results in the release of climatically relevant trace gases including CO 2, CO, CH4, NOx, and aerosols. One way of assessing the impacts of global wildfire on centennial to multi-millennial timescales is to use process-based fire models linked to dynamic global vegetation models (DGVMs). Here we present an update to the LPJ-DGVM and a new fire module based on SPITFIRE that includes several improvements to the way in which fire occurrence, behaviour, and the effects of fire on vegetation are simulated. The new LPJ-LMfire model includes explicit calculation of natural ignitions, the representation of multi-day burning and coalescence of fires, and the calculation of rates of spread in different vegetation types. We describe a new representation of anthropogenic biomass burning under preindustrial conditions that distinguishes the different relationships between humans and fire among hunter-gatherers, pastoralists, and farmers. We evaluate our model simulations against remotesensing-based estimates of burned area at regional and global scale. While wildfire in much of the modern world is largely influenced by anthropogenic suppression and ignitions, in those parts of the world where natural fire is still the dominant process (e.g. in remote areas of the boreal forest and subarctic), our results demonstrate a significant improvement in simulated burned area over the original SPITFIRE. The new fire model we present here is particularly suited for the investigation of climate–human–fire relationships on multimillennial timescales prior to the Industrial Revolution.


Introduction
Fire is one of the most important disturbance processes affecting the terrestrial biosphere. Fires affect most ecosystems from tropical forests to tundra (Bond and van Wilgen, 1996;Dwyer et al., 2000), and the evolution of fire over the Phanerozoic is likely to have had a major control on both the global carbon budget and the present distribution of plants on earth (Bond and Keeley, 2005;Pausas and Keeley, 2009;Bond and Scott, 2010;Bond and Midgley, 2012). Wildfires alter vegetation composition, structure, and distribution, biomass productivity, plant diversity and biogeochemical cycles (Johnson et al., 1998;Moreira, 2000;Ojima et al., 1994;Wan et al., 2001;Neary et al., 2005;Bond et al., 2005). Biomass burning (both wildfires and intentional combustion of biofuels) influences the spatial and interannual variability in the emissions of climatically relevant trace gases and aerosols, including CO 2 , CO,CH 4 , NO x , and black carbon (Crutzen and Andreae, 1990;Penner et al., 1992;Andreae and Merlet, 2001;Jain et al., 2006). Because of the close relationship between fire, vegetation, and climate, understanding the causes and consequences of fires is critical for any assessment of the past, present, and future state of the Earth system (Bowman et al., 2009).
Quantitative observations of wildfire are severely limited both in time and space; coverage in global data sets based on satellite observations began in the last decades of the 20th century. To improve our understanding of the drivers of fire and fire-related trace gas and aerosol emissions, it is therefore imperative to use numerical models of fire occurrence, behaviour, and impacts. But process-based modelling M. Pfeiffer et al.: Global biomass burning in preindustrial time: LPJ-LMfire (v1.0) of fire and fire emissions is challenging. Fire occurrence is influenced by climate, vegetation structure and composition, and human activities; fire behaviour is affected by weather, topography, and the characteristics of the fuel; fire disturbance alters vegetation composition and structure, and ultimately climate (Archibald et al., 2009;Bowman et al., 2009;Spessa et al., 2012). Thus, modelling fire requires representations of vegetation, fire, and climate that interact and feed back upon one another.
Mathematical models of wildfire dynamics have existed for over 40 yr (Rothermel, 1972). The original models of fire behaviour were motivated by needs for operational fire forecasting for firefighting and forest management applications. These models were applied at relatively small spatial scales of 10 0 − 10 3 ha, and have been extensively revised and updated over subsequent years (Burgan and Rothermel, 1984;Andrews, 1986;Burgan, 1987;Andrews and Chase, 1989;Reinhardt et al., 1997;Finney, 1998;Andrews et al., 2003;Andrews, 2007;Andrews et al., 2008;Heinsch and Andrews, 2010). Fire modelling at field scale is an essential part of fire management and mitigation worldwide, and modern operational fire models such as BehavePlus (Heinsch and Andrews, 2010) can be used for a wide range of fire management applications, including projecting the behaviour of ongoing fire, planning prescribed fire, assessing fuel hazard, and training.
More recently, fire models have been developed for application at larger spatial scales, e.g. for integration into dynamic global vegetation models (DGVMs) in order to simulate the fundamental ecosystem disturbance process that fire represents and, in some cases, to estimate the emissions of climate-relevant trace gases and aerosols at continental to global scale. Depending on the goals for application of the particular DGVM, the detail with which fire is represented varies, but all large-scale fire models include a representation of three key processes: 1. fire occurrence, 2. fire behaviour, and 3. fire impacts on vegetation.
The most complex representations of fire currently adapted for DGVMs incorporate and generalize many of the concepts and equations developed for operational fire forecasting models into a large-scale framework. The RegFIRM fire model (Venevsky et al., 2002), originally developed as an embedded module within the Lund-Potsdam-Jena DGVM (LPJ, Sitch et al., 2003), was one of the first global fire models that contained explicit representations of climatic fire danger and lightning-and human-caused wildfire ignitions. Building on RegFIRM, the SPITFIRE (SPread and InTensity of FIRE) fire model (Thonicke et al., 2010) included a more complete process representation of fire ignitions and behaviour, and further contained new representations of the im-pacts of fire on vegetation including plant mortality as a result of crown scorch and cambial damage, and routines for estimating trace gas and aerosol emissions. SPITFIRE was designed to overcome many of the limitations in previous fire models set within DGVM frameworks, and be flexible enough to permit simulation analyses at sub-continental to global scales with minimal input data requirements.
SPITFIRE is one of the most comprehensive fire modules for DGVMs currently available, and has been the focus of numerous studies on the role of fire in terrestrial ecosystems and the Earth system. Thonicke et al. (2010) presented the SPITFIRE model description and global assessments of simulated burned area and wildfire trace gas emissions. Gomez-Dans et al. (2013) used SPITFIRE in combination with MODIS burned area and tree cover data to improve the model's predictions of burned area at selected sites in different biomes using parameter calibration-optimization techniques. SPITFIRE has also been driven with L3JRC burned area data (Tansey et al., 2008) and MODIS burned area data (Roy et al., 2008;Roy and Boschetti, 2009) as part of the LPJ-GUESS vegetation model (Smith et al., 2001;Hickler et al., 2006) in a study examining emissions from biomass burning in Africa (Lehsten et al., 2009). Using LPJ-GUESS-SPITFIRE, Lehsten et al. (2013) examined how changes to fire frequency, including fire exclusion, affect tree-grass ratios in Africa. Recently, Spessa et al. (2012) benchmarked LPJ-GUESS-SPITFIRE against remote-sensing-based tree biomass data for pan-tropical forests and savannas (Saatchi et al., 2011;Baccini et al., 2012). The model was driven by a combination of monthly burned area from the Global Fire and Emissions Database (GFEDv3.1, Giglio et al., 2010; and long-term annual fire statistics (Mouillot and Field, 2005).
In addition to LPJ and its variants, SPITFIRE has been incorporated into other vegetation models. Spessa and Fisher (2010) coupled SPITFIRE to a global version of the Ecosystem Demography (ED) vegetation model (Moorcroft et al., 2001). ED has been incorporated into the MOSES2.2 land surface model (Met Office Surface Exchange Scheme, Essery et al., 2001;Fisher et al., 2010) and the Community Land Surface Model (CLM, Oleson et al., 2010). SPITFIRE is currently being integrated into ED-CLM (Spessa and Fisher, in preparation). With minor modifications, SPITFIRE has also been incorporated into the LPX-DGVM (Prentice et al., 2011) and applied in global experiments to quantify the contribution of wildfires to the global land-atmosphere CO 2 flux.
In the following sections, we describe LPJ-LMfire, which is a revised version of LPJ-SPITFIRE that we designed for simulating global fire and vegetation-fire interactions on centennial to multi-millennial timescales, primarily during prehistoric and preindustrial time. The purpose of this manuscript is to present a complete description of our current model code to facilitate referencing of the model in M. Pfeiffer et al.: Global biomass burning in preindustrial time: LPJ-LMfire (v1.0) 645 our future publications, and promote easier dissemination of our methods to other researchers who may be interested in using our model. We perform a detailed evaluation of the new model based on simulations and observations of fire in Alaska, and compare the results of a global simulation over recent decades to data sets of observed burned area. We conclude with recommendations for future model development.

Rationale for modifying SPITFIRE
We were motivated to modify SPITFIRE for two main reasons: (1) in some parts of the world with very little human impact on the landscape, most notably in boreal and subarctic North America, both LPJ-SPITFIRE and LPX simulated little or no burned area where observations show that large fires do occur, however infrequently. This indicated to us that the fundamental behaviour of the model and/or the data sets used to drive the model could be improved.
(2) We wanted to describe a scheme for simulating anthropogenic fire during the preindustrial period. The formulation for anthropogenic fire ignitions based on population density and a single spatially variable parameter a(N d ) did not seem appropriate to us based on what is known about the way humans used fire during preindustrial time. In updating SPITFIRE to tackle these goals, we had to make several changes to the fire module and to LPJ itself. In addition to these changes, we introduce new formulations for lightning occurrence, rate of spread in herbaceous fuels, and anthropogenic burning. A detailed description of our changes from the original SPITFIRE follows.

Methods
Here we present a new fire module, LPJ-LMfire, that is designed to be used with LPJ and similar DGVMs. The module is largely based on SPITFIRE (Thonicke et al., 2010), but has been substantially altered in a number of important ways. We made changes that improved the simulation of daily lightning ignitions, fuel bulk density, fire rate of spread, and fire mortality. In order to simulate human fire during preindustrial and prehistoric time, we replace the simple population-density-based formulation for anthropogenic ignitions with a classification of humans by their subsistence lifestyle and introduce specific goals for each group in terms of fire management of their landscape. We further introduce a new scheme to track the progression of individual fires over the entire fire season and simulate smoldering ignitions. Fires in LPJ-LMfire continue burning for multiple days once ignited and are extinguished only by changes in weather, by merging with other active fires, or by running out of fuel when encountering previously burned area. Finally, we account for passive fire suppression as a result of landscape fragmentation from anthropogenic land use. These new methods for calculating wildfire occurrence, behaviour, and impacts required changes not only to SPITFIRE, but also to LPJ, which we detail below.
In each section we detail the representations in LPJ-LMfire that are different from the original SPITFIRE, followed by any changes we needed to make to LPJ to accommodate the requirements of the fire model. The description below is intended to stand alone (i.e. the entire model can be reconstructed on the basis of the equations and parameters presented in this paper without relying on earlier published descriptions). A comprehensive list of abbreviations is provided in Table 1, a flowchart illustrating the structure of LPJ-LMfire depicted in Fig. 1, and a table listing the plant functional type (PFT)-specific parameters presented in Table A1. The remaining equations that were unchanged from original SPITFIRE are detailed in Appendix A, along with a table of supplementary symbols and abbreviations (Table A2).
As a note on random numbers, LPJ-LMfire, as with SPITFIRE and some versions of LPJ (e.g. Gerten et al., 2004), uses random numbers to calculate certain processes, including precipitation occurrence and daily precipitation amount. In LPJ-LMfire we additionally use random numbers in the calculation of lightning fire ignitions. In this paper when we describe the use of random numbers, we are referring to values drawn from a pseudo-random sequence that displays statistical randomness. To guarantee reproducibility of simulation runs in LPJ-LMfire across platforms, rather than using a built-in function, we include random number generators in the model code for sampling uniform distributions (Marsaglia, 1991), and for other distributions based on the uniformly distributed sequence (Dagpunar, 1988). We seed the random sequence at the beginning of each model run using a four-byte integer hash that is calculated from the geographic coordinates of the grid cell and is unique to at least 30 arc seconds of precision. The state of the random number sequence is stored separately for each grid cell, so the sequence of random numbers is preserved even if the model runs grid cells in parallel or a different order. This procedure ensures that every grid cell run with the same longitude and latitude will have exactly the same sequence of random numbers every time the model is run.

Factors excluding fire
As with SPITFIRE, the LMfire routine is designed to operate on a daily timestep. However, to save computation time  we implemented several checks to ensure that the fire routine is only called when fires are possible. We exclude fire when there is snow cover in the model, assuming that a snow layer will not allow the ignition and spread of surface fires.
As the current version of LPJ updates living biomass and the litter pools annually, we further skip calling the fire routine if the total vegetation foliar projected cover (FPC) of the grid cell is less than 50 %, or if the total amount of fuel, including live fuel, all four dead fuel classes and the soil surface carbon pool, is less than 1 kg m 2 . These thresholds, similar to those used in LPX (Prentice et al., 2011), are based on the assumption that if fuels are discontinuous or insufficient in quantity, a fire might start but will not be able to spread far enough from the starting point to cause a significantly large wildfire. We calibrated our thresholds by running the model for individual grid cells and evaluating the modelled fireline intensity (I surface ) in environments with low vegetation cover and/or total fuel load. These minimum fuel load and continuity thresholds are almost always met except in hot and polar deserts where vegetation reaches its bioclimatic limits.

Calculation of daily lightning ignitions
Lightning ignitions in SPITFIRE are calculated from a satellite-based climatology of monthly lightning flash density (Christian et al., 2003) that is interpolated between months and scaled to yield a quasi-daily climatology of lightning strikes (Thonicke et al., 2010). This daily number of lightning strikes is further reduced to fire ignitions based on a constant scaling factor. This approach takes into account neither the observation that lightning can be highly variable from year to year, particularly in regions where the total amount of lightning strikes is comparably low, nor that lightning occurrence is clustered in time (i.e. it is linked to precipitation events and times of atmospheric instability), nor that observations of fire ignitions suggest that a certain amount of stochasticity characterizes lightning-caused fires.
Here we describe our new approach for estimating the interannual variability of lightning, its daily occurrence, and a representation of the stochastic nature of lightning fire ignitions. Thonicke et al. (2010) argued that they expected the model sensitivity to inter-annual variability in lightning ignitions to be small compared to the overall model outcome and thus neglected interannual variability in lightning. However, we found that in places where fires are infrequent but important in terms of ecosystem impacts, and are generally caused by lightning (e.g. in boreal and subarctic North America), interannual variability in lightning occurrence is a key component of fire occurrence. In these regions, between 72 % and 93 % of all fires observed at present day are attributed to lightning ignitions (Stocks et al., 2003;Boles and Verbyla, 2000), and large interannual variability in burned area is visible in the GFEDv3 data set . Using the SPITFIRE or LPX formulations for lightning ignitions results in simulated burned area that is much smaller than observations in boreal and subarctic North America and Siberia, even though FDI is nonzero (Thonicke et al., 2010, Fig. 3c, Prentice et al., 2011. This inconsistency can be explained by the very low density of lightning strikes in the input climatology, which leads to an estimation of lightning ignitions that is well below one event per grid cell per month. We therefore believe that it is essential to capture interannual variability in lighting activity in order to simulate fire in boreal and subarctic regions that is consistent with observations. The only globally homogenized observation of lightning occurrence that is currently freely available is the LIS/OTD satellite-based data set (Christian et al., 2003), though other data sets, e.g. WWLLN (Virts et al., 2013) and GLD360 (Holle et al., 2011), are under development and could be applied in the future. The LIS/OTD data are available at the 0.5 • spatial resolution we use for LPJ-LMfire, but only as a climatology (the HRMC data set). Lower resolution LIS/OTD data are available as a multi-year monthly time series. However, for the extratropics (north and south of 42 • latitude) this time series, and the climatology, is based on only 4 yr of satellite observations. Because of the limited temporal coverage and low spatial resolution of available global lightning data, we developed a method of imposing interannual variability on climatological mean lightning frequency using ancillary meteorological data. Peterson et al. (2010) describe the correlation between convective available potential energy (CAPE) and cloud-toground lightning flashes for Alaska and northern Canada, indicating that lightning strikes are more common at times with positive CAPE anomalies. Based on this observation, we produce an interannually variable time series of lightning by scaling the climatological mean lightning flash rate with monthly anomalies of CAPE. The magnitude of the imposed variability is based on observed lightning strikes from the Alaska Lightning Detection System (ALDS, Alaska Bureau of Land Management, 2013).
To estimate the range of interannual variability in lightning amount, we analysed ALDS strike data for the time period between 1986 and 2010 for June, the peak lightning month in most of Alaska. Point observations of lightning strikes in the ALDS were aggregated on a 0.5 • grid, and grid cells with more than 5 yr of lightning strike observations (approx.  1750 valid cells) were analysed with respect to the minimum, maximum, and mean number of observed lightning strikes over all available years. For each grid cell, the minimum and maximum observed values were set into a ratio to the temporal mean. The two boxplots in Fig. 2 show the minimum-to-mean ratio and maximum-to-mean ratio distribution for all grid cells. The total range in interannual variability spanned four orders of magnitude, from 1 % of to 10times the mean. We used this range to scale climatological mean lightning strikes based on CAPE anomalies. Using CAPE from the 20th Century Reanalysis Project (Compo et al., 2011), we determined monthly anomalies on a grid cell level compared to the 1961-1990 mean CAPE value for a given month. The largest positive or negative CAPEanomaly value within the time series for a specific grid cell is used to normalize CAPE anomalies to a range between −1 and +1 for the entire time series available for a given grid cell. Applying the normalized CAPE anomaly with the scaling factor described above, the monthly number of lightning flashes is estimated as With the lightning flash density given by Eq.
(1), we disaggregate the monthly values to a daily amount and scale lightning flashes to cloud-to-ground lightning strikes. Noting that lightning and precipitation are closely correlated (e.g. Jayaratne and Kuleshov, 2006, and references therein;Michaelides et al., 2009;Katsanos et al., 2007), we allow lightning strikes to occur only on days with precipitation. Daily precipitation occurrence is simulated with a weather generator following the original SPITFIRE formulation (Thonicke et al., 2010). Simultaneous observations show that the quantity of lightning strikes is further positively correlated with precipitation amount (Piepgrass et al., 1982;Rivas Soriano et al., 2001;Zhou et al., 2002;Lal and Pawar, 2009). Therefore, to estimate the number of daily lightning strikes, we scale the total monthly lightning amount by the daily fraction of monthly total precipitation as simulated by the weather generator. With daily lightning flashes, we estimate ground strikes by using a flash-tostrike ratio of 20 % as in the original SPITFIRE. We confirmed this flash-to-strike ratio as realistic through a qualitative comparison of satellite-derived lightning flash density in the LIS/OTD LRMTS monthly time series with lightning ground-strike observations from the ALDS and from an extract of the North American Lightning Detection Network (NALDN, Orville et al., 2011) data set covering the southeastern United States.
With an estimate of lightning ground strikes, SPITFIRE calculates fire starts as a function of a fixed ignition efficiency of 4 %, yielding a total lightning flash-to-ignition ratio of 0.8 %. In contrast, the LPX fire model specifies a 3 % flashto-ignition ratio, and further reduces the number of fire starts 0.01 0.1 1 10 ratio of strikes to temporal mean Fig. 2. Maximum-to-mean ratio (top box plot) and minimum-tomean ratio (bottom box plot) for ALDS strike data in June between 1986 and 2010, based on approx. 1750 grid cells with more than 5 yr of observations. using the factor P +, which reduces the effectiveness of ignition events in wet months (Prentice et al., 2011, Eq. 1). Both of these methods result in a deterministic simulation of fire starts on any given day that is directly linked to lightning amount. The initiation of lighting-ignited fires is, however, also influenced by other factors, including the spatial distribution of lightning on the landscape, the temporal evolution of burned area during the fire season, and by a component that is observed but cannot be explained by large-scale variables, something that we term stochastic ignition efficiency.
These additional controls on fire starts are apparent when analysing patterns of lightning strikes and burned area in boreal and subarctic regions where lightning is rare but large fires develop; these are places where human impact is low, but both SPITFIRE and LPX fail to simulate burned area in agreement with observations. In attempting to improve our ability to model lightning-caused fire in the high latitudes, we made a series of changes to the way fire starts are calculated in LPJ-LMfire. Our new formulation accounts for the differential flammability of different plant types, fuel moisture, the spatial autocorrelation of lightning strikes, and previously burned area. All of these terms are combined to an estimate of ignition probability, against which we compare a uniformly distributed random number that represents the stochastic component of wildfire ignition.
Plant types differ in their intrinsic flammability as a result of leaf and stem morphology, typical canopy hydration status, and presence of phenols and other flammable compounds in the fuel (Diaz-Avalos et al., 2001). We noticed that treating all PFTs the same way with respect to ignition efficiency was problematic, especially when comparing the tropics (where lightning strikes are extremely frequent) to the extratropics (where fewer strikes appear in some cases to cause equal or more amounts of fire). In assigning PFT-specific ignition efficiency parameters, we took a top-down approach, where we qualitatively optimized the ignition efficiency parameter to match the performance of the model with respect to satellite-based observations of mean annual burned area fraction at the level of a few grid cells in areas where we judged human impact to be low (see Sect. 4.5, Fig. S9). This optimization of the parameters led to a large range of values between 0.05 and 0.5 (ieff pft , Table A1). The individual ignition efficiencies are combined into an FPC-weighted average ieff avg = npft pft=1 fpc grid ieff p ft npft pft=1 fpc grid . (2) Lightning strikes display a large degree of spatial autocorrelation, tending to cluster on mountaintops and other high terrain, tall buildings, water bodies, etc. Mazarakis et al., 2008;Uman, 2010). Because of this autocorrelation, successive thunderstorms over the course of a fire season become less likely to start new fires because lightning will strike places that have already burned. As such, we decrease the likelihood of lightningignited fires as a function of the area already burned to date: This equation is based on an empirical evaluation of NALDN data for Florida where we investigated the spatial autocorrelation of lightning strikes in relation to strike density. Similarly to LPX, the probability that a lightning strike will result in an ignition also depends on fuel moisture. LPX uses an additional parameter, β, based on a single transect across the Sahel and applied globally, to influence the relationship between fuel moisture and ignitions. Given the uncertainty in this formulation, and to avoid using another parameter, in LPJ-LMfire we use the fire danger index (FDI) as an indicator of fuel moisture. The overall ignition probability on a given day is therefore calculated as ieff = FDIieff avg ieff bf .
As explained above, this probability is compared with a uniformly distributed random number that represents the stochastic component of wildfire ignitions that helps to explain why in certain cases a single lightning strike can be sufficient to cause a fire, whereas in other cases many lightning strikes within one thunderstorm do not cause a single fire (Nickey, 1976;Keeley et al., 1989;Kourtz and Todd, 1991;Jones et al., 2009;Hu et al., 2010). The net effect of this approach is that lightning will sometimes cause a fire even though conditions are not very favourable, and vice versa. By allowing either zero or one ignition per grid cell and day, we account for the fact that lightning ignitions are discrete events.

Anthropogenic ignitions
Humans have used fire since the Palaeolithic as a tool for managing landscapes, optimizing hunting and gathering opportunities, cooking, hunting and defense, and communication (Pyne, 1994;Anderson, 1994;Pyne, 1997;Carcaillet et al., 2002;Tinner et al., 2005;Roos et al., 2010). The relationship beween humans and fire has changed over history, particularly after the Neolithic revolution when people began cultivating domesticated plants and animals (Iversen, 1941;Kalis and Meurers-Balke, 1998;Lüning, 2000;Rösch et al., 2002;Kalis et al., 2003), and during the 20th century following the widespread mechanization of agriculture and institution of industrial fire suppression. Since our goal is to develop a model capable of simulating fire in prehistoric and preindustrial time, we attempt to quantify the way in which humans in the past used fire. For us the main question is not simply how much fire people can cause, as it only takes a few dedicated individuals to cause significant amounts of fire (e.g. Eva et al., 1998), but rather -how much fire would humans want to cause, given certain environmental conditions and subsistence lifestyles? We further account for the physical limits to anthropogenic fire ignitions. Subsistence lifestyle is a very important factor determining why humans light fires and to what extent they light fires in order to manage their environment (Head, 1994;Bowman, 1998;. Hunter-gatherers use fire to promote habitat diversity and grass for game, keep landscapes open to ease their own mobility, and help prevent highintensity wildfires late in the season that could completely destroy vegetation resources. They accomplish these goals by lighting low-intensity fires early in the fire season that remove only understorey vegetation and prevent dangerous build-up of fuels (Lewis, 1985;Pyne, 1997;Williams, 2000;Kimmerer and Lake, 2001;Stewart et al., 2002). Pastoralists use fire to kill unpalatable species and stop woody encroachment, to promote the growth of fresh grass, to control parasites and animal movements, and to increase visibility while mustering (Crowley and Garnett, 2000;?). Farmers will burn crop residues after harvest and pastures for domesticated grazers and, depending on population density and availability of unused land, may use fire to prepare new cropland while old areas are abandoned, e.g. in systems of shifting cultivation.
Thus modelling human burning in preindustrial time is complex, as different groups of people had different goals for fire management and these probably changed in space and time, and because few quantitative observations exist that enable us to directly calibrate our model. It is therefore necessary to make assumptions on the relationship between humans and fire based on qualitative information, e.g. from ethnographic, anthropological, and archaeological studies. Theoretically, the only limit to how much people can burn depends on population density, average daily walking range of people, fire weather conditions, and fuel availability and structure. In most cases, people will not fully exploit the potential maximum amount of fires they can cause, as they will try to use fire in a constructive way to manage their habitat rather than destroying it by overburning (Head, 1994;Bowman, 1998;. We define this constructive use of fire in terms of burn targets for the three subsistence lifestyle groups described above. For foragers, we assume that their goal is to use fire to create and maintain semi-open landscapes, as this was the habitat most preferred by prehistoric people because habitat diversity and foraging opportunities increase with moderate disturbance, but decrease again if disturbance becomes too severe (e.g. Grime, 1973;Connell, 1978;Huston, 1979;Collins, 1992;Roxburgh et al., 2004;Perry et al., 2011;Faivre et al., 2011). We therefore link the annual amount that foragers will try to burn to the simulated degree of landscape openness, i.e. tree cover, and the effectiveness of fires to open up forest, i.e. the rate of change of vegetation cover over time. The annual burn target for foragers is calculated as with the change in grass cover being estimated as These equations imply that foragers living in an area with high forest cover will initially try to use fire to open the landscape. As the forest cover is reduced, the annual amount of anthropogenic fire will be reduced to maintain an equilibrium level of openness of the landscape. Alternatively, if anthropogenic burning has little effect on forest cover, e.g. in wet environments, humans will "give up" trying to burn their landscape after a short period of time. This quantification of hunter-gatherer fire use is based on suggestions that native North Americans repeatedly made controlled surface burns on a cycle of 1-3 yr, broken by occasional catastrophic fires that escaped the area intended to burn and periodic conflagrations during times of drought (Pyne, 1982;Williams, 2002b).
Pastoralists are assigned a constant burn target of 20 % (equal to a 5 yr fire return interval) that they will try to reach before they stop igniting fires, assuming that their interest in causing fires is less pronounced as they will try to preserve biomass for their domesticated grazers, while at the same time trying to maintain good pasture quality and avoid fuel accumulation in fire-prone environments. Present-day recommendations for prescribed fire maintenance of prairies and pastures suggest that a fire return interval target of 5 yr may even be on the more conservative side of estimates (Prairiesource.com, 1992; Government of Western Australia, Department for Agriculture and Food, 2005).
Farmers may burn unused land to expand their area under cultivation or prepare new fields as old ones are abandoned, e.g. in shifting cultivation systems. They may also light fires to control fuel build-up and mitigate the possibility of devastating wildfires in areas adjacent to their cultivated land, or use fire to maintain pastures. To account for these processes, we assign farmers an annual burn target of 5 % on land not used for agriculture, corresponding to a fire return interval of 20 yr.
Given the assumption that people burn purposely to achieve a certain goal, it is unlikely that all people who are present in a grid cell will cause fire. When 10 or more people are present in a grid cell, we therefore allow only every 10th person present to purposely ignite fires. Among all groups of people, cognitive, genetic, and economic factors mean that human social organization leads to hierarchies of group sizes. Numerous archaeological and ethnographic studies have demonstrated that these relationships are remarkably stable over time (e.g. Hamilton, 2007;Whiten and Erdal, 2012). Marlowe (2005) suggests that the optimal size of a hunter-gatherer group is 30 persons. We assume that three members of this group, e.g. able-bodied young males, will be responsible for fire management in the territory of the group. We allow for the possibility that the total number could be smaller at times, e.g. during colonization of new territory; if less than 10 people are present in a grid cell, then one person is responsible for fire ignitions. This 10 % scaling factor on active human agents of fire is most important when calculating ignitions among forager populations. In agricultural and pastoral groups, population density will nearly always be high enough to ensure that an overabundance of potential arsonists is available to aim for the burn targets we specify.
Anthropogenic ignitions are determined after the calculation of the average size of single fires and their geometry on a given day. The number of individual ignitions per firelighting person is calculated as where The area that one fire-lighting person potentially can burn in one day is given by the equation where the average distance that one person lighting fire walks in one day is limited to 10 km. How much fire people will start on a given day will depend on the environment in which they live. People who live in an environment that naturally has a lot of fire will take into account that some part of the landscape will burn naturally and adjust their burn target accordingly in order to avoid overburning. In order to take into account that people have a collective memory of the fire history in their habitat, we keep track of the 20 yr running mean of the burned area fraction in a given grid cell, and define the daily burn target for a given lifestyle group as target d,group = A gc max target y,group − bf 20 − burnedf , (10) with A gc being the grid cell area in ha. This function serves to reduce the target over the course of the year as people approach it. Once the target has been reduced to zero, people will stop igniting fires. The 20 yr-average burned area fraction is subtracted to let people stay conservative with their burning by taking into account that there can be some baseline amount of lightning-caused fire as well, thereby avoiding overburning of their target.
Ethnographic and historical studies have shown that preindustrial humans lit fires for landscape management purposes when fires were not likely to become severe, i.e. when meteorological conditions allowed burning but the overall fire danger was not too high. To represent this observation, we restrict anthropogenic burning to days when the average size of single fires,ā f , will not become larger than 100 ha. Additionally, the number of fires started by people on a given day is linked to the FDI via a multiplication factor that reduces the ignitions as FDI increases.
The decline of the risk factor, rf, follows a log-normal distribution with a maximum value of 1 at an FDI of 0.25 that then declines toward zero as FDI increases, which therefore makes it increasingly unlikely that people will keep causing fires when conditions for causing out-of-control fires become more risky. We developed this equation based on ethnographic studies from Australia showing that Aborigines preferentially cause fires at the beginning of the dry season when fire danger is still moderate, and decrease their ignition activities as FDI increases (Bowman, 1998;Yibarbuk et al., 2002;. We chose a log-normal curve to describe the relationship between anthropogenic ignitions and FDI because, even with high fire risk, the chance that someone causes a fire will not be completely zero. In cases where enough fire-lighting people are available to reach or exceed the burn target for the given day, the number of human-caused ignitions is derived from and in cases where the burn target of the day cannot be achieved due to a lack of enough fire-lighting people from Anthropogenic ignitions can be optionally specified for any given model run, but are always excluded in the model spinup before year 800 of the simulation in order to allow the development of a stable vegetation cover.

Burning of cropland
All of the equations presented in Sect. 3.1.3 concern anthropogenic burning on the fraction of the grid cell where potential natural vegetation is simulated by LPJ. We prescribe additional burn targets to account for anthropogenic burning on the part of the grid cell that is occupied by cropland. Evidence suggests that the usage of fire in cropland management was widespread in preindustrial times (e.g. Dumond, 1961;Sigaut, 1979;Otto and Anderson, 1982;Johnston, 2003;Williams, 2002a), and even nowadays is common in parts of the world where agriculture is largely unmechanized, e.g. in Sub-Saharan Africa, and parts of South and Southeast Asia, Indonesia and Latin America (Conklin, 1961;Seiler and Crutzen, 1980;Dove, 1985;Smittinand et al., 1978;Unruh et al., 1987;Kleinman et al., 1995;Van Reuler and Janssen, 1996;Cairns and Garrity, 1999;Akanvou et al., 2000;Fox, 2000;Rasul and Thapa, 2003).
Depending on agricultural practices, crop residues may be burned in situ or collected and burned throughout the year, e.g. as a fuel (Yevich and Logan, 2003). Fields that are burned may be burned immediately after harvest or shortly before planting, and in some places where double or triple cropping is practised possibly even several times per year. Cropland burning can be achieved largely independently of fire weather; for example, managed fire was historically important in places with hypermaritime climate such as the uplands of northwestern Europe (Mather, 2004;Dodgshon and Olsson, 2006).
In LPJ-LMfire, 20 % of the total simulated crop biomass produced within 1 yr remains on the fields as residues, and this remaining biomass becomes potential fuel for agricultural burning. Farmers are assumed to burn 20 % of the total cropland area within a grid cell every year. We derived this value from a qualitative comparison between total annual area burned observed in GFEDv3 and our simulated burning on natural land for regions in Africa where agricultural burning is commonly practised after harvest. It is a conservative first approximation for the past when people did not have modern-day technology available to prepare fields for the next crop planting after harvest, and likely could be much higher in places where, for example, multi-cropping is practised and all fields are burned after every harvest.
As described above, cropland and crop residue burning practices vary with space and time. We therefore make no attempt to estimate the seasonality of cropland burning, aside from excluding cropland burning when snow cover is present or temperatures are below 0 • C, and assume that burning is evenly distributed across all other days of the year. Future improvements to the model could attempt to resolve the temporal pattern of cropland burning by using a more sophisticated crop module for LPJ (e.g. Bondeau et al., 2007). For studies that focus on fire seasonality or trace gas emissions from biomass burning on a sub-annual scale, the timing of anthropogenic activities affecting seasonal patterns of fire cannot be neglected and will need to be accounted for explicitly.

Fire behaviour
As described above, boreal and subarctic regions are characterized by infrequent lightning ignitions that may still lead to large amounts of burned area because individual fires persist over the course of several weeks or months (Alaska Fire Service, 2013). On the other hand, both SPITFIRE and LPX (Prentice et al., 2011) allow fires to burn for a maximum duration of 241 min, after which individual fire starts are extinguished. Combined with the fractional occurrence of lightning ignitions described above, this representation of fire duration may be one of the main reasons why these models simulate burned area that is inconsistent with observations. The largest change we made from the original SPITFIRE was the implementation of a scheme for multi-day burning and the coalescence of fires. After making this fundamental change to the model, we had to revise other SPITFIRE formulations to make them consistent with our new approach. These revisions included changes to the representation of fuel composition and amount, to meteorological influences on fuel moisture and rate of spread, and the introduction of representation of the role of topography in influencing fire size. The new functionality and changes are detailed below.

Multi-day burning and coalescence of fires
Once a wildfire is started, it typically continues burning as long as fire weather conditions and availability of fuel do not restrict the progress of the fire (e.g. Todd and Jewkes, 2006;Desiles et al., 2007;Jones et al., 2009). Wildfires display a characteristic diurnal cycle, with the most active period being around midday and early afternoon when humidity is at a minimum and wind speeds are higher (Pyne et al., 1996). To account for these observations, we remove the 241 min limitation on fire duration specified in SPITFIRE, but maintain this value as an active burning period on any given day in calculating daily burned area. Individual ignitions persist from one day to the next until they are extinguished due to (1) merging with other fires, (2) running out of fuel from burning into areas already burned during the current year, or (3) as a result of sustained precipitation.
In LPJ-LMfire, the total number of fires burning on a specific day is therefore defined as the number of fires that were started on previous days that have not yet been extinguished, plus any potential additional ignitions on the current day. As individual fires grow in size, the likelihood of one fire burning into another or into an area that has already burned increases. To take this into account, we reduce the number of fires burning on any given day by the product of the grid cell fraction that has already burned in the current year and the total number of fires on this day. Thus, the total number of fires on any given day is calculated as In allowing fires to burn for multiple days, we needed to define threshold amounts of precipitation above which ongoing fires will be extinguished. Field observations have shown that while small amounts of precipitation will impede fire spread, fires may keep smoldering and start spreading as soon as conditions dry out again, and that the amount of precipitation required to slow or stop wildfires differs depending on the type of fuel that is burning (Latham and Rothermel, 1993;Hall, 2007;Hadlow, 2009;Pyne et al., 1996). LPJ-LMfire extinguishes burning fires when the precipitation sum over consecutive days exceeds 10 mm for grid cells that have a grass cover of less than 60 %, and 3 mm for grid cells with more than 60 % grass cover (i.e. fires are extinguished after as many rain days in a row as it takes to reach the extinction threshold).

Fuel quantity and density
While testing development versions of LPJ-LMfire, we noticed that simulated burned area greatly exceeded GFEDv3 observations in parts of Siberia and the seasonal tropical forests of South America. We diagnosed the cause as very high simulated fuel loads that in turn propagated extremely large fires. High fuel loads in the tropics were the result of unrealistic accumulation of biomass in living vegetation, whereas in the boreal regions slow decomposition of litter with low bulk density led to an unrealistically deep and loosely packed fuel bed. To improve the simulation of fire, we therefore made several changes to the way LPJ simulates biomass and fuel bed density.
In LPJ, the amount of live woody biomass in a grid cell is determined by the PFT state variables of the average individual that represents the mean of the PFT population with respect to all state variables describing the PFT, and by the individual density that represents the number of individuals in a unit area (Sitch et al., 2003). Accumulation of biomass in the average individual is limited by the maximum crown area parameter. Density is limited by space in the grid cell, with the assumption that individuals do not overlap in space (packing constraint). Thus, at equilibrium, individual density stabilizes as the size of the average individual approaches maximum crown area. In our tests, simulated biomass accumulated to very high levels in areas where disturbance is rare and growth rates are high, such as the perennially humid parts of the Amazon Basin.
To reduce biomass in LPJ-LMfire, we allow trees to reach a maximum crown area of 30 m 2 , instead of the 15 m 2 used in the original LPJ parameterization. At the same time, we increased the maximum sapling establishment rate from 0.12 individuals m −2 to 0.15 individuals m −2 . As leaves have less biomass per unit area than stems, increasing the maximum crown area parameter in the model decreases density and therefore simulated total biomass. Adjusting these two parameters leads to an overall decrease in total biomass between 5 % and 15 % for the area shown in Fig. 3, with highest reduction percentages in areas of high biomass such as the upper Amazon Basin. As described above, the reduction effect caused by the increase of maximum crown area is most relevant for the wet tropics where trees experience little disturbance and optimal growth conditions. In most extratropical regions, the new limit for maximum crown area is usually not reached due to climate-induced mortality and disturbance.
In boreal regions where we noticed very high amounts of burned area in our development simulations, we traced this back to high rates of fire spread simulated in an unrealistically deep and loosely packed fuel bed. In LPJ, litter decomposition is controlled by temperature and moisture so that under cold, dry conditions, very slow effective decomposition rates are simulated and litter tends to accumulate for decades to centuries. In boreal regions, particularly in the drier parts of Alaska and Siberia, the model therefore simulated large accumulations of aboveground litter, with values as high as 7 kg C m −2 . Following the original SPITFIRE parameterization, fuel bulk density is relatively low: 2 kg m −3 for herbaceous litter and 25 kg m −3 for woody litter. Large accumulations of litter therefore lead to the formation of a deep, loosely packed fuel bed. This problem is exacerbated when frequent fires result in widespread tree mortality and shift the vegetation cover towards being dominated by herbaceous PFTs.
Cold, dry climates lead to the accumulation of large amounts of organic matter, but the assumption that these would not be mechanically and chemically altered with time is unrealistic (Berg, 2000;Berg et al., 2001;Akselsson et al., 2005).To account for changes in the physical properties of the fuel bed with time, we introduce an aboveground organic matter pool in LPJ that schematically represents an O horizon. After having calculated decomposition in the three litter pools (fast litter, slow litter and belowground fine litter) following Sitch et al. (2003), the remaining carbon in the fast litter pool is transferred to the O horizon where it decomposes with a nominal turnover time of 2 yr at a temperature of 10 • C. This way, an organic layer can build up in cold places where litter decomposition is slow, and unrealistically large accumulations of litter are avoided. Carbon that was transferred to the O horizon does not contribute to the rate of spread calculations as it is considered to be densely packed compared to the fuels in the regular fuel size classes, but it is included into the overall fuel combustion term. As shown in Table 2, reducing the amount of dead fuel by transferring older litter into the O horizon strongly affects the simulated rate of spread, and therefore fire size and burned area. We also noticed that our implementation of the original SPITFIRE resulted in high rates of fire spread in tundra ecosystems and consequently simulation of burned area that exceeded observations (GFEDv3, Alaska Fire Service, 2013). As the standard version of LPJ does not have a tundra shrub PFT, subarctic vegetation is primarily represented by the C 3 -grass PFT, for which SPITFIRE assigns a constant fuel bulk density of 2 kg m −3 . In tundra ecosystems, herbaceous plants and shrubs grow close to the ground and typically have a dense life form, e.g. as tussocks, as an adaptation against damage from frost and snow burden (Bliss, 1962;Sonesson and Callaghan, 1991;Sturm et al., 2000). To account for the dense growth form of tundra and the general tendency of herbaceous vegetation to grow more densely and closer to the ground with decreasing temperatures, we introduced a dependency between the bulk density of the two herbaceous PFTs and the 20 yr running mean of the annual sum of degree-days on a 5 • C base (GDD20, Sitch et al., 2003): In the tropics, the annual GDD sum can be as high as 10 000, whereas in high latitudes values are typically 1000 or less.
With fewer GDDs, we decrease bulk density from typical values in tundra areas of 10-12 kg m −3 to 1-2 kg m −3 in warm tropical regions where tall grasses grow. These endpoint values are estimated based on abundant field evidence demonstrating that tropical grasses are typically tall, whereas herbaceous tundra is short and often grows in dense tussocks (e.g. Breckle, 2002;Gibson, 2009). We use GDD20 because grass bulk density should not be influenced by interannual variability in climate, as individual species have a relatively stable growth habit over time. The modification of grass fuel bulk density affects simulated rate of spread. For example, given a fuel load of 1 kg m −2 , a wind speed of 3 m s −1 , and a fuel bulk density of 2 kg m −3 , the resulting ROS is 2.36 m s −1 at an rm of 0.1 and 1.22 m s −1 at an rm of 0.5. With a fuel bulk density of 12 kg m −3 , ROS is reduced by roughly one order of magnitude to 0.27 m s −1 and 0.14 m s −1 .

Fuel moisture
For herbaceous fuels, we set the relative moisture content of the fuel to be equal to the ratio where ω nl is the mean relative moisture content of the 1 h fuel class and the live grass, and me nl is the mass-weighted average moisture of extinction for live grass and 1 h fuel. ω nl and me nl are calculated as follows: As discussed above, the implementation of multi-day burning in LPJ-LMfire led to simulations of fires that were overly large and frequent compared to observations. This overburning was partly solved by introducing the O horizon for surface litter, and by adjusting the bulk density of live herbaceous fuels. However, in drier boreal and subarctic regions, we also noticed that herbaceous live fuel moisture was very low in the middle of the growing season. This low moisture was a result of LPJ's standard representation of soil hydrology where all soils are considered to be free draining. In reality, much of the boreal and subarctic regions are underlain by permafrost, which acts as a barrier to water drainage (Kane and Stein, 1983;Niu and Yang, 2006). To approximate the effects of permafrost on soil moisture, and therefore herbaceous live fuel moisture, we impede all drainage of soil water in LPJ where permafrost is present. We define permafrost as occurring in any grid cell where the 20 yr running mean annual temperature is less than 0 • C. For woody fuels, relative moisture content is calculated as Instead of resetting the relative daily litter moisture to saturation as soon as daily precipitation exceeds 3 mm, i.e. when the Nesterov Index (NI) is set to zero, we calculate ω o as a mass balance between drying and wetting of the fuel, assuming that at a threshold of 50 mm precipitation all fuel will be completely wet, and lesser amounts of rain will partially wet the fuel according to the amount of precipitation. The drying term is estimated as a function of daily maximum and minimum temperature similar to the way the Nesterov Index is calculated in original SPITFIRE, based on the difference between the day's minimum and maximum temperature, the fuel water content, and a fuel drying parameter integrated over the α-parameters given in Thonicke et al. (2010) according to fuel composition: wet = 1, prec > 50 mm prec 50 , prec ≤ 50 mm , with 50 mm of daily precipitation being the threshold definition for heavy rain given by the World Meteorological Organization (http://severe.worldweather.org/rain/), at which we assume all fuel to be water-saturated, independent of its previous water status. The water balance between drying and wetting is calculated as follows: which is essentially a simple water bucket approach similar to the way the soil water balance is calculated in LPJ. The fuel moisture on the current day is defined as The variable caf represents α combined over all fuels, and is calculated as The mass-weighted average moisture of extinction over all fuels, me avg , is calculated as Depending on the grass cover fraction of the grid cell, FDI is calculated as

Fire rate of spread
In contrast to SPITFIRE, we assume that fires will be mostly carried in light fuels as these are easily ignited due to their high surface area-to-volume (SAV) ratio and low fuel bulk density, whereas heavier fuel components will sustain burning once fire has started at a given place. As each PFT in LPJ occupies an exclusive space on the grid cell, the possibility that their fuels are spatially collocated is also excluded. Our Monte Carlo simulations on the continuity of natural land depending on the fraction that is occupied by agricultural land (Sect. 3.2.6, Eq. 33) revealed that, in a randomly distributed spatial arrangement of two differing entities, the fractional occupation ratio has an influence on the continuity of both entities. This result also applies to the distribution of herbaceous versus woody PFTs and, thus, fuels. For example, if a herbaceous PFT occupies more than 60 % of the grid cell, fire rate of spread is determined by the properties of the herbaceous fuel because it is not possible to arrange the remaining 40 %, i.e. the woody PFTs, in a way that interrupts the continuity of the herbaceous fuel. Below 60 % herbaceous cover, the average contiguous size of patches of herbaceous vegetation rapidly decreases as long as areas occupied by grass or trees are assumed to be distributed more or less randomly, and the influence of woody fuels on the overall rate of spread becomes more dominant. We therefore calculate rate of fire spread for herbaceous and woody fuel components separately and then average the two calculated rates of spread according to the coverage of the herbaceous and woody PFTs on the landscape.
To calculate rate of spread in grass, we use a modified form of the equation given in Mell et al. (2012), setting the fuel bulk density for these light fuels equal to the ρ livegrass value calculated in Eq. (15).
Equation (28) accounts for the variable density of live grass depending on GDD20 as calculated in Eq. (15). Compared to SPITFIRE, the rate of spread in this new equation requires fewer parameters (wind speed, ratio of relative fuel moisture to its moisture of extinction, and fuel bulk density) and typically results in slower rate of spread when all other conditions are equal. The rate of spread in woody fuel is calculated as in SPITFIRE, with the exception that we use a fixed value of 5 cm 2 cm −3 for SAV assuming that fire will be carried primarily by the finest component of the fuel bed. For details on the calculation of rate of spread, see the equations in Appendix A.
We determine the surface forward rate of spread as the weighted average of the rate of spread in the woody and herbaceous fuel according to the cover fractions of tree-and grass-PFTs on the landscape: ROSf s = ROSf sw treecover + ROSf sg grasscover treecover + grasscover .
In addition, we introduced a wind multiplier for high-wind conditions: at a wind speed of 10 m s −1 and above, the calculated ROS will be doubled, as the BEHAVE-based ROS is increasingly too low at higher wind speeds (see Fig. 13 in Morvan et al., 2008):

Effect of terrain on average fire size
Terrain can be an important factor influencing the spread of fires (Pyne et al., 1996). We argue that areas with high relief energy should have smaller average fire sizes compared to areas that are completely flat as dissected topography will inhibit fire propagation. Although fire rate of spread is usually faster upslope due to more fuel surface being exposed to the flames than on flat terrain and additional upslope wind effects, at 0.5 • spatial resolution no individual grid cell of ∼1000-3000 km 2 represents one single slope. Rather, all upslopes will be accompanied by downslopes on the opposing side where fire spread will be slowed or impeded. Terrain with high relief energy is also characterized by varying slope exposures. A dry sun-exposed slope will be opposed by a shady slope with wetter fuel conditions, different vegetation, and in some cases a sparsely vegetated crest that separates both slopes and impedes the spread of fires from one catchment into a neighbouring one (Guyette et al., 2002). Fuel continuity also can be broken by areas of unvegetated rock and cliffs, which are more likely to occur in complex terrain. Our qualitative observations of remotely sensed burned scars (Alaska Fire Service, 2013), databases of individual fire size (National Interagency Fire Service, 2013), and previous modelling studies (Parks et al., 2012) show that very large fires, i.e. those that would consume an entire 0.5 • grid cell, are rare in mountainous regions. To capture this effect, we calculate a terrain impedance factor slf = which affects mean fire sizeā f as a downscaling factor: We determined the median slope angle γ of a 0.5 • grid cell by aggregating the maximum D8 slope (Zhang et al., 1999) at 1 arc minute resolution from the ETOPO1 global digital elevation model (Amante and Eakins, 2009). Median slope angle at this scale ranges roughly from 0 • to 17 • from horizontal. A world map of slf is shown in Fig. S2.
With the size of individual fires scaled according to the average slope angle, more fires will be required to burn an equivalently sized total area in more complex terrain as compared to flat terrain.

Passive fire suppression through landscape fragmentation
For the first time in human history, modern technology allows people to actively suppress and extinguish wildfires to protect their lives and properties. In the past, possibilities to actively suppress and extinguish wildfires were limited (Skinner and Chang, 1996;Pausas and Keeley, 2009). Nevertheless, increases in population densities and parallel increases in land use eventually contributed to landscape fragmentation and thereby indirect suppression of wildfires. Following Archibald et al. (2009) we simulate the effect that anthropogenic landscape fragmentation has on fire spread and therefore burned area. In order to estimate the effects of anthropogenic landscape fragmentation, here defined as the fraction of cropland vs. unused land, we performed a Monte Carlo simulation on a grid of 100 × 100 pixels where we increased the fraction of cropland by 1 % increments from 0 to 1. For each step, we randomly assigned pixels within the grid to either be cropland or unused land and calculated the average contiguous area size of natural patches based on an 8-cell neighbourhood. To estimate the final average contiguous area size of natural patches, we performed 1000 repetitions of the experiment at each land use fraction. The resulting relationship between the cropland fraction of a grid cell and the average contiguous area size of unused patches can be approximated by the following equation: ac area = 1.003 + e (16.607−41.503f nat ) −2.169 A gc , with A gc being the grid cell area in ha. The equation accounts for changing land use as fragmentation is recalculated every year based on the information on how much land within a grid cell is agricultural land. The average contiguous area size of natural patches is used to set an upper limit toā f , the size of individual fires, in the fire routine. At very high land use fractions, we limit the minimum allowed averaged patch size to a kernel size of 10 ha, not allowing any fragmentation that causes natural patches smaller than this size. The concept of connectivity and fragmentation being related to the proportions of two different phases, in our case agricultural land and unused land, is well known in other scientific contexts, e.g. in soil science where unsaturated soil water conductivity depends on the ratio between water-filled and air-filled pore space (Richards, 1931;Newman and Ziff, 2000). For a detailed depiction of the Monte Carlo simulation results, see Supplement, Fig. S1.

Fire mortality
Fire mortality in the original version of SPITFIRE was simulated through a combination of cambial damage and scorching of tree crowns following Peterson and Ryan (1986), where tree kill is a function of fire intensity, bark thickness and tree height. Thus, to simulate realistic amounts of tree kill, it is essential to have a representation of the size and shape of trees in the model that is realistic. However, the population averaging of the allometric equations in LPJ leads to the simulation of average individuals that are much shorter and thinner than mature trees in nature. To overcome this limitation, SPITFIRE applied an unpublished scheme to disaggregate the biomass represented by the average individual into a series of size classes with height and diameter that are relative to the height of the average individual simulated by LPJ. We use an adaptation of this scheme to approximate realistic tree heights in LPJ-LMfire. We begin by prescribing a PFT-specific relationship between the simulated range in height for the average individual and the typical range in height from sapling to mature tree of a real individual of that PFT as it is observed in the field. Thus any given height of the average individual can be mapped to a mean real height (H real ) for the PFT. Recognizing that the average individual represents a range of tree ages and sizes, we disaggregate the biomass of each average individual into seven height classes following a skewnormal distribution centred onH real estimated above. The heights of each height class are equally spaced and range from 50 % ofH real for the shortest class to 125 % ofH real for the tallest class. Stem diameter is calculated separately for each height class based on the observed relationship between maximum tree height and diameter for each PFT. Bark thickness is calculated using the PFT-specific bark thickness parameters given in Thonicke et al. (2010) (par1, par2, Table A1). As in SPITFIRE, mortality resulting from cambial kill is calculated separately for each height class, and the total mortality over all classes is summed up across all classes per PFT. Apart from bark thickness, the probability of mortality due to cambial damage also depends on the residence time of the fire, τ l , in relation to the critical time for cambial damage. Thonicke et al. (2010) do not provide the exact equation used in SPITFIRE to calculate τ l , but refer to Peterson and Ryan (1986). In LPJ-LPMfire we calculate τ l using Eq. (8) of Peterson and Ryan (1986): With our revised height class scheme, we needed to reparameterize the PFT-specific RCK-and p values that describe the probability of mortality due to crown damage. When we used the SPITFIRE RCK parameters, close to 1 for all woody PFTs with the exception of the tropical broadleaf raingreen PFT, an undesired result of our multipleday burning scheme was that excessive crown kill resulted in much of the simulated global vegetation cover being converted to grasslands in places with frequent fire occurrence. Observational data, e.g. from vegetation maps and the Global Land Cover Facility (GLCF) tree cover data set (DeFries et al., 2000), showed that many of these places clearly should be forested. While we acknowledge that using parameters from observed plant traits is a good strategy, given the unrealistic allometry simulated for LPJ's average individual and the simplification presented by our height class scheme, direct representation of the characteristics of individual trees is not strictly possible. Future model development should include better representation of the size and shape of trees in the model, e.g. by using a cohort-based approach such as that used in LPJ- GUESS (Smith et al., 2001). In LPJ-LMfire, we set RCK to a constant value of 0.5 for all tree PFTs and p to a constant value of 0.3. We further add the restriction that deciduous trees can only be killed by crown scorch if green leaves are present at the time of fire occurrence. In nature, most grasses grow quickly enough to finish their life cycle within one growing season (Gibson, 2009). Some herbs and grasses are annual species that sprout from seeds every year, while for many perennial herbaceous plants the entire aboveground biomass dies back after the growing season and then resprouts from the root mass during the next growing season (Cheney and Sullivan, 2008;Gibson, 2009). In LPJ however, herbaceous PFTs take 3-10 yr to reach equilibrium potential aboveground biomass under constant climate, soil and CO 2 forcing in part because establishment and allocation are updated only once annually. In SPITFIRE, herbaceous biomass is removed as a result of combustion. In areas with frequent fire, LPJ-SPITFIRE simulates herbaceous biomass and FPC that are lower than observations. This inconsistency affects not only fire behaviour but also general biogeochemical cycling in ecosystems where herbaceous vegetation is present.
To avoid an unrealistic reduction in herbaceous biomass in LPJ-LMfire as a result of fire, we convert combusted live grass biomass to carbon, but do not remove the grass biomass from the live biomass pool at the end of year, similarly to the scheme used by Kaplan et al. (2011) to simulate the harvest of agricultural crops. This correction results in more realistic biomass and coverage of grasses when simulating fire. In the future, a new and more realistic implementation for the development and senescence of grasses within LPJ should be implemented, which will require moving to a daily time step for grass allocation, as, for example, has been done for crops in LPJ-ML (Bondeau et al., 2007).

Data sets and model runs used for model evaluation
Evaluating a complex DGVM and fire model such as LPJ-LMfire requires suitable input data for driving the model, including information on climate including lightning, soils, topography, atmospheric CO 2 concentrations, and human population density and anthropogenic land use. Unfortunately, not all parts of the world where fire is observed are equally well represented in terms of quality data for driving and testing DGVMs with fire. In the simulations described below, we prepared a standard, global driver data set for LPJ-LMfire using the data sets listed in Table 3. To drive the model with the best possible approximation of actual climate conditions, we use a baseline long-term mean climatology with a native spatial resolution of at least 0.5 • to which interannual variability is added in the form of anomalies from a lower resolution reanalysis climate simulation that covers the period 1871-2010. We calculated anomalies in the reanalysis data relative to a 1961-1990 standard period, and linearly interpolated the 2 • reanalysis grid to 0.5 • using the CDO software (Schulzweida et al., 2012).
In all of the simulations presented in this paper, the model was spun up for 1020 yr with a detrended version of the 20th Century Reanalysis climatology, with the atmospheric CO 2 concentrations of 1871, and then run in a transient simulation from 1871 to 2010. For the Alaska case study we replaced LIS/OTD with the ALDS data set for the time period of record that overlapped with our experiments .
Since we focus on the overall performance of the model in simulating fire behaviour and impacts on ecosystems, and since the development of the demographic history data sets is the subject of a separate publication, we exclude anthropogenic ignitions from the simulations presented here.
We needed model-independent data to evaluate simulated fire frequency and behaviour, e.g. satellite-derived or groundbased data of annual burned area. To evaluate LPJ-LMfire's performance in Alaska, we compared simulated area burned between 1986 and 2010 with the AFS historical burned area polygon data set (Alaska Fire Service, 2013). For global model evaluation, we used GFEDv3  and the global burned area data set published by Randerson et al. (2012).

Model results and evaluation
In the following sections, we first present and discuss LPJ results for simulated aboveground biomass and the O horizon.
We then present our case study for Alaska where we evaluate LPJ-LMfire simulation results with reference to the highquality data sets on lightning strikes that we used to drive the model, and detailed maps of annual burned area that we used to test model output. We present and discuss a world map of potential natural fire return interval that could be used for ecosystem management and restoration and, finally, compare a global fire scenario to global observations of burned area.

Aboveground biomass
As noted in Sect. 3.2.2, living aboveground biomass simulated by LPJ was consistently overestimated compared to values reported in literature, especially in places with high biomass such as the Amazon Basin, where simulated values reached a maximum of more than 30 kg C m −2 . After the modifications we made to maximum crown radius and maximum establishment rate, aboveground biomass simulated in the central Amazon Basin ranged between 18 and 21 kg C m −2 (Fig. 3a). Comparisons of our simulated biomass with satellite-derived observations (Saatchi et al., 2009) show that even after the modifications, LPJ's estimates of aboveground live biomass are likely to be still on the high end of estimates. Aboveground biomass carbon estimates collected by Malhi et al. (2006) for old-growth Amazonian forests range between 8.5 and 16.7 kg C m −2 . Estimates of biomass carbon for tropical moist forests in the Brazilian Amazon collected by Houghton et al. (2001) range between 10 and 23.2 kg C m −2 , with a mean of 17.7 kg C m −2 . In regions with generally lower biomass, e.g. in the Caatinga of northeast Brazil or in the Andes, simulated and satellitederived biomass values reported by Saatchi et al. (2009) are generally in good agreement, although the model underestimates biomass in parts of the Andes. Figure 4 shows the global amount of carbon stored in the new LPJ O horizon. The highest values are found in northeastern Siberia and northern North America, with values ranging between 2 and 3.5 kg C m −2 . In northern Europe, simulated values range between 1 and 2 kg C m −2 . These values do not capture the high end of values reported in literature, but are well within the observed range. For example, Mäkipää (1995) reported a range of 0.5 to 3 kg C m −2 for the organic layers of forest soils in southern Finland, depending on nutrient status and site wetness. For the arctic tundra of North America, Ping et al. (2008) reported values as low as 0.7 kg C m −2 for mountain sites, and reaching 15.1 kg C m −2 for lowland sites. Pregitzer and Euskirchen (2004) summarize organic soil horizon stocks from a number of studies, giving a range between 0.2 and 19.5 kg C m −2 for boreal forests. The values simulated by LPJ are therefore within a realistic range, although site-specific variability cannot be reproduced at 0.5 • spatial resolution.

Fire in boreal ecosystems: the Alaska case study
Fire is an important process in the boreal region and controls a variety of different ecosystem processes such as succession, tree recruitment, vegetation recovery, carbon storage, soil respiration and emission of atmospheric trace gases (Landhaeuser and Wein, 1993;Kurz and Apps, 1999;Johnson, 1992;Harden et al., 2000;Turetsky et al., 2002;Bergner et al., 2004;Kasischke et al., 2005). Alaska was particularly suitable for our model evaluation, first because neither SPITFIRE nor LPX was able to simulate adequate amounts and realistic variability of burned area in boreal and subarctic environments, and also because the availability of data to drive and evaluate the fire model is excellent for this region.
Because sufficiently dry conditions occur comparatively rarely, fire is highly episodic in boreal and subarctic Alaska and northern Canada (Kasischke et al., 2002), and hence the observational record is dominated by relatively few big fire years. Lightning is the main source of ignitions for large fires in boreal ecosystems. For the period 1950-1969, Barney (1971) showed that ∼ 24 % of all fire ignitions in Alaska were caused by lightning, but fires started by lightning accounted for more than 80 % of total area burned. Todd and Jewkes (2006) provide an extensive year-by-year overview from 1950 to 2005, listing the total number of wildfires per year caused by humans and lightning, and the corresponding number of acres burned by these wildfires. A total of 89 % of all burned area between 1950 and 2005 can be attributed to lightning-caused fires (Todd and Jewkes, 2006). From 1986 to 2005, 11 yr had more than 95 % of the total annual area burned attributed to lightning fires, 13 yr more than 90%, and 16 yr more than 80 %. One of the reasons why the highly variable fluctuations in burned area could not be reproduced by the original version of SPITFIRE could be because interannual variability in lightning occurrence was neglected as described in Sect. 3.1.2 above. Furthermore, smoldering fires are an important part of fire behaviour in boreal and subarctic environments. For example, the recent Anaktuvuk River tundra fire smoldered for nearly two months as the tundra dried out before spreading rapidly at the end of the summer (Jones et al., 2009). With the high-quality data sets that are available on fire in Alaska, we set out to see if the improvements we made to LPJ-LMfire substantially improved the model performance in this ecologically important region.

Simulated and observed area burned
Since the majority of burned area in Alaska is due to lightning-ignited fires (Todd and Jewkes, 2006), we set the model up only to simulate ignition and spread of natural, i.e. lightning-ignited, fires on land not subject to human land use. We distinguish the following seven major ecoregions (Fig. 5)  Depending on the ecoregion in consideration, the simulated and observed area burned on average over the time period from 1986 to 2010 varies considerably. In the following sections, we compare and discuss simulated fire occurrence with observed burned area by ecoregion.

Intermontane Boreal ecoregion
The Intermontane Boreal ecoregion, situated between the Alaska Range and the Brooks Range, is the most important region of Alaska for fire. On average, 93 % of the total area burned in Alaska is located in this area. Both the observational data and the simulation results identify this area as the region most affected by fire. In this region, observations show an average annual burned area of 4834 km 2 over 25 yr and a standard deviation of 6285 km 2 or 0.96 ± 1.25 % of the total area of the region (Table 4). Our simulated annual burned area of 4736 ± 5654 km 2 , or 0.94 ± 1.13 %, agrees well with observations, slightly underestimating both the total amount and the magnitude of the interannual variability in burned area. The absolute range of area burned in this region is approximately the same for both the observations and simulation, with a minimum of 136 vs. 0 km 2 and a maximum of 26 464 vs. 25 500 km 2 , respectively (Fig. 6). For both observations and simulation, the annual mean burned area is larger than the median, indicating that the annual fire regime is characterized by relatively low area burned, occasionally interrupted by extreme years during which large areas burn. In contrast to the mean, where simulated burned area is slightly less than observations, the median and 75 % percentile burned area are slightly higher in the simulation than in the observations (Fig. 6).  In Fig. 7 we show the simulated and observed time series of burned area in the Intermontane Boreal region. LPJ-LMfire reproduces observations of burned area well not only in terms of the average area burned over the 25 yr period, but also in terms of the interannual variability.

Arctic Tundra
Compared to the Intermontane Boreal ecoregion described above, burned area in the other six ecoregions is very small in terms of total area burned as well as percent of the ecoregion burned (Fig. 6, Table 4). Our simulations therefore correctly identify the location of the most important ecoregion for fire in Alaska. However, our simulations overestimate the mean annual area burned as well as the maximum annual area burned for ecoregion AT (Arctic Tundra) compared to the observation data. This is due to 2 yr within the simulated time series, 2008 and 2009, for which we largely overestimate the total area burned, whereas in most other years we simulate low amounts of burning that match the observational data in magnitude and variability. Exceptional years with very large single tundra fires are known to occur, e.g. the Anaktuvuk River fire in 2007 (Jones et al., 2009). Although LPJ-LMfire is capable of simulating years with exceptionally large amounts of fire in Alaska's arctic tundra, we are not able to reproduce burned area in exactly those years when large burned area was observed.

Bering Taiga and Bering Tundra
Burning in the westernmost part of Alaska (ecoregions BTA and BTU) is generally low in the observational data (Fig. 6, Table 4), with a maximum of 675 km 2 burned during the period 1986-2010, with an average of 86 km 2 yr −1 , and a median of 27 km 2 yr −1 for the Bering Taiga, and a maximum of 367 km 2 yr −1 , an average of 48 km 2 yr −1 and a median of 0 km 2 yr −1 for the Bering Tundra. This implies that an average of 0.03 % of the Bering Taiga and 0.05 % of the Bering Tundra region burned over the 25 yr period. Our simulations underestimate burning in these regions, especially for the Bering Taiga, where the simulated maximum burned area is 329 km 2 yr −1 , with an average of 22 km 2 yr −1 and a median of 0 km 2 yr −1 . For the Bering Tundra, we simulate a maximum of 148 km 2 yr −1 , an average of 15 km 2 yr −1 , and a median of 0 km 2 yr −1 , therefore also underestimating observations.

Ecoregions ART, CR and AM
For ecoregion ART (Alaska Range Transition) LPJ-LMfire simulates a mean annual burned area of 134 ± 393 km 2 yr −1 and a median of 4 km 2 yr −1 compared to an observed mean annual burned area of 91 ± 109 km 2 yr −1 and a median of 37 km 2 yr −1 (Fig. 6, Table 4). We therefore underestimate the median while overestimating the mean, with the latter again being augmented due to one single fire year, 2007, for which we simulate a maximum of 1907 km 2 yr −1 against an observation value of only 299 km 2 yr −1 . All other 24 yr for ecoregion ART are within the range of observation concerning total area burned and interannual variability. Ecoregions CR (Coastal Rainforest) and AM (Aleutian Meadows) are ecoregions with extremely low amounts of burned area, both observed and simulated, in total as well as percentage of region's area. For ecoregion CR, an average of 13 ± 38 km 2 yr −1 in the observation data compares to a simulated average of 10 ± 47 km 2 yr −1 . In ecoregion AM, burned area is recorded in 4 out of the 25 yr of observation compared to 2 yr of fire simulated by LPJ-LMfire. These results reveal that though we may not be able to reproduce exact numbers for area burned at the very low end of fire observations, we are still able to simulate fire occurrence behaviour realistically even in areas where burning is rare and reproducing any fire at all in the simulations is challenging.

Discussion of Alaska burned area results
While overall mean simulated burned area was close to that observed, peak fire years in our simulated time series did not always match observed peak fire years (Fig. 7). The cause for this mismatch may be linked to the uncertainty in daily weather conditions resulting from the usage of a weather generator and monthly climate data. Using monthly climate forcing constrains total precipitation amount and number of wet days, but the timing of rainy days within a given month may be very different in the simulation compared to the true weather situation, e.g. if simulated wet days all come clustered at the beginning or end of the month, whereas in reality they had been more equally distributed over the month. In such a case, the consequences for fuel wetting and drying are different between observation and simulation, with simulation overestimating fuel dryness and FDI and therefore leading to higher amounts of area burned. Moreover, the timing and amount of precipitation matters for simulating fire extinction in LPJ-LMfire, as either one day with more than 10 mm precipitation (3 mm precipitation with more than 60 % grass cover) or several consecutive days with a sum of more than 10 mm precipitation are required to extinguish fires in our simulation. If, for example, a fire is burning in a given month and the simulated clustering of rainy days within this month is less pronounced than the clustering that occurred in reality, the fire may continue burning although in reality it was extinguished. This may also be true for the opposite case, where fires are extinguished although they should have kept burning. Another uncertainty is linked to wind speed: as we lack the capability in our weather generator to disaggregate wind speed to daily or hourly values, we use climatological mean wind speed, which may underestimate the infrequent, high-wind events that are responsible for the largest episodes of fire spread. Finally, LPJ-LMfire does not simulate the feedback mechanism between fire and wind; for example, large, intense fires such as those observed in boreal forests may produce strong convection that increases wind speeds in the vicinity of the fire, which in turn enhances fire spread.
Correct simulation of fires in tundra regions is challenging for several reasons. The most significant problem leading to a general overestimation of simulated burned area on the Alaska North Slope is the simple soil water scheme of LPJ that is not able to explicitly simulate permafrost or wetlands. Detailed analyses of grid pixels in northern Alaska revealed that soils dry out very quickly as soon as all snow has melted in May or beginning of June, and, because it is linked to soil moisture, the water content of the live grass drops quickly. Summers in northern Alaska are dry, while at the same time day length is long; therefore simulated evapotranspiration is high and helps to draw down soil moisture in combination with surface runoff and drainage. Overall, this leads to simulation of environmental conditions that are far drier than in reality where thawing of the active layer proceeds slowly down the soil column over the course of the summer and, by limiting evapotranspiration, keeps soils and vegetation wetter than would otherwise be the case. If lightning occurs in the period between May and July, simulated fires spread very fast and therefore lead to an overestimation of burned area. In most of the cases where we overestimate burning, fires are ignited early in summer when in reality conditions are likely still too wet; the simulated fires spread quickly due to the fuel being dry and keep burning through summer due to the lack of precipitation. In addition to the poor representation of wetlands and permafrost in LPJ, the tundra on Alaska North Slope is characterized by a high density of water bodies including many lakes, peatlands, streams and rivers, which is not taken into account in LPJ. In reality, these water bodies will limit the spread of fires, as can be observed for the Anaktuvuk River fire which is bordered by rivers on its western and eastern margins. Future improvements to LPJ and the fire model therefore should focus on the implementation of adequate permafrost and wetland simulation modules (e.g. Wania et al., 2009;Koven et al., 2009;Ringeval et al., 2010) and the incorporation of some spatial statistic representing water body distribution on a grid cell level as a limiting factor to the spread of fires. This could be accomplished similarly to the way in which we account for the effects of landscape fragmentation on fire size as a result of topography (Sect. 3.2.5) or land use (Sect. 3.2.6). As LPJ-LMfire has no PFT that specifically represents it, tundra vegetation in the model is simulated with the C 3 -grass PFT. As described in Sect. 3.2.2, we tried to improve the representation of tundra vegetation with respect to fuel conditions by scaling the density of live grasses to the number of growing degree-days and by accounting for permafrost-impeded drainage of soil water. Eventually, woody shrub vegetation and tussocks could be represented by one or more separate tundra PFTs (e.g. Kaplan et al., 2003;Wania et al., 2009) as each of the constituent tundra vegetation plants have different density, height, and flammability that would affect fire spread.
Comparing the Bering Taiga and Bering Tundra ecoregion to the Arctic Tundra in northern Alaska reveals that all three ecoregions are characterized by generally very low amounts of lightning. They can therefore all be classified as ignition-limited fire regimes. In contrast to the Arctic Tundra region, the two western regions have their precipitation maximum in summer, which coincides with the potential fire season. As a consequence of frequent rainfall events with oftensubstantial daily precipitation amounts, fuels stay wet and soil water status is high (Fig. 8). In the already rare case of a lightning ignition, fires therefore tend to spread slowly, stay small and are soon extinguished, especially when compared to fires started in the Arctic Tundra. Rare but important fires in boreal and subarctic environments develop during particular conditions, e.g. an exceptionally long string of dry weather. As LPJ-LMfire uses a weather generator to disaggregate monthly climate variables to daily values, it is possible that the specific circumstances that in reality led to a fire, i.e. having an ignition while at the same time simulating a sufficiently long dry period after the ignition so that the fire can spread, are not captured by the model simulation. With only few lightning sensors located in the far west of Alaska, it is also possible that the actual amount of lighting occurring in these two ecoregions is underestimated and not all lighting is recorded.
Apart from the limitations discussed here, using daily and interannually variable lightning as described in Sect. 3.1.2 allows us to simulate fire in boreal regions, with results showing considerable interannual variability in total burned area. Although we may not be able to reproduce observed annual area burned exactly on a year-to-year basis because of the limitations highlighted above, with LPJ-LMfire we capture the overall behaviour of boreal fires well in terms of being able to simulate long-term averages and variability that are consistent with observations.

Simulated fire return intervals in Alaska
Fire return interval (FRI), i.e. the number of years between successive fires in an area, is widely used to characterize natural fire regimes and assess the changes in fire frequency caused by climate change. For the recent past, efforts to reconstruct FRIs based on fire scar data sets have been performed by Balshi et al. (2007), who present maps of fire return intervals in boreal North America and Eurasia using historical fire records for the second half of the 20th century. In places where fire is infrequent, however, FRIs may ex-12 25 50 100 200 300 400 500 700 1000 2000 fire return interval (years) Fig. 9. Simulated fire return intervals in Alaska for a 1000 yr run with detrended 20th century climate. To facilitate comparison, the colour schemes used here and in Fig. 11 are the same as those used in Balshi et al. (2007). ceed the period of modern observations. Detailed historical records of burned area in the boreal forest in the best case hold a little more than 70 yr of data in Alaska and Canada and even less than that in Eurasia. Short records may be not representative of the overall average fire regime as by chance they may, for example, represent a time of relatively high or low fire activity and therefore lead to an overestimation or underestimation of average FRIs over longer time scales. The need to perform spatial interpolation of FRIs over large spatial scales introduces further uncertainty.
Analysis of charcoal accumulation rates from sedimentary archives has been applied successfully on local to regional scales to reconstruct FRIs over longer time scales (e.g. Higuera et al., 2009;Lynch et al., 2004;Tinner et al., 2006;Higuera et al., 2008;Brubaker et al., 2009). However, centennial to millennial scale climate variability probably affected FRIs as ecosystems adjusted to changing climate. It is therefore difficult to characterize steady-state equilibrium FRIs or estimate how future climate changes could affect burning, based solely on palaeo-archives. The advantage of DGVMs containing fire models is that they can be run for long time periods using detrended steady-state climate, allowing vegetation and fire regime to equilibrate so that conclusions can be made as to what the equilibrium FRI would be if climate at any given time stayed constant.
To estimate FRIs for Alaska, we made a model run over 1000 yr with steady-state climate after vegetation and fire regime had equilibrated. Following Balshi et al. (2007), we define FRI as the time required to burn an area equal to the entire 0.5 • grid cell. The FRI within a grid cell is consequently calculated as the ratio of 1000 yr and the number of times a grid cell area burned during these 1000 yr. We present our simulated fire return intervals in Fig. 9, using the same colour scheme as in Balshi et al. (2007), but without applying any smoothing. Agreeing with Balshi et al. (2007), we simulate frequent burning with return intervals between 12 and 50 yr in eastern Alaska located in the Intermontane Boreal ecoregion between Brooks Range and Alaska Range. Towards the west of ecoregion IB, the FRIs predicted from our simulation become more heterogeneous, from less than 50 yr to more than 500, therefore being slightly lower than the FRIs estimated by Balshi et al. (2007). Towards the extreme west of mainland Alaska, we simulate FRIs between 900 and 2000 yr for some grid cells, but mostly FRIs are longer than 2000 yr. Compared to Balshi et al. (2007) we estimate significantly longer FRIs in some grid cells, especially for ecoregion BTU (Bering Tundra). This may be linked to the possibility that the already low amounts of lightning are underestimated in the LIS/OTD lightning climatology used for this experiment, due to the limited 4 yr length of record of the lightning climatology and the low detection efficiency at high latitudes. In contrast, we simulate shorter fire return intervals for the Arctic Tundra, which typically fall in the 100-200 yr and 500-700 yr categories. Given the model shortcomings related to the simulation of tundra vegetation and permafrost (see Sect. 4.3.2), these results may be biased somewhat towards shorter FRIs than are actually observed.

Global fire under natural conditions
To characterize the behaviour of LPJ-LMfire globally and place it in the context of previous fire modelling work, we performed an experiment analogous to that presented by Bond et al. (2005), contrasting global biomass in a "world without fire" to one where natural fires are simulated. The global effects of fire on aboveground live biomass are shown in Fig. 10. Both panels represent a world with potential natural vegetation and no anthropogenic land use. Panel (a) shows biomass with natural fires caused by lightning ignitions, while panel (b) shows a world without fire. Panel (c) shows the difference in biomass between a world with and without fire. The maps clearly reveal the parts of the world that are mostly affected by fire disturbance and therefore have less biomass than they potentially could have in a world without fire. On a 100 yr basis, the total amount of global carbon stored in aboveground living biomass is 208 ± 2 Pg less for the simulation with fire compared to the simulation without fire, totaling 948 ± 3 Pg C with fire. No impact of fire on biomass is simulated for the wet tropics where very little fire is simulated, such as the Amazon and Congo basins or in Indonesia, all places that naturally store large amounts of carbon in forests. Most of the biomass loss related to fire disturbance is simulated in the seasonal tropics and subtropics: in the Miombo woodland region south of the Congo Basin, in the east and southeast of the Amazon Basin, in the Sahel, in India and Southeast Asia, and in northern and southern Australia. The impact of fire on biomass is also clearly visible in the grassland regions of central and western North America, the western Mediterranean, southwestern Russia, Kazakhstan and Uzbekistan. Fires in the boreal regions can be extensive, but the return interval is too long to have a discernible impact on carbon storage in aboveground live biomass compared to ecosystems with short fire return intervals.
The results we present here are broadly consistent with those in Bond et al. (2005), who showed, in a series of experiments running a DGVM with and without fire, that the largest reductions in tree cover as a result of natural fire are in the seasonal subtropics. Bond et al. (2005, Fig . 6) also show a large reduction in forest cover in central Europe and the eastern United States, areas where fire impacts in LPJ-LMfire are more muted. In contrast, LPJ-LMfire shows a large reduction in biomass in the grassland areas of central North America, on the Eurasian steppe, in central and southern Australia, and in southern South America when comparing "fire on" with "fire off" scenarios. Bond et al. (2005) state that FRIs simulated by their model in these natural grassland areas are much too long with respect to observations (75-200 yr modelled where 2-5 yr are observed). LPJ-LMfire shows much shorter FRIs (Fig. 11) of 1-5 yr in much of these natural grassland regions that are more consistent with field observations. The map of global FRIs in Fig. 11 shows that fires are most frequent in places where three factors are coincident: a. enough biomass to sustain frequent burning; b. sufficient amounts of lightning ignitions; c. seasonally varying meteorological conditions, specifically a pronounced dry season that allows fuel drying.
If any of these three conditions is not present, wildfires are unlikely to occur. As noted above, fire is rare in the Amazon and Congo basins and on the Indonesian archipelago. In these regions, lightning ignitions and biomass are not limiting, but meteorological conditions are typically too wet for the development of wildfires, with the exception of relatively infrequent severe drought events, e.g. in extreme El Niño years (Page et al., 2002(Page et al., , 2012. In the desert and high-mountain  regions of the world, e.g. in the Sahara desert, the southern part of the Arabian Peninsula, and on the Tibetan Plateau, the absence of biomass is the limiting factor for fire. Large parts of the world's boreal and subarctic ecosystems have enough biomass to support frequent burning, but the number of lightning ignitions generally tends to be low compared to lower latitudes, with snow and temperatures below 0 • C occurring for half a year or more, and the summer season is frequently the wettest time of the year.
In contrast, in any part of the world where all three factors are met, fire return intervals are short, e.g. in the Sahel, the western Mediterranean, the Near East, in the Miombo woodlands south and east of the Congo Basin, in most of Australia, and in the xerophytic Caatinga shrublands of northeastern Brazil.

Comparison to contemporary observations of burned area
While LPJ-LMfire has been primarily designed to simulate fire behaviour during preindustrial time, we compared the results of a global model run with satellite-based estimates of burned area that cover recent decades. In our model experiments we did not attempt to account for either anthropogenic ignitions or active suppression of wildfires, but we did account for passive fire suppression through landscape fragmentation as a result of agricultural land use. The differences between simulated and observed burned area may therefore, in certain regions, highlight the importance of human influence on the geographic distribution of fire at present. In a few parts of the world where human impact is minimal, we were further able to identify potential shortcomings of the current version of LPJ-LMfire and priorities for future model development.
As described in Sect. 3.4 above, we ran LPJ-LMfire with climate and soils data that reflect the late 20th and early 21st centuries (Table 3). The model was spun up for 1020 yr with 1871 CO 2 concentrations and land use, and then run in a transient climate, CO 2 , and land use scenario for the period 1871-2010. Used land was defined as the sum of the agricultural and urban fractions and was specified from the HYDE v3.1 anthropogenic land cover change scenario (Klein Goldewijk et al., 2010). In our simulations, fires were only allowed to burn on the unused fraction of each grid cell, and the only ignition source was lightning. We compare our model results with the global burned area products GFEDv3.1 (Giglio et al., 2010, hereafter GFED) and the data set presented by Randerson et al. (2012, hereafter JR12). GFED provides complete annual coverage for the years 1997-2011, while JR12 covers the period 2001-2010. The main difference between the two observational burned area products is that JR12 accounts for numerous additional small fires not included in GFED, which results in an increase in mean annual burned area of up to 30 % in some regions, mainly in the tropics and subtropics.
We compare modelled with observed burned area on the basis of a multi-year mean of the annual total burned area fraction of each 0.5 • grid cell. We extracted the time periods from our LPJ-LMfire run overlapping with the period covered by the observational data sets, summed the monthly values in the observational data sets to create annual totals, and calculated average burned area over the number of years of record. In comparing LPJ-LMfire with GFED, we masked the difference between model and observation where the differences were less than the aggregate uncertainty specified in the GFED database. For comparison with JR12, we masked areas where the model-data mismatch was less than 1 %.  The differences between LPJ-LMfire and GFED are shown in panel a of Fig. 12; differences with JR12 are in Fig. S8. Overall, the spatial pattern and magnitude of the residual between model and observations are similar regardless of the observational data set we used. The greatest differences between model and observations are found in the seasonal tropics of Africa, both north and south of the Equator, where LPJ-LMfire shows substantially less burned area than the observations. Further large negative residuals are seen in northern Australia, along the steppe belt of Eurasia from Ukraine to Kazakhstan, in Southeast Asia particularly in Cambodia, in the Amur region of the Russian Far East, and in the lowlands of Bolivia and Paraguay. In contrast, the model shows relatively more burned area compared to observations in several regions, notably in the Caatinga region of north-eastern Brazil, in Iran and western Turkmenistan, in most of southern Australia, in the western United States, and in the Chaco dry forest region of northwestern Argentina.
In panel b of Figs. 12 and S8, we place these differences between model and observations in the context of the anthropogenic imprint on the global land surface by means of a simple classification of the residual based on human impact. We specified human impact based on the GLOBIO methodology (Ahlenius, 2005, Fig. S9), which identifies the presence of anthropogenic features on the ground including urban areas, open cast mines, airports, roads, railroads, canals, and utility lines. Half-degree grid cells covered 1 % or more by anthropogenic features were classified as being substantially influenced by human activities. On the basis of this classification, 75 % (347 out of 464 Mha) of the mean annual global burned area in JR12 occurs on land influenced by human impact (for GFED the amount is 71 % or 258 of 363 Mha). The areas of largest disagreement between model and observations are even more concentrated in regions of human impact. In grid cells where the difference between LPJ-LMfire and the observational data sets is greater than 1 %, human impact is visible on 93 % of the area.
As discussed above, where LPJ-LMfire underestimates observed burned area in much of the seasonal tropics, this is most likely the result of regular, intentional burning of cropland and pastures, and in some cases deforestation fires. In the Sahel, extensive agricultural burning is common practice and occurs annually for several months during the Northern Hemisphere autumn and winter when people ignite fires to remove crop residues, and to renew pasture grasses and stop woody encroachment on pastures, for hunting, and to control pests and wildfires (Menaut et al., 1991;Klop and Prins, 2008;Kull and Laris, 2009;NASA, 2011). This is also the case for the Miombo woodlands south of the Congo Basin and on Madagascar, where intentional burning plays an important role (Eriksen, 2007;Le Page et al., 2010). Van Wilgen et al. (1990) estimate that humans cause 70 % of all annual fires in African savannas. Likewise, the large underestimate in burned area in the Eurasian steppe and Amur valley is related to agricultural burning (Tansey et al., 2004;Warneke et al., 2009).
In places where LPJ-LMfire overestimates burned area relative to observations and human impact is considered important, three processes that are not included in LPJ-LMfire may explain the differences: (1) removal of biomass from grazing and land degradation, (2) industrial fire suppression, and (3) landscape fragmentation from roads and other anthropogenic features that are not classified as agricultural land use by HYDE. Eastern Brazil, northern Argentina, southern Australia, East Africa, northern Mexico, and the western Great Plains of the United States are occupied by extensive open rangelands (Klein Goldewijk et al., 2010;Ramankutty et al., 2008) where livestock grazing leads to a reduction of fine fuel load that is not accounted for in our model simulations. Furthermore, the semiarid regions of the northern Sahel, Central Asia, and the Near East are characterized by extensive soil degradation as a result of overgrazing and millennia of land use (Dregne, 2002). This degradation is not accounted for in our simulations, and the model simulates more biomass than is actually present and, therefore, more fire. Industrial fire suppression is common in Europe, North America and Australia, and may explain much of the additional LPJ-LMfire overestimate relative to observations in these areas. In the US alone, expenditure on wildfire suppression has increased continually over the last 70 yr (Stephens and Ruth, 2005;Calkin et al., 2005;Westerling et al., 2006;Nazzaro, 2007;Gebert et al., 2007Gebert et al., , 2008. Finally, while we accounted for landscape fragmentation as a result of agricultural land use in our simulations, additional fragmentation effects caused by the presence of human infrastructure such as roads were not included. The combination of industrial fire suppression with a high magnitude of human impact is the likely cause for the overestimate in burned area in developed countries of the temperate regions.
On the remaining 7 % of land area that may be classified as having insignificant human impact, we show overestimates in burned area in subarctic western Canada and eastern Siberia, in a small area along the southern margin of the Sahara in Mali and Niger, and markedly in the southeastern Amazon Basin in the transition zone between tropical forests and the Cerrado savanna. In contrast, the model underestimates fire in boreal Canada, in the eastern Central African Republic, central Australia and in central Brazil. These residuals are useful for understanding the limitations of both our model and the observational data sets.
The unprojected maps in Figs. 12 and S8 exaggerate the area, but the overestimate in burned area in the subarctic may be caused by several factors that were already discussed in our analysis of LPJ-LMfire results for Alaska. These include an inadequate representation of permafrost that influences soil hydrology and therefore fuel moisture, an overall overestimate in modelled aboveground biomass also caused by permafrost and/or lack of soil, and fine-scale landscape fragmentation caused by the rivers, lakes, and wetlands that are extensive in this region . Permafrost is important in all of these northern areas where LPJ-LMfire overestimates observations . In areas of boreal Canada further south where LPJ-LMfire underestimates burned area, several factors not included in our simulations may influence the model results, including tree kill events as a result of insect infestations, and remote industrial activities including logging, and hydroelectric and oil and gas development. Furthermore, in the boreal and subarctic areas where modelled burned area is both under-and overestimated, the very short 4 yr period upon which our lighting climatology is based for the extratropics means that we may misestimate the number of ignitions in these regions. The large interannual variability and differences in spatial pattern in lightning we observed using the ALDS data for Alaska shows that the LIS/OTD climatology is at best a rough approximation of the actual amount of lighting strikes, and that, in certain years, our LIS/OTD-based estimates could result in a substantial over-or underestimate in the actual number of potential ignition events.
As described above, most of the temperate regions of the world are so extensively impacted by human activities, both at present and historically, that it is impossible to disentangle the natural fire regime from anthropogenic influences at the 0.5 • spatial resolution used in our model simulations. A better test of our process fire model in mid-latitude settings could be to perform case studies in protected regions such as national parks or wilderness areas at very high (∼ 1 km) spatial resolution.
In the subtropics and tropics, the largest area of disagreement between model and observations is in the transition zone between Cerrado savannas and tropical forests in the southeastern Amazon in northern Mato Grosso and southern Pará states of Brazil. In this region, a pronounced dry season of three to four months combined with high temperatures leads to rapid fuel drying and high FDI in our simulations. Combined with lightning activity in all seasons as seen in the LIS/OTD climatology, modelled fire frequencies and burned area in this region are high. In reality, several factors preclude the development of large fires in this region and should be included in future improvements to the model. In this region tropical forests tend to develop on lowland environments with deep soils; the maximum rooting depth in LPJ-LMfire is 2 m, but much more deeply rooted trees have been commonly observed in the Amazon (Kleidon and Heimann, 2000). More deeply rooted trees would extend the period of greenness for tropical raingreen vegetation, the dominant PFT in this region in our model simulations, effectively limiting the length of the dry season for the vegetation.
Furthermore, green forest vegetation will effectively shade the litter and other fuel in the understorey, reducing the rate at which it will dry out. In experiments in the seasonal tropics of the southeastern Amazon, Uhl and Kauffman (1990) showed that land cover controls fuel moisture, and that under equivalent climate conditions, fuels in intact forests would never become dry enough to burn, while in grasslands only 24 h of dry weather was required to support sustained burning. In this sense, the Nesterov Index approach used for estimating fuel moisture in SPITFIRE may be inadequate. We suggest that future models would benefit from an energy balance approach to estimating fuel moisture, particularly in forest understoreys. On the other hand, it is possible that the observations of burned area in this region are underestimates of the actual situation. Randerson et al. (2012) suggest that burned area in their data set may be particularly underestimated in regions where small fires occur in forest understorey. The combination of these limitations in both the model and the data sets probably leads to the large positive residual observed in this region.
Adjacent to this region of overestimated modelled burned area in central Brazil is a discontinuous region of underestimated burning in areas shown on our map to be largely free of human influence. This region in Mato Grosso, Goiás, and Tocantins states is an area where rapid land cover change, in the form of deforestation and conversion to agriculture and pasture, has been important in recent decades (de Souza et al., 2013). This recent deforestation has been documented as being associated with an increase in burned area (Lima et al., 2012). Our human influence map is based on the VMAP0 data product (NIMA, 2000) that was largely assembled from data collected during the period 1974-1994. We suggest that the negative residual in burned area in this region is a result of recent human activities not currently captured by our human impact database.
Similarly, the large areas of underestimated burned area in the easternmost Central African Republic and northern Australia have been attributed to human action, though not as a result of anthropogenic land cover change. In Africa, large savanna fires are intentionally lit to facilitate hunting in sparsely populated areas (Eva et al., 1998). Likewise, in northern Australia frequent intentional human burning is an important part of traditional landscape management that is widespread in sparsely populated areas at present (McKeon et al., 1990;Dyer, 1999;Yibarbuk et al., 2002;Crowley and Garnett, 2000).
In summary, our simulations of burned area over recent decades caused only by lightning ignitions and only passively suppressed through agricultural and urban land use show substantial differences with observational data sets of fire. We expect these differences, because the complex human relationship with fire at present is well known, and we have made no attempt to prescribe this in our simulations. In parts of the world where human impact is limited, modelled mean burned area fraction often agrees within 10 % of observations on a decadal average. In those areas of low human impact where we do show important disagreement between model and observations, we can identify limitations in our model and driver data sets. In boreal and subarctic Canada and Russia, we may overestimate fire because of our over-simplistic treatment of permafrost, and because of the presence of lakes, wetlands, and barren ground that are not accounted for in our model input. In the tropics and subtropics, we may overestimate burning because of an inadequate representation of the effects that canopy shading and deeply rooted vegetation have on fuel moisture. Future developments to the model should address these issues by improving the soil hydrology scheme, by using recently developed methods for simulating permafrost (Wania et al., 2009) or deep tropical soils (Poulter et al., 2009).

General discussion
Realistic simulation of global vegetation dynamics requires the inclusion of disturbance regimes that influence vegetation development, alter vegetation structure and composition and affect global carbon budgets. Simulation of fire, arguably the most important disturbance process that affects the terrestrial biosphere, is of crucial importance for a complete model representation of terrestrial vegetation dynamics. Starting with SPITFIRE, we developed LPJ-LMfire to overcome some of the shortcomings of the original model. LPJ-LMfire includes major changes to the process representations of fire occurrence, fire spread and fire impact. In boreal and subarctic regions in particular, LPJ-LMfire results are in much better agreement with observations compared to SPITFIRE or LPX. In other parts of the world, the changes that we made to SPITFIRE in developing LPJ-LMfire are harder to distinguish at the coarse resolution at which we run the model because of the pervasive nature of human impact on fire at present, both through ignitions and fire suppression.
Under a natural fire regime excluding human interference, lightning is the most common ignition source for wildfires. Accounting for interannual variability in lightning and the occurrence of lightning on a daily timescale is important, especially in regions where the total amount of monthly lightning strikes is low, and therefore an equal distribution of lightning strikes on all days within a month may result in a significant underestimate of lightning ignitions. By correlating the occurrence of lightning strikes with the occurrence of precipitation, we provide a more realistic way to simulate lightning ignitions. In boreal and subarctic environments where lightning ignitions are rare, SPITFIRE and LPX show unrealistically little fire compared to observations. With meteorological forcing that is similar to that used by these models, but by accounting for intra-monthly and interannual variability in lighting ignitions in LPJ-LMfire, we show that simulation of realistic fire behaviour is possible. This is illustrated by our time series of burned area for central Alaska (Fig. 7) that includes single years with significant amounts of area burned, while many other years have no or only very little fire. It is possible that the meteorological forcing used by all three models is unrealistic, but the inclusion of observed variability in lightning strikes in the model provides a parsimonious solution that results in a better match with observations.
By allowing the ignition of smoldering fires during wet conditions and simulating fires that persist over the course of multiple days instead of extinguishing each fire after a length of time that is limited to 241 min, LPJ-LMfire more closely reflects the true behaviour of fire. Likewise, the calculation of fuel wetness as a mass balance function of wetting and drying in LPJ-LMfire, rather than relying on a precipitation threshold of 3 mm as suggested by Nesterov (1949) and used in SPITFIRE, made a substantial improvement in the agreement of the model results with observations. Eventually, the introduction of additional shrub PFTs as intermediates between herbaceous vegetation and tree PFTs should be considered, especially for an appropriate representation of tundra and xerophytic vegetation. Introduction of shrub PFTs will help ameliorate the current tendency of the model to overestimate herbaceous vegetation cover in fireprone areas and the strong positive feedback between fire and vegetation that results in an overestimate of fire frequency and the prevalence of grasses, a problem sometimes still observed, for example, in the Arctic tundra of northern Alaska, or in southern Spain and central Australia. Further improvements should also focus on the inclusion of a scheme to simulate wetlands and permafrost in order to capture the way in which permafrost keeps tundra organic matter wet, even under dry meteorological conditions. Since our version of LPJ does not represent permafrost dynamics, soil and fuel drying and hence fire occurrence are overestimated in tundra areas such as northern Alaska where wetlands and permafrost are common. Other future improvements to LPJ-LMfire should include development of a scheme to simulate crown fires in addition to the surface fires simulated by the current version of the model.
By introducing a slope factor related to the median slope angle of each 0.5 • grid cell, we present a simple way to account for the role that topographic complexity plays in limiting fire size and rate of spread. Eventually, a representation of other natural firebreaks such as rivers and lakes should be built into the fire module. An approximation of the number of rivers that could act as fire breaks could be handled by using drainage density information extracted from a digital elevation model. That rivers constrain the spread of fires can be observed, for example, in case of the large Anaktuvuk River fire from 2007 in the Alaskan tundra that was ultimately constrained by the two rivers: Nanushuk to the west and Itkillik to the east (Jones et al., 2009). A measure of fragmentation by water bodies could be indirectly accounted for using Eq. (14), which links the numbers of fires burning at any time to the degree to which the landscape has been fragmented due to previous burns in the fire season.
Grass PFTs should be implemented such that they are able to reach full cover and complete their life cycle within one growing season. Our overall simulation results indicate that this would be particularly important for mesic tropical savannas, where fire is a prevalent feature of the ecosystem and most species of grass have annual lifecycles (Scholes et al., 1997). To accomplish this it would be necessary to run the entire model at a monthly or shorter timestep. Calculating processes such as allocation, turnover, and mortality that are currently updated annually in LPJ on a shorter timestep would also provide the additional advantage that burned area could be tracked continuously over time rather than resetting calculated burned area at the beginning of each calendar year. While the Northern Hemisphere summer fortunately roughly corresponds with the fire season for a large part of the world, it is not correct for southern South America, southern Australia, and parts of Southeast Asia and the Sahel.
Our comparison of the LPJ-LMfire results with observational data sets of burned area shows that anthropogenic impact on controlling the spatial pattern of fire observed in recent decades may be important in many parts of the world. About three-quarters of the earth's surface where fire is observed at present occurs on lands influenced by human activities, and human interactions with fire range from rigorous suppression, e.g. in the US, Europe, or parts of Australia, to liberal intentional burning on agricultural and natural land, e.g. in Sub-Saharan Africa. The sheer variety of examples for human-fire interactions makes it clear that modelling the spatial and temporal pattern of fire at the present day would require detailed parameterizations of the human relationships with fire at sub-national level. Prescription of human behaviour with respect to fire could likely overcome most of the differences between observed and simulated results, but would no longer be process-based fire modelling. Moreover, such a detailed parameterization of anthropogenic behaviour would reflect present-day customs and policies with respect to fire, but would not have predictive power for the future or for the past, as human preferences and subsistence strategies change over time. Our main interest in developing LPJ-LMfire was to model fire, including anthropogenic fire, in the preindustrial past, when humans had neither present-day technology to suppress and control wildfires, nor modern-day agricultural technologies that allow people to abstain from using fire for agricultural purposes.
We therefore decided to develop a more general scheme for representing human-fire interactions in the past that is based on observations and knowledge on people who still use fire in traditional ways at present day, such as Australian Aborigines or subsistence farmers in developing countries, or on historical ethnographic observations of fire usage, e.g. of the native North Americans. Literature research, evidence from palaeoproxies such as charcoal preserved in sediments, and discussions with anthropologists and archaeologists led us to the conclusion that humans in the past used fire for a variety of different reasons, depending on their lifestyles and habitat, and that terrestrial biomass burning related to human activity must have been very common. By developing a method of representing the way in which people with different subsistence lifestyles interact with fire, with LPJ-LMfire we are able to perform quantitative estimates on the impact of anthropogenic burning on vegetation, carbon pools and trace gas emissions on a global scale during preindustrial time. We realize that this approach may be too simplistic to address specific local-scale peculiarities of human burning behaviour, but believe nonetheless that our approach to classifying people's relationship with fire based on their subsistence lifestyle is a more appropriate way of addressing the implementation of human burning in the past than prescribing the patterns of human-influenced fire observed in the 21st century.

Conclusions
Beginning with LPJ-SPITFIRE (Thonicke et al., 2010), we made improvements to several aspects of the original formulation and achieved a more realistic process representation of fire occurrence, fire behaviour, and fire impacts, particularly in boreal and subarctic ecosystems. With our updated model, LPJ-LMfire, we were able to simulate realistic fire regimes in Alaska, one of the key regions of the world where SPITFIRE results did not agree with observations. We also developed a scheme to distinguish among the ways in which preindustrial people with different subsistence strategies interact with fire to achieve their land management goals. LPJ-LMfire is a major improvement on past global fire models and will be particularly useful for studying changes in global fire on millennial timescales, providing a basis for further improvements, modifications and model development.

Appendix A
In this appendix we provide the equations used in LPJ-LMfire that were not changed from the original SPITFIRE. With these, we provide a complete documentation of LPJ-LMfire. Variable and parameter abbreviations used in addition to those in Table 1 are provided in Table A2.
Total dead fuel mass within the first three fuel size classes, plus mass of the live grass:
Reaction intensity: Ratio of propagating flux to reaction intensity: Wind coefficient: Heat of pre-ignition: The rate of spread is then calculated as follows: Backward rate of spread (decreases with increasing wind speed): ROSb s = ROSf s · e −0.012·U forward . (A37)
The total distance travelled by a fire within a day is estimated as DT = t fire · (ROSf + ROSb) .
The mean area burned by one single fire is calculated as a f = min π 4 · LB · DT 2 · 0.0001, ac area .

A4 Combustion of dead fuel
Fraction of live grass consumed by surface fire: To calculate how much fuel has been consumed in total within one grid cell over the course of a year, FC(class) needs to be multiplied with the annual area burned (in m −2 ).
Calculation of surface fire intensity: If the surface fire intensity is less than 50 kW m −1 , it is considered to be too low for burning and fires are extinguished.

A5 Fire mortality and combustion of live fuel
Crown scorch is calculated per PFT. For seasonally leafbearing trees, crown scorch is relevant as long as there are leaves on the tree.
The probability of mortality due to cambial damage is given by where τ l /τ c is the ratio of the residence time of the fire to the critical time for cambial damage (Peterson and Ryan, 1986). The critical time for cambial damage τ c (min) depends on the bark thickness (BT) (cm): τ c = 2.9 · BT 2 , (A61) (Peterson and Ryan, 1986;Johnson, 1992), which is calculated from the diameter at breast height (DBH, cm) using where par1 and par2 are PFT-specific constants (Table A1).
The total probability of mortality due to crown damage P mCK and cambial damage P m (τ ) is calculated as P m = P m (τ ) + P mCK − P m (τ ) · P mCK . (A63)
These are calculated on a daily basis. To calculate the annual total, the daily sum is accumulated over the course of the year: annBB dead(PFT) = annBB dead(PFT) + BB dead (PFT) (A69) Biomass is burned from live fuel by fuel type and PFT. For tree-type PFTs: The annual running sum of mortality probability is calculated as ann kill(PFT) = ann kill(PFT) + P mCK(PFT) · AB frac . (A75) Updating of the litter pools is done once at the end of the year: laf(PFT) = max laf(PFT) − annBB dead(PFT,1) , 0 , For the tree-type PFTs, live biomass that was killed but not consumed by burning is transferred to the litter pools, and the individual density is updated based on the fraction of individuals that were killed over the course of the year: N ind−kill(PFT) = ann kill(PFT) · N ind(PFT) , In case of a PFT being killed off completely by fire, reset presence to "false" and set all biomass pools of that PFT (lm ind(PFT) , sm ind(PFT) , hm ind(PFT) , rm ind(PFT) ) to zero.

A7 Trace gas emissions
Total carbon emissions from burning, across all PFTs: