A fire model with distinct crop , pasture , and non-agricultural burning : Use of new data and a model-fitting algorithm for FINAL . 1

This study describes and evaluates the Fire Including Natural & Agricultural Lands model (FINAL) which, for the first time, explicitly simulates cropland and pasture management fires separately from non-agricultural fires. The nonagricultural fire module uses empirical relationships to simulate burned area in a quasi-mechanistic framework, similar to past fire modeling efforts, but with a novel optimization method that improves the fidelity of simulated fire patterns to new observational estimates of non-agricultural burning. The agricultural fire components are forced with estimates of cropland 5 and pasture fire seasonality and frequency derived from observational land-cover and satellite fire datasets. FINAL accurately simulates the amount, distribution, and seasonal timing of burned cropland and pasture over 2001–2009 (global totals: 0.434× 10 and 2.02× 10 km yr−1 modeled, 0.454× 10 and 2.04× 10 km yr−1 observed), but carbon emissions for cropland and pasture fire are overestimated (global totals: 0.295 PgC yr−1 and 0.706 PgC yr−1 modeled, 0.194 PgC yr−1 and 0.538 PgC yr−1 observed). The non-agricultural fire module underestimates global burned area (1.91× 10 km yr−1 modeled, 10 2.44×10 km yr−1 observed) and carbon emissions (1.14 PgC yr−1 modeled, 1.84 PgC yr−1 observed). The spatial pattern of total burned area and carbon emissions is generally well reproduced across much of sub-Saharan Africa, Brazil, central Asia, and Australia, whereas the boreal zone sees underestimates. FINAL represents an important step in the development of global fire models, and offers a strategy for fire models to consider human-driven fire regimes on cultivated lands. At the regional scale, simulations would benefit from refinements in the parameterizations and improved optimization datasets. We include an 15 in-depth discussion of the lessons learned from using the Levenberg-Marquardt algorithm in an interactive optimization for a dynamic global vegetation model.

Abstract.This study describes and evaluates the Fire Including Natural & Agricultural Lands model (FINAL) which, for the first time, explicitly simulates cropland and pasture management fires separately from non-agricultural fires.The non-agricultural fire module uses empirical relationships to simulate burned area in a quasi-mechanistic framework, similar to past fire modeling efforts, but with a novel optimization method that improves the fidelity of simulated fire patterns to new observational estimates of non-agricultural burning.The agricultural fire components are forced with estimates of cropland and pasture fire seasonality and frequency derived from observational land cover and satellite fire datasets.FINAL accurately simulates the amount, distribution, and seasonal timing of burned cropland and pasture over 2001-2009 (global totals: 0.434×10 6 and 2.02×10 6 km 2 yr −1 modeled, 0.454 × 10 6 and 2.04 × 10 6 km 2 yr −1 observed), but carbon emissions for cropland and pasture fire are overestimated (global totals: 0.295 and 0.706 PgC yr −1 modeled, 0.194 and 0.538 PgC yr −1 observed).The non-agricultural fire module underestimates global burned area (1.91×10 6 km 2 yr −1 modeled, 2.44 × 10 6 km 2 yr −1 observed) and carbon emissions (1.14 PgC yr −1 modeled, 1.84 PgC yr −1 observed).The spatial pattern of total burned area and carbon emissions is generally well reproduced across much of sub-Saharan Africa, Brazil, Central Asia, and Australia, whereas the boreal zone sees underestimates.FINAL represents an important step in the development of global fire models, and offers a strategy for fire models to consider human-driven fire regimes on cultivated lands.At the regional scale, simulations would benefit from refinements in the parameterizations and improved optimization datasets.We include an in-depth discussion of the lessons learned from using the Levenberg-Marquardt algorithm in an interactive optimization for a dynamic global vegetation model.

Introduction
Vegetation fire is an important force for the Earth system at local, regional, and global scales.It can shape ecosystems (Bond and Kelley, 2005;Staver et al., 2011), affect human health (Johnston et al., 2012;Marlier et al., 2012;Hahn et al., 2014), exacerbate or mitigate anthropogenic climate change (Ward et al., 2012;Ciais et al., 2013), and cause direct economic damage (Doerr and Santín, 2013;Bryant and Westerling, 2014).Fire occurrence can even affect the likelihood of more burning, through positive and negative feedbacks resulting from fire's impact on weather, climate, and vegetation (Laurance and Williamson, 2001;Balch et al., 2008;Zhang et al., 2008).Anthropogenic climate change and increases in atmospheric carbon dioxide concentrations have already increased -or can be expected to increase -the frequency and Published by Copernicus Publications on behalf of the European Geosciences Union.severity of burning in some parts of the world, while other regions could see decreased burning (Gillett et al., 2004;Westerling et al., 2006;Flannigan et al., 2009;Krause et al., 2014) A full accounting of the importance of vegetation fire to the Earth system at present as well as historically and into the future requires the use of dynamic global vegetation models (DGVMs).These simulate processes of vegetation establishment, growth, mortality, disturbance, and competition at large scales using varying levels of mechanism, which allows the regional-and global-level biogeochemical implications of ecosystem dynamics to be fully estimated.When DGVMs are coupled with models of the soil, atmosphere, and oceans, the resulting Earth system models (ESMs) even simulate how these major components of our planet interact with and feed back upon one another.To understand the complex nature of fire's role in the Earth system, realistic models of vegetation burning must consequently be designed and incorporated into DGVMs.
However, fire does not exist solely at the interface of climate and vegetation.Humans play an important role in regulating the fire regimes of many regions around the world (Flannigan et al., 2009;Bowman et al., 2011;Archibald et al., 2013).This can come about as a result of many processes, one of which is fire's use as a tool to manage agricultural lands.Croplands can be burned to facilitate planting or harvest; for example, sugarcane is typically burned before being harvested, and farmers in many parts of the world burn their crop wastes in the field after harvest (Yevich and Logan, 2003).Pastures and rangelands often see regular burning to reinvigorate the soil and control non-palatable weeds (Uhl and Buschbacher, 1985;Laris, 2002).
The way people burn croplands and pasture in a given region can differ from how the ecosystems there would burn in the absence of humans, in terms of both frequency and seasonal timing (Le Page et al., 2010;Magi et al., 2012;Rabin et al., 2015).This is significant for modeling efforts because it suggests a decoupling of agricultural fire from the mechanisms governing non-agricultural fire.For example, whereas the fire regime of southern Mali might naturally be dominated by large burns late in the dry season, humans have imposed a regime of small, scattered early burning to avoid such hard-to-control fires (Laris, 2002(Laris, , 2011)).
Unfortunately, previous development of global fire models has mostly glossed over the distinction between agricultural management burning and other burning.Anthropogenic effects on fire are most commonly modeled as dependent solely on population density, not land use (e.g., Venevsky et al., 2002;Arora and Boer, 2005;Pechony and Shindell, 2009;Thonicke et al., 2010;Li et al., 2012;Melton and Arora, 2016;Hantson et al., 2016;Rabin et al., 2017).Moreover, the effect of population density is only to increase or decrease the amount of fire relative to that which would occur naturally -not to affect the intra-annual timing of fire.There are a few exceptions.The LPJ-LMfire model (Pfeiffer et al., 2013) includes functions to simulate how pre-industrial societies could manage cropland and pasture using fire, but these depend on assumptions that may not apply as well to modern agricultural practices.A fire model developed for the Community Land Model (CLM) by Li et al. (2013) simulates cropland fire, with annual burned area based on socioeconomic data (population density and gross domestic product) and timing based on observations, but pasture is not simulated as a land cover/use type distinct from grassland.The HESFIRE model (Le Page et al., 2015) accounts for how the amount of human land use (cropland and urban areas) affects burning, but again pasture is not considered.Neither of these latter two models, moreover, take into account how human activity can affect the timing of fire.
To some extent, the neglect of pasture burning in particular -or its convolution with non-agricultural burning -can be attributed to a lack of data.Cropland and a number of other vegetation types can, like fire, be algorithmically mapped using medium-resolution satellite imagery.Overlaying maps of vegetation type and burned area allows the generation of observational datasets of fire activity on different land covers (e.g., Giglio et al., 2010).However, no such map of global pasture distribution exists -only maps at relatively coarse resolutions describing the fraction of each grid cell that is pasture (e.g., Ramankutty et al., 2008;Klein Goldewijk et al., 2010).When developers of global fire models have designed and parameterized models of non-agricultural burning, they have thus been limited in their choice of observational data with which to constrain their models.The options have been to either focus on regions with low fractions of cropland and/or pasture (thus potentially biasing their parameterization towards parts of the world inhospitable to agriculture) or to use a dataset "contaminated" with signals from cropland and/or pasture burning.Recently, however, Rabin et al. (2015) used a statistical method to estimate burned area associated with cropland, pasture, and non-agricultural lands at regional scales based on observations of total burned area and estimated land use/cover distributions.This presents an opportunity to create a fire model that not only explicitly simulates burning practices on cropland and pasture, but also to develop a model of non-agricultural burning based on a purer observational signal.
However, the choice of reference data is only the first step in model development.Model fitting, also referred to as optimization or parameterization, is also critical, and many different methods can be used.Empirical fire models have often been fitted against observations of weather, climate, vegetation state, and anthropogenic factors using regression-type methods (e.g., Archibald et al., 2009;Lehsten et al., 2010) or multidimensional search algorithms (Knorr et al., 2014).Because these methods treat fuel availability as an independent variable, however, they ignore how fire affects the fuel available for future burning.This fire-biomass feedback can be accounted for by running the fire model interactively with vegetation for parameterization purposes.This process is performed in combination with data from the literature when possible, but it is rather manual and based on trial and error.Ideally, model fitting would combine the best parts of these two approaches to algorithmically search parameter space for the "best" set of values based on how the model actually performs.Le Page et al. (2015) recently used the Metropolis Markov Chain Monte Carlo method to do just this in fitting the HESFIRE model.This stand-alone model accounts for fuel availability indirectly, with parameterizations based on precipitation and time since fire.Unfortunately, because of the need for high numbers of iterations, this method cannot be feasibly applied in fire models that are coupled with computationally expensive DGVMs.
Here we describe the development and performance of a DGVM-coupled fire model that uses the new disentangled estimates of burned area associated with cropland and pasture (Rabin et al., 2015) to enable true separation of fire patterns and processes between non-agricultural and agricultural land.A module for non-agricultural fire is fit against the purer, non-agricultural burning data -i.e., observational estimates excluding fire on cropland and pasture -using an algorithm that explores parameter space interactively with the fire and vegetation model.Cropland and pasture fire are explicitly simulated -for the first time, in the case of modern-day pasture fire -by a different module using derived climatologies.

Fire model
The Fire Including Natural & Agricultural Lands (FINAL) model comprises two different sub-models, simulating fire on agricultural and non-agricultural land separately.Here we describe the model's structure, beginning with the land and vegetation model within which it has been developed, then detailing the separate setups used for simulating nonagricultural and agricultural fire, and finally explaining the simulation of fire's effects on vegetation.

Land and vegetation model
The land model LM3, run by the National Oceanic & Atmospheric Administration Geophysical Fluid Dynamics Laboratory (NOAA-GFDL), is a state-of-the-art global dynamic vegetation and land surface model that can be run either offline or interactively with atmosphere and oceans in the GFDL Earth System Model (Shevliakova et al., 2009;Dunne et al., 2013).It simulates five different live plant biomass pools: leaves, heartwood, sapwood, labile carbon, and fine roots.The "stem" biomass pool is comprised of the heartwood, sapwood, and labile carbon pools.One of five different plant "species," representing biome types with different physiological properties, is assigned to each point based on bioclimatic envelopes and amount of biomass.Here, LM3 is run at a spatial resolution of 2 • latitude by 2.5 • longitude.
One of the most interesting features in LM3 is that it uses sub-grid cell units called tiles, which allow land in different land use types (and in different stages of recovery from land use) to have distinct simulated vegetation and soil.Grid cells can have one each of "natural," cropland, and pasture tiles, along with several "secondary" tiles representing land in different stages of recovery from wood harvesting or agricultural abandonment.Other, non-vegetated tiles represent glaciers and lakes.Tiles are not spatially arranged, instead existing effectively as a list within each grid cell.Wood harvest and land use transitions occur once per year.At the same time, secondary tiles are merged together if they have similar amounts of heartwood biomass; this prevents the computational burden from becoming unreasonable.
The tiled structure of LM3 could allow it to simulate the heterogeneity of vegetation that fire can create across a landscape, and cropland and pasture tiles could have fire occur in a completely different way than non-agricultural tiles.The original LM3 fire model did not burn cropland and pasture at all; elsewhere, fire happened once per year based on fuel loading, drought, and historical fire frequency (Shevliakova et al., 2009).The next two sections will describe the structure of the new fire models developed for non-agricultural (natural and secondary; Sect.2.2) and agricultural (cropland and pasture; Sect.2.3) tiles.

Burned area: non-agricultural land
The fire model for non-agricultural lands is based on that developed for the Community Land Model (CLM) by Li et al. (2012Li et al. ( , 2013)).Total burned area (BA) in the natural and secondary fire model is calculated as the product of the number of fires (N fire ) and burned area per fire (BA pf ): (1)

Number of fires
Lightning and humans both serve as sources of ignitions, some fraction of which actually become fires.Li et al. (2012) modeled their equation for the density of lightning ignitions after that elaborated by Prentice and Mackerras (1977).At each time step, the number of ignitions from lightning (I n , ignitions km −2 ) is a function of latitude ( , radians) and the density of lightning flashes (L, flashes km −2 ): The number of anthropogenic ignitions (I a , ignitions km −2 ) is a function of population density (people km −2 ):  3.) The second part of Eq. ( 3) is intended to represent the fact that each person can be expected to light fewer fires as population density increases (Venevsky et al., 2002).
To calculate the number of ignitions actually becoming fires (N fire ), the total number of ignitions (A T [I n +I a ], where A T is the area of the tile in km 2 ) is multiplied by five functions that vary from zero to one, representing the suppressive effects of relative humidity (f RH ), soil moisture (f θ ), aboveground biomass (f AGB ), temperature (f T ), and population density (f P D ): (4) Li et al. (2012) calculate the effect of relative humidity on number of fires as where RH (range 0-1) is the relative humidity in the tile.
Relative humidity ceases limiting fire (i.e., f RH = 1) below RH = 0.3, and it suppresses all fire above RH = 0.7.However, the artificial limitation of this formulation to the range [0, 1] would cause problems during our parameterization, which requires a continuously differentiable function.Instead we used the Gompertz function: This function also varies between zero and one, with the parameter β RH,1 controlling the location of the curve along the x axis and β RH,2 determining the steepness of the function as it decreases from one to zero.Li et al. (2012) formulate the effect of soil moisture on number of fires as where θ is relative soil moisture over the top 5 cm and θ e is a parameter determining the soil moisture level where approximately 95 % of fires are suppressed.This is a continuously differentiable function, but for consistency we used (like f RH ) a Gompertz function: In addition to flammability as determined by fuel moisture, Li et al. (2012) calculate the effect of above-ground biomass on number of fires as where AGB (kgC m −2 ) is the sum of aboveground biomass in the heartwood, sapwood, labile carbon, live leaf, and leaf litter pools.(80 % of the total biomass carbon in the heartwood and sapwood pools is assumed to be in the aboveground stem, with the remainder in coarse roots.)The parameters (kgC m −2 ) determine the levels of aboveground biomass below which fire is impossible (AGB lo ) and above which biomass is no longer limiting (AGB up ).However, as with f RH , the fact that this function is not continuously differentiable would create problems for parameterization, so we used a Gompertz function instead: The effect of temperature on number of fires is calculated as where T ( • C) is the temperature of the canopy.The T * parameters ( • C) serve the same purpose as the parameters in the original formulation of f AGB (Eq.9); that is, no fire can occur (f T = 0) at or below T lo and temperature does not limit fire (f T = 1) at or above T up .After Li et al. (2013), we set T lo to −10 • C and T up to 0 • C. Because we did not include this function in the optimization, we did not convert it to a Gompertz function as we did with f RH and f AGB .
The suppressive effect associated with increasing population density on all potential fires (as opposed to just anthropogenic ignitions, as accounted for in Eq. 3) is calculated as where P D is human population density (people km −2 ).f P D → 0.01 as P D → ∞, and f P D = 0.99 where P D = 0, after Li et al. (2012).β P D determines the shape of the function between these limits.Li et al. (2013) also included a suppressive effect of per capita gross domestic product (GDP) on number of fires.This was based on the idea that relatively wealthy parts of the world might have more valuable property to protect and a better capacity for suppression than less developed regions.However, for several reasons, we chose not to include this function.First, although globally gridded maps of GDP exist for the past 25 years or so (van Vuuren et al., 2007), no existing data sets describe the distribution of economic status before 1990.Second, the functions elaborated by Li et al. (2013) are somewhat ad hoc, not taking into account other variables that might be responsible for the observed relationships.Bistinas et al. (2014), for example, showed that an apparent relationship between GDP and burned area (Aldersley et al., 2011) can be better explained as an emergent property resulting from the effect of population density.That result does not deal with GDP per capita, of course, but it does indicate the care that must be taken to avoid confounding variables when modeling fire.We thus declined to include GDP effects on burning in our model.

Burned area per fire
Burned area per fire is calculated in the CLM fire model based on an approximation of individual fires having elliptical shapes, with the point of ignition being one focus and the fastest spread occurring along the major axis (Fig. S1; van Wagner, 1969).It is made up of three main components: duration, shape, and rate of spread.Up to a certain point, fires become more elongated with increasing wind speed.That is, higher winds increase the length-to-breadth ratio (LB; Fig. S1): where W is wind speed (m s −1 ) at 10 m above ground level.High winds also increase rate of downwind spread relative to the rate of upwind spread, which can also be thought of as increasing the head-to-back ratio (HB; Fig. S1).HB is related to LB as Forward rate of spread (ROS f , m s −1 ) -i.e., spread rate downwind from an ignition -is a function of wind speed, fuel moisture, and vegetation type.Vegetation type ("species" sensu LM3) determines the maximum possible rate of spread in a tile.We initially defined maximum rate of spread for each species (ROS max,sp ) based on values used by Li et al. (2012 and Corrigendum) for similar plant functional types (PFTs): 0.4 m s −1 for C 3 and C 4 grass, 0.3 m s −1 for tropical and evergreen trees, and 0.22 m s −1 for temperate deciduous trees.However, we included maximum rate of spread for tropical tree and C 3 and C 4 grass in the optimization (β ROStt and β ROSgr , respectively; Sect.2.6), so 0.4 and 0.3 m s −1 represent their starting values.Their final values can be found in Table 3.
Note that although Li et al. (2012 and Corrigendum) actually used 0.22 m s −1 for all forest types other than needleleaf, we increased the initial value of maximum rate of spread in tropical tree tiles closer to that given by Li et al. (2012 and Corrigendum) for shrub PFTs (0.34 m s −1 ).This was done because the rate of spread in tropical savannas is much higher than that in tropical closed forests (especially moist forests), but LM3 has no "shrub" or "savanna" species, with the result that much of the world's tropical savannas are classified as "tropical tree." The rate of spread realized by any given fire increases with wind speed towards the limit of ROS max,sp according to the function g(W ): where Here, LB max = 11 and HB max ≈ 482 are the limits of LB and HB as W → ∞ (Eqs.13 and 14).
Fires spread more slowly in wet conditions, so fuel moisture is considered in rate of spread.Li et al. (2012) multiplied rate of spread by f RH (Eq.5) as well as f RH (θ ), the latter being identical to f RH except with soil moisture (θ ) replacing relative humidity (RH).However, we substituted f RH (θ ) with f θ for simplicity and transparency.Thus, the complete equation for forward rate of spread in FINAL is as follows: The final component of burned area per fire is the length of time between ignition and extinction.After Li et al. (2012), we set fire duration (d, seconds) to 24 h (86 400 s).
Li et al. ( 2013) also include functions that reduce burned area per fire based on population density and GDP per capita.We did not include either of these.The issues with using GDP per capita are described in Sect.2.2.1.Population density might be considered a more trustworthy and meaningful statistic, but as with the GDP functions, the method used by Li et al. (2013) to describe the effect of population density on fire size was somewhat ad hoc and did not take into account possible confounding factors.Moreover, our model optimization (Sect.2.6) would have essentially seen the functions relating population density to number of fires and burned area per fire as one large, complicated function.For simplicity and parsimony, we consequently did not include an effect of population density on burned area per fire.
Several limits are imposed on BA pf .If the burned area calculated at a time step (i.e., BA pf × N fire ) is greater than the area of the tile that has not yet burned that day (A t,un ), BA pf is adjusted for consistency: Moreover, we add a limitation to fire size based on landscape fragmentation, based on the idea that fragmentation of the landscape into burnable and unburnable patches tends to prevent fires from reaching their maximum possible size (Archibald et al., 2009;Hantson et al., 2015).Maximum possible fire size as a function of tile size and fraction unburnable area in the grid cell is modeled after the function described by Pfeiffer et al. (2013): .
Here, A g refers to the area of land (including non-vegetated "land" such as glaciers or lakes) in the grid cell, and A g,burnable refers to the area of vegetated land in the grid cell other than cropland.BA pf,max is calculated at the end of each model day -after burning, tile splitting, and land-use transitions have occurred -and applied to the following day.Burned area is calculated at every fast time step (30 model minutes) and accumulates throughout each day.At the end of each model day, burning occurs (Sect.2.4).

Burned area: cropland and pasture
Burned area on cropland and pasture tiles is estimated in a simpler way than that on natural and secondary tiles.At the beginning of each month, some fraction of each cropland and pasture tile burns according to a mean monthly climatology of burned fraction of cropland and pasture.These gridded climatology maps are based on results from the "unpacking" analysis of Rabin et al. (2015), which provided monthly estimates of burned area associated with cropland, pasture, and non-agricultural ("other") land.For each of 134 regions around the world, using the GFED3s burned area data (Randerson et al., 2012), Rabin et al. (2015) estimated the F k,m parameters in the following equation: where the summation is over all N grid cells in the region, A k,i,m represents the area of each land use type (cropland c, pasture p, and non-agricultural land/"other" o) in grid cell i in month m, and BA m is the total burned area in the region in that month.This calculation was performed for each of the 108 months in 2001-2009.Each parameter F k,m thus represents the net influence of land use k on fire in the average grid cell in the region that month.In some instances, F k,m can be negative, which was interpreted to represent a suppressive influence of k on fire on other land use types.
Here we use the climatological mean results for F c and F p , constrained to non-negative values in order to focus on how much burning actually occurs on cropland and pasture, rather than including their suppressive influences: where k ∈ {c, p} and M ∈ [1, 12] is the month of the year corresponding to time step t.Note that, in Rabin et al. (2015), forcing F k,m ≥ 0 resulted in estimates of total burned area (i.e., burned area summed across all three land cover/use types) slightly greater than the value from GFED3s: 4.93 Mha yr −1 as opposed to 4.68 Mha yr −1 .Because the land cover distributions used in the unpacking (Rabin et al., 2015) differ slightly from those used in this study, burned fraction for each grid cell in the unpacked data was adjusted here so that the model output would match the burned area from the unpacking.

Fire effects
Carbon in the leaves, stems, and aboveground litter of a burned tile is combusted (i.e., transferred to the smoke pool; Sect.2.5) according to species-specific fractional combustion completeness (CC) values, based on those used by Li et al. (2012).The remaining non-combusted biomass in leaves, stems, and fine roots is subjected to species-and poolspecific fractional mortality (M; i.e., transferred to aboveor belowground litter), again based on values from Li et al. (2012).Combustion completeness and mortality values used here can be found in Table 1.Note that although the heartwood and sapwood pools are assumed to be 80 % aboveground ("stems") and 20 % belowground ("coarse roots"), CC stem and M stem are the same for both above-and belowground pools.This was necessary because LM3 assumes a constant 80-20 % split.However, fire-killed heartwood and sapwood is transferred to aboveground or belowground litter proportionally.
If less than 1 km 2 of a tile burns, the tile's biomass is reduced according to CC × BF and (1−CC) ×M× BF, where BF is the burned fraction of the tile.This is the method that has been used by every other global fire model previously developed.However, it does not reflect the reality that an actual fire results in a mosaic where only part of the landscape has been burned.To better represent this process, when ≥ 1 km 2 burns in a given day, FINAL splits the tile into two new tiles -one burned and one unburned.Biomass on the burned tile is reduced by CC and (1−CC) ×M, while the unburned tile is not affected.This "fire tile splitting" occurs on all land cover types except cropland.The 1 km 2 threshold was set to reduce computational demand and avoid calculation errors associated with small tiles.

Other changes
The implementation of daily fire and associated tile splitting necessitated many adjustments to parts of the LM3 code base not dealing with fire directly.Previously, tiles would only be created and/or merged once per year, and secondary vegetation was the only land type allowed to have multiple tiles within a single grid cell.The code for land transitions needed to be reworked to allow daily splitting and merging.We also changed the code to allow all vegetation types, instead of just secondary land, to have multiple tiles.The criteria for merging tiles were also altered to be based on aboveground biomass available for fire (AGB in Eq. 9) instead of heartwood.Moreover, we changed the binning structure by which tiles are determined to have similar-enough biomasses to be merged.Previously, bin edges were located at 0.5, 1, 2, 3, 4, 5, 6, 8, 10, and 1000 kgC m −2 .To better sample ranges of biomass where fuel is limiting, we replaced the first two bin edges with 0.1, 0.3, 0.5, 0.7, 0.9, and 1.1 kgC m −2 .Finally, various aspects of carbon accounting throughout the model needed to be adjusted for daily tile splitting and merging.
More frequent fire also required other changes.The original LM3 fire module burned once annually at the end of each year, with the burned carbon being emitted gradually over the course of the next year to avoid sudden unrealistic pulses of emissions.With the new fire model operating daily, burned carbon from one day is now emitted over the course of the next day.Previously, grazing of pasture happened once per year, but in order to more reasonably simulate emissions from pasture fire we made grazing occur daily.We also boosted the fraction of live leaf biomass removed by grazers from ∼ 0.07 to 4 % day −1 for the main runs (FINAL.0 and FINAL.1;Table 2).This resulted in more realistic estimates of aboveground biomass in pasture, and of annual global consumption of biomass by grazers.
Finally, the original LM3 model did not explicitly simulate aboveground dead biomass, which is an important component of the fuel bed in some ecosystems, affecting fire spread and/or emissions.We thus used the version of LM3 with the Carbon, Organisms, Rhizosphere, and Protection in the Soil Environment model (CORPSE; Sulman et al., 2014), which in addition to simulating the dynamics of soil organic matter also simulates leaf litter and coarse wood litter pools.The

Parameter optimization
Simply copying exact parameters from the model described by Li et al. (2012Li et al. ( , 2013) ) was not possible for a number of reasons.First, here we separately model cropland, pasture, and non-agricultural burning.Li et al. (2013), in comparison, included special modules for cropland, deforestation, and peat fire -pasture burning being convolved with all other fire.Now that we have extracted the influence of pasture (a significant source of fire activity that often differs from what might be expected under a totally "natural" fire regime) from non-agricultural burning, we expect to find different relationships between fire and its driving variables.Second, CLM is of course a different model than LM3, with its own idiosyncrasies and biases distinct from those of LM3.
Although Li et al. (2012Li et al. ( , 2013) ) strove to parameterize their equations based on independent data as much as possible, some functions were entangled with how their model itself worked.Third, as described in Sect.2.5, we added some processes and removed others.Fourth, Li et al. (2012Li et al. ( , 2013) ) tested their model against version 3 of the Global Fire Emissions Database (GFED3) burned area dataset (Giglio et al., 2010), whereas we used the GFED3s dataset (Randerson et al., 2012), which includes an additional estimate of burning from small fires and thus has significantly more burned area than GFED3.Finally, Li et al. (2012Li et al. ( , 2013) ) used different climatic forcing data than we did.All these differences meant that we needed to reparameterize at least some parts of the non-agricultural fire model.Here we begin by briefly walking through the algorithm used to carry out the optimization, and then describe the parameters that we chose to optimize.

The Levenberg-Marquardt algorithm
We used the Levenberg-Marquardt method as the basis of our optimization routine.This algorithm uses the first derivatives of a performance metric with respect to each parameter to iteratively move through parameter space in search of a local minimum of the sum of squared errors.It starts with an initial guess, then evaluates the sum of squared errors S in non-agricultural burned area between the unpacked data and the estimates generated by the model: (Here, the summation is performed across all M months in the parameterization run period and all N sample grid cells selected for the optimization.)The algorithm then generates a new parameter set guess and the model is rerun.If the new guess decreases the sum of squared errors, it is "accepted," with a new guess then being generated based on it.
If not, it is "rejected," and a new guess is generated based on the original guess.Guesses are adjusted by interpolating between steps that would be generated by either the gradient descent method or the Gauss-Newton algorithm, leaning more towards the former when far from a minimum and the latter when near a minimum.More detailed information on the Levenberg-Marquardt algorithm, including its derivation, can be found in Levenberg (1944), Marquardt (1963), and Transtrum and Sethna (2012).
Briefly, we ran the model for 1991-2009 in a sample of 241 grid cells.A Python script evaluated the model performance and suggested a new parameter set, which was fed back into the model.The Python script then checked the performance of the new parameter set, accepted that set if its performance was improved relative to the previous set, and generated a new guess.This process continued until the routine encountered at least five rejected parameter set guesses.We did not optimize over all grid cells because of computational limitations; even with all 241 grid cells being run in parallel, each iteration of the optimization took around two hours.More details on our implementation of the algorithm, including how the grid cells in the sample were selected, can be found in Appendix A.
Four optimizations were performed, with slightly varying initial conditions (Table 3) in order to enhance the robustness of the results.Optimization 1 was performed using the parameter values from the literature or from Gompertz curve fitting; Optimizations 2-4 used parameter values sampled from a ±25 % uniform distribution around the Optimization 1 values.Note that we began the optimization runs in 1991 even though only the 2001-2009 data would be used for comparison to observations; the idea was to allow for the vegetation and fire regime in at least some of the grid cells (especially in regions where frequent fire is the norm) to equilibrate given the fire frequency of each new iteration of the model.

Parameters chosen
From the equation for anthropogenic ignitions (I a , Eq. 3), we optimized β I a , which can be thought of as controlling a sort of "baseline" value for how many ignitions each person can be expected to provide at each time step.Technically, we optimized β I a ,m , which is describes the baseline number of ignitions per person per month instead of per time step (of which there are 48 per day): All other things being equal, higher values of β I a ,m result in more fires.
We also optimized β P D from the function describing human suppression of all non-agricultural fires as a function of population density (f P D , Eq. 12).All other things being equal, a higher value of this parameter would result in a faster approach of the fraction suppressed towards its upper limit.
Because the LM3 definition of a "species" to describe vegetation type is so broad, we thought it would be especially important to pay attention to several biome-specific maximum rate of spread parameters in FINAL.The "tropical tree" type in LM3 encompasses a wide range of real-world biomes, from tropical rainforests to semiarid shrublands.The rates of spread for fire in these systems are quite different, and so we included maximum rate of spread in tropical tree regions (β ROStt ) in the optimization.We also included the rate of spread in C 3 and C 4 grasslands (β ROSgr ), because preliminary testing showed strong overestimates in regions dominated by the C 4 grass species in particular.
Finally, we optimized parameters from f RH (β RH,1 and β RH,2 , Eq. 6), f θ (β θ,1 and β θ,2 , Eq. 8), and f AGB (β AGB,1 and β AGB,2 , Eq. 10).We generated initial guesses for these parameters by fitting Gompertz functions, with the upper asymptote set at 1, to the corresponding functions from Li et al. (2012).Fitting was performed using the MATLAB Curve Fitting Toolbox (MATLAB and Curve Fitting Toolbox Release 2014b, The MathWorks, Inc., Natick, Massachusetts, United States.) Note that we did not optimize all possible parameters in the model.For example, we did not include the parameters affecting the upper and lower asymptotes of f P D (Eq.12) in the interest of limiting the degrees of freedom with regard to the combined population density functions.Given that we were already optimizing two parameters governing the effect of population density on number of fires (β I a ,m and β P D ), we decided to exclude the other parameters in Eq. ( 12).The rest of the parameters that we did not optimize were excluded in the interest of somewhat limiting the scale of the optimization procedure.This is especially true with regard to the parameters in Eq. ( 13) (governing the effect of wind speed on fire length : breadth ratio) and Eq. ( 20) (governing the effect of decreasing burnable area on maximum fire size).The parameters in these equations are generally based on phe-nomena external to global vegetation modeling -Eq.( 13) is derived from empirical equations used by the Canadian Forest Service (Arora and Boer, 2005), and Eq. ( 20) is derived from an experiment performed by Pfeiffer et al. (2013) independent of any fire or vegetation model.Because tropical savanna grid cells, with the highest initial sums of squared errors, were expected to exert the most influence on the optimization procedure (Fig. A3), we focused on optimizing parameters regarding variables that are known to be influential there.Other parameters -such as temperature, or the rate of spread in boreal forests -might not have been wellconstrained in this procedure because of their low importance in cells with high initial error.
3 Experimental setup and analysis

Experimental runs
Spinup of the land to pre-industrial conditions began with a "bare ground" scenario and ran for 300 years, during which climate forcings (Sect.3.2) from 1948 to 1977 were repeatedly cycled through.During spinup, atmospheric CO 2 concentration was held constant at 286 ppm and land use was turned off.Next, we simulated years 1861-1947, using repeated 1948-1977 climate forcings but historical land use and atmospheric CO 2 concentration (Sect.3.2).Finally, the model was run from 1948 to 1991 with historical climate forcings, land use, and atmospheric CO 2 .This run -referred to as LM3_ORIG (Table 2) -provided initial conditions for other model runs, including the optimization.Note that the daily grazing intensity (Sect.2.5) was set at its default value of ∼ 0.07 % for LM3_ORIG.
The new model (Sects.2.2-2.5), with new parameters as described in Sect.4.1 and  tings as initially guessed in the parameterization (FINAL.0)was also performed, which when compared to FINAL.1 would allow us to explore where the optimization improved or worsened model performance.For both FINAL.0 and FINAL.1, daily grazing intensity (Sect.2.5) was set at 4 %.

Input data
The LM3 land and vegetation model is run "offline" in this study, meaning that it is forced by a set of meteorological and radiation-related variables without any interaction between the land and atmosphere.The variables used here to force LM3 -daily precipitation, surface air pressure, specific humidity, wind vectors, and downward longwave and shortwave radiation -are taken from the observation-based dataset developed by Sheffield et al. (2006).All variables are interpolated to the spatial and temporal resolution of the LM3 fast time step, here set to 30 model minutes.Carbon dioxide (CO 2 ) concentrations are taken from Meinshausen et al. (2011).Historical data on land use transitions and wood har-vesting come from the harmonized dataset created by Hurtt et al. (2011) for use in Earth system models.The mean distributions of cropland, pasture, and non-agricultural land in this study over 2001-2009 are presented in Fig. 1.As discussed above (Sect.2.3), cropland and pasture burning is forced using climatologies from (Rabin et al., 2015).For the non-agricultural fire model, we used a gridded monthly climatology of lightning flash rate (flashes km −2 ) based on data from the Lightning Imaging Sensor (LIS) and Optical Transient Detector (OTD) remote instruments.Specifically, we used the LIS/OTD low-resolution monthly time series (LRMTS) described by Cecil et al. (2014).This dataset is provided at a 2.5 • × 2.5 • resolution, which we interpolated to match the LM3 resolution of 2 • latitude by 2.5 • longitude.The version of LRMTS that we used, v2.3, included maps of flash rate for each month in the period 1996-2014.We found the average of each month (January, February, etc.) and used these to build our climatology.
Non-agricultural burning in FINAL also requires input data on population density.We used the historical population density estimates from HYDE 3.1 (Klein Goldewijk et al., 2010), coarsened from their original 5 min resolution to the LM3 resolution (2 • latitude by 2.5 • longitude).We interpolated population density linearly between each time point in the HYDE dataset.

Evaluation
The new model's performance in terms of recreating observed patterns of burned area and fire carbon emissions is evaluated here by comparison against GFED3s and the unpacked fire data.In addition to global totals of mean annual fire activity, we assess the spatial distribution of fire using maps of mean annual burned fraction and emissions.Unfortunately, due to the short satellite record of fire occurrence, the model must be evaluated against the same time period used for calibration.The model can thus be expected to perform less well outside 2001-2009.
The accuracy of seasonal fire trends is tested by comparing the difference between the intra-annual timing of burned area simulated by the model with the timing as estimated by the unpacking analysis.This is quantified using mean phase difference, as described by Kelley et al. (2013).Each grid cell's annual pattern of fire can be described as a vector in the complex plane: where x m,i is the mean burned area in month m for grid cell i, and θ m is an arbitrary angle unique to month m and calculated for all grid cells as The mean vector L i for each grid cell has end points that can be described in Cartesian coordinates as the origin and  x m,i sin (θ m ) .
The phase (P i ), defined where fire occurrence is not distributed evenly across all months, describes the mean timing of the fire season: The phase in terms of the day of the year can be calculated as 2π × 365.Mean phase difference (MPD), which is used here to describe the difference in timing of the fire season between model results and observations, is calculated as where modeled and observed phases are designated with the subscripts mod and obs, respectively.MPD varies from zero to one, with MPD = 0 if all modeled phases correspond exactly to observed phases and MPD = 1 if all modeled phases differ from observed phases by the maximum possible amount (6 months).

Optimization
Of the four optimization runs performed, only three completed successfully (Appendix A).In Optimization 1, the algorithm repeatedly increased β RH,2 , resulting eventually in model crashes.Optimizations 3 and 4 resulted in similar final functional forms; we chose to discard Optimization 4 since its final sum of squared errors (SSE; 3.667 × 10 9 ) was higher than that of Optimization 3 (3.657× 10 9 ).We were thus left with Optimizations 2 and 3; we used the final parameter sets from both of these for global model runs.Optimization 2 initially seemed like it might be the better candidate, since the SSE of its final parameter set (3.240 × 10 9 ) was lower than that of Optimization 3.However, although Optimization 2's final guess performed better in the selected grid cells during optimization, it actually performed worse than Optimization 3's best guess -and indeed, worse than Optimization 2's initial guess!-when run for the entire globe.This suggests that using SSE as the sole criterion for model selection is not sufficient.This issue, and the specifics of the Optimization 2 results, will be discussed further in Sects.5.3 and 5.4.For the remainder of this paper, except where specified, results will refer to those from Optimization 3 and the global model runs using its final parameter set.In this section, we discuss only the raw results from  S2) and Optimization 2 (Fig. S3).Optimization 3 resulted in final parameter values and functional shapes broadly similar to the initial guesses.Figure 2 shows the progression of the parameter guesses, along with the sum of squared errors associated with each parameter set guess, through Optimization 3.After an initial drop in SSE over the first six guesses, subsequent guesses did not result in much improvement, with SSE not differing by more than 0.001 % between accepted guesses after the 19th iteration (Fig. 2a-b).The optimization was stopped after the 42nd iteration, at which point seven consecutive guesses had been rejected.The functions resulting from the new parameter set are visualized, in comparison with how they were in the Li et al. (2012Li et al. ( , 2013) ) model as well as in the initial optimization guess, in Fig. 3.
In Optimization 3, the density of anthropogenic ignitions I a decreased at all positive levels of population density (Fig. 3a).Moreover, the parameter β P D -which controls anthropogenic suppression of burning f P D -increased (Fig. 2e), meaning that a larger fraction of ignitions (both lightning and anthropogenic) are suppressed wherever population density is greater than zero, though most noticeably between densities of ∼ 10-100 people km −2 (Fig. 3b).The net effect is to reduce unsuppressed anthropogenic ignitions (i.e., I a × f P D ) relative to the initial guess: the peak dropped from 3.6×10 −5 to 1.8×10 −5 ignitions day −1 , with the location of the peak shifting from 18.6 to 9.1 people km −2 (Fig. 3e).
The f RH and f θ functions in Optimization 3 do not differ much from the initial guess to the final accepted guess (Fig. 3c-d).(It should be noted, however, that the initial guesses for parameters in f RH resulted in a less suppressive function than in Li et al. (2012Li et al. ( , 2013)), while the initial f θ was more suppressive.)Altogether -i.e., taking into account moisture effects on both ignition success probability and rate of spread -burned area in FINAL is proportional to (f RH × f θ ) 3 .The bulk of this net impact on flammability caused by the changes to f RH and f θ is concentrated in the range of 0-20 % soil moisture and 0-50 % relative humidity, with grid cells in this zone seeing a reduction in flammability (i.e., fraction of unsuppressed ignitions becoming fires) of around 0.1 between the initial and final guesses.The impact of the changes to the moisture functions is most clearly seen in the Sahara Desert (Fig. S4d).
β AGB,1 increased and β AGB,2 decreased (Fig. 2).These changes resulted in a rightward shift of the function and a decrease in the slope from low to high biomasses (Fig. 3).Biomass is thus more limiting in Optimization 3's final parameter set than in its initial one.Whereas the original parameter set gave f AGB = 0.99 at AGB = 1.67, the final function does not reach that value until AGB = 2.52.
Optimization 3 saw maximum rate of spread decrease more than 35 % for grassland (Fig. 2l) between the initial and final guesses, a result which likely has to do with the model overestimating fire in these low-biomass systems.This parameter decreased sharply for most of the optimization, but as f AGB appropriately began to take on more of the responsibility for regulating fire there, grassland maximum rate of spread began to increase back towards its initial guess.Maximum spread rate nearly doubled for the "tropical tree" vegetation type (Fig. 2k), due to a tendency towards underestimation of burned area in that biome.
Comparing the results of FINAL.0 with FINAL.1, we can see that much of the improvement came in regions where the initial parameter set severely overestimated burned area (Fig. 4).A map of sum of squared errors (Fig. 4d) can be used to visualize performance improvement as would be "seen" by the optimization algorithm for included grid cells.Arid regions tended to see the most improvement, as evident in the Sahara, parts of the western United States, the dry savannas and shrublands of Africa and Australia, and the west and central Asian steppes.Moister savannas, as well as the Caatinga, were most negatively impacted by Optimization 3; the boreal zone and Southeast Asia also suffered but to a lesser degree.A map showing model-output SSE change of only the optimized grid cells (Fig. S2c) suggests that, contrary to our expectations, tropical savanna regions did not dominate the optimization.Instead, relatively lower-burning dry subtropical savannas and temperate steppes saw the largest improvements, with tropical savannas often seeing worsened performance.

Burned area
Figure 5 compares, over 2001-2009, maps of mean annual burned fraction (i.e., fraction of land area) from run FINAL.1 with those from GFED3s (Randerson et al., 2012) and the unpacking analysis.Figure 6a shows the difference in mean annual burned fraction between the model and the unpacked observations, against which the non-agricultural model was parameterized.Considering all land cover types combined, the new fire model recreated the general pattern of annual fire activity well compared with both GFED3s (Randerson et al., 2012) and the unpacked data (Figs.5a, b, f; 6a).The largest modeled overestimates relative to the unpacked data occurred in the grasslands and shrublands of western South America, the western Caatinga of northeast Brazil, and at various points throughout the African savannas (Fig. 6a).Most of the severe model underestimation relative to the unpacked data occurred in the African tropical savannas, as well as (to a lesser extent) the tropical savannas of northern Australia (Fig. 6a).The modeled burned fractions of cropland and pasture match the unpacked numbers almost exactly (Figs.6c,d), which is not surprising considering that the unpacked data were used to force the model on cropland and pasture tiles.There are some notable discrepancies, however.Specifically, there is too much cropland fire in one European grid cell and too little in several grid cells in northern Australia (Fig. 6c).Pasture fire did not experience such severe error in burned fraction anywhere (Fig. 6d).
The strong correspondence of modeled cropland and pasture fire with the unpacked observations (as expected since the latter were directly used to drive the former) suggests that the majority of the error seen in total burning must be associated with fire on non-agricultural lands.Indeed, although the non-agricultural fire model generally captured the worldwide distribution of fire -with tropical savannas, grasslands, and shrublands generally dominating burned area -the fit is by no means perfect (Fig. 6b).There are a number of regions where the model simulates little to no non-agricultural burning but the unpacked data show significant amounts of fire (Figs.5e,i) .This phenomenon is especially noticeable in the eastern African savannas, the shrublands of western Australia, and throughout the tropical and temperate grasslands, savannas, and shrublands of South America.
Worldwide, the non-agricultural fire model underestimated burned area, with 1.91 × 10 6 km 2 yr −1 simulated as having burned -an underestimate of 22 % relative to the unpacked estimate (Table 4).Unsurprisingly given the spatial results presented above, global averages for cropland and pasture were much better -0.434 × 10 6 km 2 yr −1 (4 % underestimate) and 2.02 × 10 6 km 2 yr −1 (1 % underestimate), respectively.Mean annual global burned area across all land covers over 2001-2009 was modeled as 4.36×10 6 km 2 yr −1 , an underestimate of 6.7 % relative to GFED3s and an underestimate of 12 % relative to the unpacked total.The time series of annual burned area over 2001-2009 for each land cover from the model (i.e., FINAL.1) are compared with the GFED3s and unpacked estimates in Fig. 7a.
The non-agricultural fire model performed well in terms of simulating the within-year timing of burned area (Figs.S5e,  i).This was reflected in the results for combined burning across all land cover types, which corresponded well with both GFED3s and unpacked burned area (Figs.S5a-b, f); the phase of model-estimated fire was 32 days later than observed for all fire combined as compared with total unpacked fire (mean phase difference MPD = 0.18), and 49 days later than observed for non-agricultural fire specifically (MPD = 0.27).

Carbon emissions
Just as the model tended to underestimate total global burned area, it also underestimated carbon emissions from fire (Table 4).The 2.14 PgC yr −1 simulated by the model represents an underestimate of 14 % relative to GFED3s and of 17 % relative to the unpacking data.This is again principally due to non-agricultural fire, for which the model simulated 1.14 PgC yr −1 as opposed to the unpacked estimate of 1.84 PgC yr −1 -an underestimate of 38 %.Agricultural fire emissions were actually overestimated, with 0.295 PgC yr −1 for cropland and 0.706 PgC yr −1 for pasture -overestimates of 52 and 31 % compared to the unpacked values of 0.194 and 0.538 PgC yr −1 , respectively.
The spatial distribution of errors in total fire carbon emissions (Fig. 6e) generally reflects the distribution of errors in simulated burned area (Fig. 6a).As with burned area, there are sizable regions where the model simulates little to no  (2015).Solid blue lines: observations of total emissions from GFED3s (Randerson et al., 2012).Other solid lines: model-estimated total and by-land cover fire emissions.
non-agricultural fire carbon emissions but the unpacked data show otherwise (Fig. 8e, i).Cropland fire emissions, as with burned area, are underestimated in northern Australia; there are also two regions in central Africa where cropland fire emissions are overestimated despite essentially correct annual burned fraction (Fig. 8c, g).The areas of slightly underestimated pasture burned fraction are not apparent on the map of pasture fire emissions error; large overestimates of emissions from pastures in the tropical savanna biome are instead the most apparent aberrations (Fig. 8d, h).

Model performance in context: burned area
In terms of spatial distribution, the model tends to overcluster non-agricultural burned area relative to the unpacked estimate.That is, it tends (especially in savanna regions) to simulate a highly spatially heterogeneous distribution of nonagricultural burned area, with some areas burning very little and others burning far too much (Fig. 5).It is important to consider, however, that although the unpacking method generates accurate estimates of total burned area at the level of each analysis region, the burning tends to be too evenly distributed within each region (Rabin et al., 2015).This results in an overly smooth map, as can be seen by comparing maps (a) and (b) in Fig. 5. Non-agricultural burning in the real world might thus exhibit more spatial clustering than is apparent in Fig. 5e.The burned area smoothing resulting from the use of relatively large unpacking regions in the boreal zone (and especially in Russia; Fig. 1 in Rabin et al., 2015) may have contributed to the model's poor performance there.
To get a sense of the spatial clustering of real-world nonagricultural fire, we have constructed a map of mean annual "GFED3s non-agricultural" burned fraction by subtracting unpacked cropland and pasture burned fraction from mean annual GFED3s total burned fraction.(The exact numbers from this map are not very meaningful, since it is possible to have values less than zero in grid cells where unpacking estimated more cropland and pasture burning than all burning observed by GFED3s; the purpose of this exercise is only to examine spatial heterogeneity.)A map of the coefficient of variation in 6 × 6 grid cell (12 • latitude × 15 • longitude) kernels across this map is compared with similar maps for mean annual modeled and unpacked non-agricultural fire in Fig. S6.As expected, the coefficient of variation is much higher in the GFED3s data than the unpacked data, indicating stronger spatial clustering of non-agricultural fire in the real world.The fact that the model simulates more heterogeneity than the unpacked estimate therefore indicates that the model is capturing heterogeneity in fire drivers that are important to actual fire patterns.This is not to say, of course, that the heterogeneous patterns simulated by the model exactly match the observations -in some places they do not, as is apparent in Fig. 5.
Although savanna regions may have shown the largest absolute difference in modeled vs. unpacked fire activity, smaller differences can be just as important in other areas.For example, the GFED3s and unpacked data show a mean annual burned fraction of 1-5 % for the boreal forests of central Alaska and northwestern Canada (Figs. 5a-b, e), which would correspond to a mean fire return interval of 20-100 years.While this is a low rate of burning relative to tropical savannas, for example, it still represents an important process for the structure and function of that ecosystem.The non-agricultural fire model captures almost no boreal forest fire whatsoever (Fig. 5i), which should hamper the ability of LM3 to accurately simulate vegetation there.One possible contribution to this deficit is the importance of multi-day fires in the boreal region.We followed Li et al. (2012) in assuming that all fires last 24 h, but this assumption is not well-supported by the literature.Korovin (1996) found that almost 60 % of forest fires in Russia over 1947-1992 lasted longer than one day, and that fires lasting longer than 10 days accounted for nearly 70 % of the burned forest area.Stocks et al. (2003) found a similar importance of very large (and thus presumably long-lasting) fires in Canada, with individual burns of more than 20 000 ha comprising over 65 % of mean annual burned area over 1959-1997.Ideally, FINAL would replicate this pattern by explicitly modeling the duration of individual fires based on evolving weather conditions.Several global fire models have introduced such a component, but with mixed results.The LPJ-LMfire model developed by Pfeiffer et al. (2013), which allows fires to burn for about four hours per day until they experience significant precipitation, actually tends to overestimate boreal forest fire.The HESFIRE model (Le Page et al., 2015) also allows fires to burn indefinitely, calculating twice per day an extinction probability based on fuel load, attempted suppression intensity, landscape fragmentation, and weather conditions.However, like FINAL, HESFIRE simulates too little fire in the boreal region (Le Page et al., 2015).A new version of FINAL, FINAL.2, does include multi-day fire, and is successfully able to reproduce the distribution of fire frequency binned by duration in boreal Canada.However, even with that and other changes impacting fire behavior in the boreal zone, FINAL.2 still does underestimate burned area there (Ward et al., 2018).

Model performance in context: emissions
The tendency of FINAL.1 to underestimate total global 2001-2009 burned area is reflected in an underestimate of the associated carbon emissions -by 7 and 14 %, respectively, relative to GFED3s (Table 4).GFED3s and the unpacking data show respective average emissions densities of 0.53 and 0.52 kgC m −2 of burning for all fire combined, whereas FINAL.1 gives 0.49 kgC m −2 (based on Table 4).The largest discrepancy in fire carbon emissions density between the modeled and unpacked estimates is on cropland, where FINAL.1 simulates 0.68 kgC m −2 but the unpacking analysis gives only 0.43 kgC m −2 (59 % overestimate; Table 4).Emissions density on pasture are also overestimated, by 33 %; non-agricultural emissions density is underestimated by 21 %.
Given how extensive pasture burning is at a global scale, it is especially important to understand why the C emissions from pasture fire were so significantly overestimated -especially in the tropics (Fig. 8d, h).Emissions from pasture fires, as with all fires, are the product of three quantities: burned area, aboveground biomass, and combustion completeness.
The model simulates burned pasture area quite accurately (Table 4), and it is unclear how incorrect combustion completeness would affect emissions over long time scales.We thus conducted a detailed examination of grazing intensity and pasture biomass simulated by LM3, which can be found in Appendix B. Briefly, it appears that excess dead woody material is to blame for overestimates of pasture (and probably cropland) fire emissions density.At least in the tropics, this is partially due to the fact that much slash wood remaining during forest clearance is simulated as being left on the ground, when in reality it is mostly burned away shortly after cutting.The fact that LM3 uses a single global value for grazing intensity may also contribute to mis-estimates of pasture burning emissions, but this is likely outweighed by the effect of dead wood.
Note that the records of fire emissions in the GFED product are not purely observation-based.GFED emissions estimates are generated by forcing a version of the Carnegie-Ames-Stanford-Approach (CASA) model with GFED burned area, using vegetation type and soil moisture to determine combustion completeness (van der Werf et al., 2006Werf et al., , 2010)).Biases may exist in that model that result in incorrect estimates of aboveground biomass and/or combustion completeness.Apparent discrepancies between GFED3s and FINAL-simulated fire emissions thus may not represent true errors by FINAL relative to reality.

What do optimization results suggest?
Optimization resulted in fewer anthropogenic ignitions and stronger anthropogenic suppression for any given value of population density (Fig. 3a-b, e).This suggests that, by grouping together non-agricultural fires with pasture fires, previous modeling efforts may have overestimated the contribution of humans to burning on non-agricultural land.That is, by extracting a "pure" non-agricultural fire signal, our study shows that pasture burning practices may have been responsible for much of what was once characterized as general anthropogenic fire, and that humans enhance fire on non-agricultural lands less than once believed.In terms of the general shape of net anthropogenic influence on nonagricultural fires -including the location and width of the peak -our results do not differ substantially from the function described by Pechony and Shindell (2009) or that used by Li et al. (2012;Fig. 3e).Knorr et al. (2014), in comparison, used the Levenberg-Marquardt algorithm to fit a simple empirical fire model in a non-interactive fashion and found that the peak was actually located closer to a population density of 0.1 people km −2 than to the value of ∼ 10 people km −2 that we found here.
When considering the results of this optimization, it is important to keep in mind that even if the Li et al. (2012Li et al. ( , 2013) ) model had been used in LM3 without modification, performance would have differed from the original CLM version.Structural differences between CLM and LM3 result in dif-ferent vegetation dynamics and micrometeorology relevant for fire.We also used a different source for climate forcing data and calibrated our model based on different burned area data.These and other differences create uncertainty about exactly why any given function's parameters shifted as they did during our optimization.The Fire Model Intercomparison Project (FireMIP; Rabin et al., 2017) could be informative in this regard.
As mentioned in Sect.4.1, Optimization 2 performed better than Optimization 3 in the 241 grid cells chosen for optimization but not when considering all grid cells.The greatest differences between these two final parameter sets have to do with anthropogenic ignitions.Whereas Optimization 3 resulted in decreased ignitions per person and increased suppression (Fig. 3), Optimization 2 decreased suppression and greatly increased ignitions per person (Fig. S7).This extreme human burning parameterization -far in excess of empirical estimates and functions in other fire models -explains the worsened performance of Optimization 2's final parameter set in Europe, South and Southeast Asia, and the eastern United States (Fig. S8).Optimization 2 also resulted in a "backwards" shape for f RH , where lower relative humidity results in less fire.Additionally, f RH is never much greater than 0.2, and f θ is never much greater than 0.5.The net result of Optimization 2 is that it performs worse than Optimization 3 across much of the United States, Central America, Europe, South and Southeast Asia, and Australia (Fig. 9).On the other hand, Optimization 2 outperforms Optimization 3 across most of the boreal zone and in the Cerrado; the two perform similarly in the southern African savannas but Optimization 2 performs slightly better in the north (Fig. 9).

Levenberg-Marquardt optimization: lessons learned
One of the limitations of the Levenberg-Marquardt algorithm is that it can only "move downhill."At every iteration, it searches for new parameters in the direction of lower sum of squared errors from the current point in parameter space, even though the set of parameters with the lowest possible sum of squared errors may be in a totally different direction.
As an analogy, imagine a person given the task of finding the lowest point in a city.Using a "downhill-only" algorithm, this person would literally walk downhill from their starting point and stop when they reach a point -the local minimumwhere continued travel in any direction would be uphill.The person might more thoroughly search the city for its lowest point by occasionally turning uphill and/or randomly taking a bus to a totally different part of the city -analogous to the behavior of the Metropolis-Hastings or simulated annealing algorithms.Levenberg-Marquardt being a downhill-only algorithm is not a fatal flaw, especially when the initial parameter set guess is well-informed based on the literature.It may well represent an improvement in methodology over the manual trial-and-error approach.But it is important to re- member that Levenberg-Marquardt should not be expected to produce the universally best possible parameter set.Another, potentially more serious limitation of the Levenberg-Marquardt algorithm is its use of the sum of squared errors (SSE) as a metric to gauge model performance.While the setup used here does account for accuracy of burned area simulations in both space and time, SSE tends to result in a bias towards improving performance in grid cells where the model simulates burned areas much higher or much lower than observations.This tendency to reduce absolute error would be fine if the goal of optimization were to produce a model that accurately simulates burned area for its own sake, but relative error can be more reflective of how well the model simulates the state of the vegetation.For example, assume two hypothetical 1000 km 2 grid cells: one dominated by tropical grassland, where observations show 100 % annual burning but the model simulates 25 %; and one dominated by boreal forest, where observations show 1 % annual burning but the model simulates 0.25 %.In both cases, the model is producing 75 % less fire than what actually occurs -a difference that could be extremely important to the simulated structure and function of both ecosystems.However, because the absolute error in the grassland grid cell (−750 km 2 yr −1 ) is so much greater than that in the boreal forest grid cell (−7.5 km 2 yr −1 ), the former will, all other things being equal, have a much greater influence on the direction and magnitude of the step towards the next parameter set guess.Our use of an equirectangular grid -with cells of constant size in terms of latitude and longitude but not physical area -means that cells from high latitudes are much smaller than cells from the tropics, which exacerbates this issue.Because the observations show that tropical savannas burn far more than any other biome, the absolute errors are highest there (Fig. 6).These regions thus likely drive most of the optimization, which could have led to the neglect of performance in, for example, the boreal region.An optimization algorithm that took relative error into account might thus improve performance in low-fire regions, while worsening it where fire is frequent.
The fact that Levenberg-Marquardt only considers the SSE of a parameter set can lead to situations as observed with Optimization 2, where the final parameter values result in functions that bear little resemblance to those derived from empirical analyses (Fig. S7).A different algorithm could penalize extreme functional forms and thus preferentially stay near more reasonable parameter values.
Simply substituting an alternative measurement for SSE in a Levenberg-Marquardt context would be less than ideal for addressing these issues.In addition to being the performance metric -i.e., the statistic by which the algorithm determines whether a parameter set has resulted in improved model performance -SSE is an inherent part of the mathematics in the Levenberg-Marquardt algorithm generating the direction and size of the step from the most recently accepted guess to the next accepted guess (Levenberg, 1944;Marquardt, 1963;Transtrum and Sethna, 2012).Using a different performance metric would still result in guesses designed to minimize SSE.This would at best reduce the efficiency of the algorithm, and at worst result in searches orthogonal to the direction of improved performance.To most effectively avoid the problems inherent with SSE, a completely different algorithm -preferably one that can use any arbitrary performance metric -would be needed.The Markov Chain Monte Carlo method (MCMC) is one such option, which has the additional benefit, as discussed above, of being a global search algorithm.It has been widely used in the Earth sciences, including by Le Page et al. (2015) to fit a global fire model.Those authors used as their performance metric a combination of (a) accuracy of classification of grid cells into burned fraction bins and (b) level of correspondence between modelsimulated and observed interannual variability.However, being a global search, MCMC requires many iterations to converge on an optimal solution -Le Page et al. (2015) reported iteration counts of hundreds to over a thousand.The deeply model-interactive setup used here -where the complete model of soil, vegetation, and fire was forced with climatic data for 19 model years -took around two hours per iteration with all 241 grid cells being run in parallel, which made MCMC and similar many-iteration algorithms computationally infeasible.
The choice of grid cells and initial conditions is also extremely important to any automated model fitting algorithm.The robustness of our results is enhanced by our use of different initial parameter set guesses (Knorr et al., 2014;Le Page et al., 2015).A more structured and informed approach to sampling grid cells for the optimization -and increasing the number of grid cells -would further improve confidence in the selected parameter set.A larger set of grid cells might have prevented Optimization 2 from traveling in the direction it did, with improved performance in the optimized grid cells but worsened performance overall.Region-specific optimizations might also be beneficial; although the general influences of different variables on fire might be consistent across biomes, vegetation structure and other factors likely mean that the relative importance of things like relative humidity or population density vary between, e.g., boreal and temperate forests, or tropical savannas in South America and Africa.
6 Conclusions: regional variations in pattern and practice FINAL represents the first attempt in a dynamic global vegetation model to separate present-day cropland, pasture, and non-agricultural burning.The importance of this can be seen, for example, in differences between pasture and nonagricultural land in the timing of the fire season -especially in Central Asia (Fig. S5).These land use/cover types also differ in fire frequency, as exhibited for example in northern Australia (Fig. 5).Overall, the combined fire model tends to perform well over much of sub-Saharan Africa, Brazil, Central Asia, and Australia.However, non-agricultural burning specifically is not well-represented in several important regions; these include eastern sub-Saharan Africa, South American savannas and grasslands, interior Australia, South and Southeast Asia, and the boreal zone (Fig. 5).The strong limitation of fire by soil moisture may have much to do with performance in those parts of the world (Sect.5.3).The apparent deficiencies of the non-agricultural fire module -the first to be tested against globally gridded estimates of nonagricultural burning -may reflect the need to more fundamentally rethink how non-agricultural fire is represented in global models.
The use of climatologies for cropland and pasture burned area is a significant limitation of FINAL.1.It allows very little interannual variability (Fig. 7) -only what results from changing agricultural area.Perhaps more importantly, however, the use of a climatology based on just nine years of observations makes it difficult to justify the use of the model very far into the past or future.Economic development can result in changes in technology, types of crops, and legislative priorities (banning crop fires, for example), all of which can affect the amount and timing of agricultural fire.Climate change has and will continue to affect the timing, length, and quality of growing seasons (Porter et al., 2014); the associated impacts on planting and harvest date will affect the timing of crop residue burning, and people will shift the timing of burns to match the shifting phenology of pasture vegetation.It is thus important to understand what information people consider in their decisions of whether, when, and how much to burn.Literature reviews and new research could shed light on traditional methods for climate forecasting based on changes in the weather and vegetation (e.g., Kagunyu et al., 2016), as well as how these cues might be tied to the timing of prescribed fire for various purposes (e.g., Laris, 2002).Advanced analytical methods could also be applied to climate and fire history observations to look for lagged, region-specific relationships of agricultural burning with weather at weekly to monthly time scales.
While temporal variation is neglected, this first version of FINAL does begin to account for regional variation in agricultural fire management practices.Other aspects of FINAL and LM3, as with many global fire and vegetation models, could be improved by representing such geographic variation.Livestock grazing intensity, as discussed above, is one important example.The shape of the population density-fire relationship also likely varies across the world.Some fire models include a spatially dependent human ignitions term (Thonicke et al., 2010) to account for this effect.Incorporating this geographic variation into FINAL could improve performance, but it would be important to do so based on independent analyses so as to avoid simply compensating for the model's errors.

Appendix A: Levenberg-Marquardt algorithm: implementation
Our implementation of the Levenberg-Marquardt algorithm (Fig. A1) began with a Bash script that set up the files and directories necessary to run the fire model at the 241 points.These points would then be run for 1991-2009 in parallel.Once this first iteration was complete, a Python script calculated the sum of squared errors (S) over each grid cell (c), year (y), and month (m): Here, E refers to the model-estimated burned area, and O refers to an observation-based estimate of burned area.Specifically, we focused on non-agricultural lands, using as our "observations" estimates generated for each month and year by the method detailed in Rabin et al. (2015) but with F k estimates restricted to non-negative values.The Python script then generated a new parameter set guess based on the initial values of the parameters and saved a flag telling the Bash script to run the model again with the new guess.
After this and subsequent model runs, another Python script would calculate the associated value of the sum of squared errors (S t ) and compare it to the sum of squared errors from the most recently accepted guess (S * ).If S t < S * , the current parameter set guess (β t ) would be "accepted" and become the new value of β * , and λ would be decreased.Otherwise, β t would be "rejected," with β * retaining its previous value, and λ being decreased.In either case, a new guess would then be generated based on β * and the new value of λ, the model would be run again, and the process would repeat (Fig. A1).
The Python script we developed was based on a MATLAB routine for Levenberg-Marquardt solutions of nonlinear least squares problems called marquardt.m(Nielsen, 2001), further documented in Nielsen (1999).Besides porting it to Python, we made a number of changes to the original code.Some restructuring was related to the fact that the new parameter sets could not be evaluated within Python.Others were to incorporate new features, such as the limited multiplicative damping based on work by Transtrum and Sethna (2012).Nielsen (2001) uses a somewhat complex method to update δ after every each iteration (Fig. A2).If S t ≥ S * , λ is multiplied by a value ν, whose initial value is 2 and is doubled after every rejected guess.If a guess is accepted (S t < S * ), ν is reset to 2, and λ is decreased.We made some changes to the original code as a result of the aforementioned restructuring, with λ being reduced as

Model& run& t
. Method for updating λ, after Nielsen (1999Nielsen ( , 2001)). where Note that there have been many methods proposed over the years for updating the damping parameter in the Levenberg- Marquardt algorithm.These impact the size of the steps the algorithm takes while searching through parameter space, with implications for efficiency.However, the math by which the algorithm determines which direction on each dimension to move is unaffected.
The algorithm has several possible stop conditions.We set a maximum of 300 iterations, which was never reached.The algorithm would also stop if the Python script detected that the gradient was decreasing very slowly: if the step size was very small: ) or there was an issue of near-singularity in one of the matrices involved in solving for the new parameter step: where is the smallest number allowed by the numerical precision of the Python environment.However, in practice, we usually ended up halting the algorithm manually.Each iteration took about two hours, and once we noticed neither the sum of squared errors nor any parameter changing by very much, we would stop the runs.This could have been avoided by choosing more appropriate threshold values for the stop conditions, but likely did not appreciably impact the results.We initially selected 250 land cells at random from the LM3 grid, but rejected 9 for various reasons (all glacier, all lake, etc.).This left us with 241 grid cells which we would use for the optimization.Preliminary tests, however, revealed a few problems with the selection: a bias towards improving model fit in grid cells with strong model underestimation was evident (i.e, grid cells where the model simulated too much fire were undersampled), and the high northern latitudes -which make up a small fraction of global land area and an extremely small fraction of global fire activitywere judged to be oversampled.We got rid of 14 of those far northern grid cells (from Greenland and the Canadian tundra), then selected 23 new cells to bring us up to 250.The new cells were specifically selected from cells where a preliminary model run either underestimated or overestimated non-agricultural burned area relative to the unpacked data.Unfortunately, the model's performance in that preliminary run did not well match how the model actually performed in our optimization run.As such, we ended up oversampling areas of underestimation, leading to a bias towards making the model burn too much.We then culled the most extreme underestimated grid cells one by one until the sums of squared errors from underestimated and overestimated grid cells generated by the initial guess were approximately equal.This left us again with 241 grid cells, whose locations and initial sum of squared errors are shown in Fig. A3a.A histogram of the mean annual error in burned area of the initial guess (Fig. A3b) shows that the positive and negative errors in this new dataset are approximately balanced.Although the global amount of grazed vegetation seems to have been simulated well (as discussed above), much variation likely exists among regions in how intensely land is grazed.This is not captured by the assumption in our model of a 4 % daily grazing rate.Combustion completeness values being too low would also lead to too-high estimates of aboveground biomass, but the possible effect of this on estimated emissions is unclear.Increasing combustion completeness would increase fire emissions in the short term, but as any individual pasture tile grew older and approached equilibrium biomass, fire emissions might be no different.That is, decreased biomass with increased combustion completeness might not change emissions density.
The amount of leafy vegetation on pastures -not just that consumed -also appears to have been simulated well.On average over 2001-2009, FINAL.1 simulated 3.4 kgC m −2 of aboveground biomass on pastures, including both live vegetation and dead material.This was broken down into live leaves (0.22 kgC m −2 ), live stems (0.94 kgC m −2 ), leaf litter (0.45 kgC m −2 ), and dead woody material (1.8 kgC m −2 ); these pools are mapped for the world's major pasture regions in Fig. A4.In their work in the Waikato region of New Zealand -a moist, temperate ecosystem dominated by C3 grasses - Hanna et al. (1999) defined active pastures as containing no more than 0.2 kgC m −2 of live leaves or 0.15 kgC m −2 of dead material.FINAL.1 simulated less than 0.1 kgC m −2 of live leaf tissue in New Zealand, and indeed the world's temperate pastures seem to satisfy the ≤ 0.2 kgC m −2 criterion (Fig. A4a).The tropics generally see much higher modeled pasture leaf biomass; in all cases, leaf biomass does not much exceed 0.25 kgC m −2 (Fig. A4a).Uhl and Kauffman (1990) describe a pasture in eastern Amazonia with 0.6 kgC m −2 of nonwoody material; this is close to the simulated value of combined live and dead leaf C (Fig. A4a,  c).Kauffman and Cummings (1998), looking at three other pastures in Amazonia, found a range of 0.8-1.5 kgC m −2 of fine fuels, which included both live and dead leaf material as well as fine woody debris.Again, this corresponds well with our results (Fig. A4a, c), although we do not simulate fine woody debris.Kauffman and Cummings (1998) also found 1.3-5.2kgC m −2 of large downed trunks remaining from the initial clearance of forest for pasture; the simulation produces levels of woody litter in that range for pastures in the Atlantic Forest region of Brazil and in southern China (Fig. A4d).Savadogo et al. (2007) found a mean of 0.045 kgC m −2 of live biomass in the tree and bush savanna of Burkina Faso S. S. Rabin et al.: A global fire model with agricultural burning practices -a value similar to that of the combined live leaf and stem pools simulated by LM3 (Fig. A4a, b).
However, LM3 does seem to have overestimated pasture biomass in tropical savanna regions when including dead woody material.In that same part of Burkina Faso, Savadogo et al. (2007) found a mean of 0.07 kgC m −2 of dead material, whereas LM3 simulated values of around 0.2-0.3kgC m −2 (Fig. A4c, d).Ottmar et al. (2001) found that land in the Cerrado with a significant herbaceous layer (campo limpo, campo sujo, and cerrado ralo) generally tended to have less than 1 kgC m −2 of aboveground live and dead biomass; LM3 simulated about 1-1.5 kgC m −2 (Fig. A4e).It is not clear whether the sites examined by Ottmar et al. (2001) were actively grazed; if not, pastures there would be expected to have even less biomass, in which case LM3's overestimate would be more pronounced.
The fact that FINAL does not explicitly simulate fire associated with land clearance likely contributes to its overestimation of pasture (and cropland) fire emissions density.In the version of LM3 used here, biomass killed during land use transitions can be either harvested or wasted.Harvested wood biomass goes to one of three long-lived virtual emissions pools, while wasted biomass is transferred to litter.But in reality, wood remaining after harvest (also known as slash) is often burned, especially in the high-biomass moist tropical forest biome.The emissions involved are significant: tropical deforestation burns were estimated by van der Werf et al. (2010) to contribute up to 15 % of global annual fire CO 2 -C emissions on average.Instead of breaking this out into a separate flux, LM3 and FINAL are conflating land clearance fire emissions with the emissions from subsequent burning of the cleared land for agricultural management.This is unfortunately not a mere accounting quirk; the use of one or two burns to get rid of most of the remaining slash wood means that fire emissions drop rapidly a few years after land clearance, whereas LM3 and FINAL simulate a gradual decrease over time.

Figure 1 .
Figure 1.Mean fractional land cover of (a) non-agricultural land, (b) cropland, and (c) pasture over 2001-2009 as simulated in model runs (after Hurtt et al., 2011).Gray cells did not contain any of the indicated land cover type.

Figure 2 .
Figure 2. Trace plots showing the progression of sum of squared errors (SSE) (a), the percent change in SSE between each accepted parameter set (b), and each of the ten parameters (c-l) over the length of the optimization.x axes show iteration number, y axes show sum of squared errors or parameter guess value, and color of points indicate whether the associated parameter set guess was accepted (blue) or rejected (red).

Figure 3 .
Figure3.Optimization 3: changes in functions that were optimized, from originalLi et al. (2012Li et al. ( , 2013) ) functions (solid gray), to initial guesses with Gompertz-style functions where necessary (dashed red), to final parameter set (solid blue).Color bar in panel (f) indicates difference in the cubed product of f θ and f RH (range 0-1) between the original and new parameterizations, with blue indicating a lower value in the new parameterization.

Figure 4 .
Figure 4. Change in non-agricultural fire model performance for Optimization 3. (a-b) Mean annual burned fraction on non-agricultural lands from the initial guess (a) and the final parameter set (b; identical to Fig. 5i) (c-d) Difference between runs FINAL.0 and FINAL.1 in correspondence of modeled to unpacked non-agricultural burning (Fig. 5e) as measured by mean annual burned fraction (c) and sum of squared errors of burned area evaluated at monthly resolution (d).For (c) and (d), blue indicates improvement by FINAL.1 over FINAL.0.

Figure 7 .
Figure 7. Annual time series of observed and model-estimated burned area (a) (km 2 ), and fire carbon emissions (b) (PgC yr −1 ), from 2001 to 2009.Dashed lines: observational estimates of total and by-land cover fire emissions fromRabin et al. (2015).Solid blue lines: observations of total emissions from GFED3s(Randerson et al., 2012).Other solid lines: model-estimated total and by-land cover fire emissions.

Figure 9 .
Figure 9. Difference in non-agricultural fire model burned area sum of monthly squared errors (SSE) between Optimizations 2 and 3. Blue represents areas where the latter performs better.

Figure A3 .
Figure A3.Summary of performance of Optimization 1 initial guess in grid cells chosen for optimization with regard to non-agricultural burning.(a) Map of sum of squared errors.(b) Histogram of error in mean annual burned area.
Rabin et al.:A global fire model with agricultural burning practices density of anthropogenic ignitions per time step.(Henceforth, β will denote parameters determined during our optimization routine as described in Sect.2.6.The final values of these parameters can be found in Table With β I a representing the rate of ignitions per person at each time step and P D representing population density (people km −2 ), the first part of Eq. (3) gives a starting value for www.geosci-model-dev.net/11/815/2018/Geosci.Model Dev., 11, 815-842, 2018 S. S.

Table 1 .
Combustion completeness and mortality values for each "species" and tissue pool.Note that "stem" refers to both aboveground and belowground stem biomass, and that "root" refers only to fine roots.n/a -not applicable.

Table 2 .
Experimental runs performed in this study."1-300" indicates that 300 years were simulated, but that these are not tied to any historical data, and thus they do not correspond to any actual historical years."Graze rate" refers to the amount of non-wasted leaf biomass consumed each day by livestock on pasture tiles.n/a -not applicable Name Rabin et al.: A global fire model with agricultural burning practices default setting for CORPSE is to simulate 15 different belowground soil cohorts (age classes); to improve computational efficiency, we set CORPSE to simulate only one.

Table 3 .
Initial and final parameter sets for each optimization.Here, values are rounded to nearest 10 −4 ; full-precision values can be found in TableS1in the Supplement.Optimization 1 did not complete successfully.
Table 3, was run from 1948 to 2009 (FINAL.1;Table2).This run began with initial conditions as produced for the beginning of 1948 by the original LM3 run described above (LM3_ORIG).An experimental run with the complete new model structure but all set- www.geosci-model-dev.net/11/815/2018/Geosci.Model Dev., 11, 815-842, 2018 Optimization 3. A discussion of what the results imply for LM3 and fire modeling generally can be found in Sect.5.3.Figures in the Supplement illustrate the difference in sum of squared errors in optimized grid cells for model runs with the initial and final parameter sets from Optimization 3 (Fig.