Improving the LPJmL4-SPITFIRE vegetation–fire model for South America using satellite data

Vegetation fires influence global vegetation distribution, ecosystem functioning, and global carbon cycling. Specifically in South America, changes in fire occurrence together with land-use change accelerate ecosystem fragmentation and increase the vulnerability of tropical forests and savannas to climate change. Dynamic global vegetation models (DGVMs) are valuable tools to estimate the effects of fire on ecosystem functioning and carbon cycling under future climate changes. However, most fire-enabled DGVMs have problems in capturing the magnitude, spatial patterns, and temporal dynamics of burned area as observed by satellites. As fire is controlled by the interplay of weather conditions, vegetation properties, and human activities, fire modules in DGVMs can be improved in various aspects. In this study we focus on improving the controls of climate and hence fuel moisture content on fire danger in the LPJmL4-SPITFIRE DGVM in South America, especially for the Brazilian fireprone biomes of Caatinga and Cerrado. We therefore test two alternative model formulations (standard Nesterov Index and a newly implemented water vapor pressure deficit) for climate effects on fire danger within a formal model–data integration setup where we estimate model parameters against satellite datasets of burned area (GFED4) and aboveground biomass of trees. Our results show that the optimized model improves the representation of spatial patterns and the seasonal to interannual dynamics of burned area especially in the Cerrado and Caatinga regions. In addition, the model improves the simulation of aboveground biomass and the spatial distribution of plant functional types (PFTs). We obtained the best results by using the water vapor pressure deficit (VPD) for the calculation of fire danger. The VPD includes, in comparison to the Nesterov Index, a representation of the air humidity and the vegetation density. This work shows the successful application of a systematic model–data integration setup, as well as the integration of a new fire danger formulation, in order to optimize a process-based fire-enabled DGVM. It further highlights the potential of this approach to achieve a new level of accuracy in comprehensive global fire modeling and prediction.


Introduction
Fire in the Earth system is an important disturbance leading to many changes in the vegetation and has substantial impact on biodiversity, human health and ecosystems (Langmann et al., 2009). Fire is responsible for ca. 2 Pg of carbon emissions, which constitutes 20 % of global carbon emissions (Giglio et al., 2013;Werf et al., 2010). Fire-induced aerosol emissions and land surface changes modify evapotranspiration and surface albedo and have therefore a crucial impact on global climate (van der 5 Werf et al., 2008;Yue and Unger, 2018). Despite a tendency for globally declining burnt area (Andela et al., 2017;Forkel et al., 2019b), more frequent and intense drought-periods lead to increasing fire-prone weather and surface conditions worldwide and therefore fire danger (Jolly et al., 2015). Growing fire danger along with land-use change are increasing ecosystem's vulnerability, which could in turn shift entire regions into a less vegetated state (Silvério et al., 2013). To account for these effects, it is extremely important to include well performing fire modules in Dynamic Global Vegetation Models (DGVMs). 10 Especially in South America, tropical forests, woodlands and other ecosystems are vulnerable to increasing fire danger and land use change (Cochrane and Laurance, 2008). This study focuses on the fire behavior in central-northern South America and especially on the Brazilian biomes Caatinga and Cerrado, which are the most fire-prone regions in South America (Fig. 1).
Together with the Amazon rainforest they form an area of very high biodiversity and have a large impact on the global carboncycle and the regional water cycle (Lahsen et al., 2016). Compared to the Amazon, the Cerrado and Caatinga are both less  Giglio et al., 2013) and the biomes Amazonia, Cerrado and Caatinga (IBGE, 2019;Harvard, 2019) 15 densely vegetated and drier biomes, but with very different vegetation and precipitation dynamics. The Cerrado is a savannalike biome with a mixture of shrubs, high grasses and dry forest parts. With a precipitation of ca. 1500 mm per year the Cerrado 1987), the Keetch Byram Drought Index (Keetch and Byram, 1968), the Angström Fire Danger Index (Arpaci et al., 2013), and the Nesterov index (Venevsky et al., 2002)). However, regional studies show that fire weather indices tend to have different predictive performances for fire occurrence (Arpaci et al., 2013). Hence, the performance of different fire weather indices should be ideally tested in order to accurately represent fire danger in DGVMs. However, not all fire weather indices can be easily adapted for global fire models because they require input variables that are not available within a DGVM framework. 5 Hence a fire danger index for a DGVM should be as complex as necessary but still relatively easy to implement. As a result, the relatively simple Nesterov index has been widely used within global fire models (Venevsky et al., 2002;Thonicke et al., 2010).
Here, we aim to improve the simulated occurrence of fire (i.e. burnt area) in the LPJmL4-SPITFIRE model for South America and in particular for the fire-prone biomes Cerrado and Caatinga. We aim to evaluate the performance of two alternative fire 10 danger indices within SPITFIRE based on the already implemented Nesterov index (Venevsky et al., 2002) and the newly implemented water vapor pressure deficit (VPD thereafter, Pechony and Shindell, 2009;Ray et al., 2005). Furthermore, we apply a formal model-data integration framework (LPJmLmdi, Forkel et al., 2014) to estimate model parameters that control fire danger, fire behavior, fire resistance and mortality against satellite-based data sets of burnt area and above-ground biomass (Fig. 2). Our approach is likely to improve the representation of spatial-temporal variations in fire behavior in different biomes 15 to enable a much better modelling of the impact of climate change on fire-vegetation interactions in the current century.  Schaphoff et al., 2018a, b), is a well established and validated 20 process-based DGVM, which globally simulates the surface energy balance, water fluxes, carbon fluxes and stocks, and natural and managed vegetation from climate and soil input data. LPJmL simulates global vegetation distribution as the fractional coverage of plant functional types (PFT), which is called foliage projective cover (FPC), and managed land as fractional coverage of crop functional types (CFT). The establishment and survival of different PFTs is regulated through bioclimatic limits and effects of heat, productivity and fire on plant mortality. Therefore, it enables LPJmL to investigate feedbacks, for exam- 25 ple between vegetation and fire. In standard settings, which are also used here, the model operates on the grid of 0.5 • × 0.5 • latitude-longitude with a spinup time of 5000 years, repeating the first 30 years of the given climate data set.
Since its original implementation by Sitch et al. (2003), LPJmL has been improved by a representation of the water balance (Gerten et al., 2004), a representation of the agriculture (Bondeau et al., 2007), and new modules for fire (Thonicke et al., 2010), permafrost (Schaphoff et al., 2013) and phenology (Forkel et al., 2014).  Thonicke et al., 2010) is a process-based fire module, used in various vegetation models (e.g. Lasslop et al., 2014;Yue and Unger, 2018), including LPJmL4. We describe here its main features, which are published in Thonicke et al. (2010). SPITFIRE calculates fire disturbance by simulating the ignition, the danger, the spread and the effects of fire separately. As ignition sources SPITFIRE considers human ignition and lightning flashes. Human ignitions 5 (n h,ig ) are calculated as a function of population density: where k(P D ) = 30.0 · exp(p h · P D ).
(2) P D is the human population density (individuals km −2 ) and a(N D ) (ignitions individual −1 day −1 ) describes the inclination 10 of humans to cause fire ignitions (Eq. 3 and 4 in Thonicke et al., 2010). p h is a parameter, which is set to -0.5 in Thonicke et al. (2010). This relationship assumes that human ignitions are lowest on very low populated regions and on high populated regions through a higher level of urbanization and landscape fragmentation. The ignition is highest for a medium-small population density. Lightning-caused ignitions are prescribed by lightning data from the OTD/LIS Gridded Climatological data set (Christian et al., 2003), assuming that 20 % of the flashes reach the ground and 4 % of cloud-to-ground strikes can start a fire.
In the study area of South America human ignitions are by far the most dominant ignition source, due to missing lightning in the dry season.
Fire danger is by default computed by using the Nesterov index which accounts for the maximum and dew point temperatures as well as scaling factors for different PFTs on a daily time step. In the following section, we describe the calculation of the 5 fire danger indices in detail (Sect. 2.2). Fire duration t f ire (min) is calculated as a function of the fire danger index, assuming that fires burn longer under a high fire danger: where p t is set to -11.06 in Thonicke et al. (2010). The maximum fire duration per day is 240 minutes.
The calculation of the forward rate of spread ROS f,surf ace (m min −1 ) is based on the Rothermel equations (Rothermel, 1972;10 Pyne et al., 1996;Wilson, 1982): where I R is the reaction intensity, ξ the propagation flux ratio, Φ w a multiplier that accounts for the effect of wind, the effective heating number, Q ig the heat of pre-ignition and ρ b the fuel bulk density (Eq. 9 in Thonicke et al., 2010). ρ b (kg m −3 ) is a PFT-dependent parameter and describes the density of the fuel, which is available for burning. It is weighted over 15 the different fuel classes. Hence, a changing PFT distribution has an impact on ROS f,surf ace .
The simulated fire ignitions, fire danger and fire spread are then used to calculate the burnt area, fire carbon emissions, and plant mortality. Plant mortality depends on the scorch height (SH) and the probability of mortality due to crown damage P m (CK).
SH describes the height of the flame at which canopy scorching occurs. It increases with the 2/3 power of the surface intensity I surf ace : where F is a PFT-dependent parameter. Assuming a cylindrical crown, the proportion CK affected by fire is calculated as: where H is the height of the average woody PFT and CL the crown length. The probability of mortality P m (CK) due to crown damage is then calculated by: where rCK is a PFT depended resistance factor between 0 and 1, and p in the range of 3 to 4. Disturbance by fire mortality has a large impact on the vegetation dynamics, which are calculated within LPJmL. SPITFIRE further includes a surface intensity threshold (10 6 , fraction burnt area per grid cell), which describes the threshold of the possible area burnt below which the surface intensity is set to zero and hence burnt area, emissions and fuel consumption is set to zero.

30
SPITFIRE considers anthropogenic effects on fire by taking into account human ignitions but does not account for fire suppression. Only wildfires occurring in natural vegetation are simulated. Fire on managed land like agriculture or pasture areas is not implemented, which has to be taken into account if simulated burnt area is compared with satellite observation.
Furthermore, we introduced a small technical change in the LPJmL4 interaction with SPITFIRE compared to the original SPITFIRE implementation: In the version 4.0 of LPJmL the fire litter routine calculates the leaf and litter carbon pools in a 5 daily time step. Since the LPJmL tree allocation works at a yearly time step, this implementation leads to an incorrect LPJmL4-SPITFIRE interaction. We now split the fire-litter routine into two parts; the first one allocates burnt matter into the litter at a daily time step without recalculating the pools and the second one calculates the leaf and root carbon pools at a yearly time step.

Fire danger indices 10
The fire danger index (FDI) is a key parameter within process-based fire models such as SPITFIRE. The FDI determines the probability and the intensity of a spreading fire, which impacts fire behavior.

Nesterov index-based fire danger index (F DI N I )
The fire danger index within SPITFIRE is based on the daily (d) calculated Nesterov Index NI(d) (Venevsky et al., 2002), which is widely used in numeric fire simulations. The NI is a cumulative function of daily maximum temperature T max (d)( • C) and 15 dew-point temperature T dew (d)( • C) and set to zero at a precipitation ≥ 3 mm or a temperature ≤ 4 • C: The resulting fire danger index has been calculated as in Schaphoff et al. (2018a) (slightly different compared to Thonicke et al. 20 (2010)) by taking into account the NI as measure for weather conditions and a PFT-dependent scaling factor α N Ii : where n is the number of PFTs and m e the moisture of extinction, which is a PFT-dependent parameter and is weighted over the litter amount. We will use the scaling factors α N Ii in the parameter optimization (Sect. 2.4). We implemented a new fire danger index, based on the water vapor pressure deficit (VPD). The VPD describes the difference of the saturation water vapor pressure e s and the actual water vapor pressure in the air. For the parameterization of the VPD we used an approach based on Pechony and Shindell (2009): where T is the air temperature, RH the relative humidity and Z the Goff-Gratch equation (Goff and Gratch, 1946) to calculate the saturation vapor pressure. The flammability F at time step t for each grid cell can then be expressed as: where VD is the vegetation density, R the total precipitation in mm/day and c R is a constant factor (c R = 2 day/mm). Here we used the simulated FPC from LPJmL4 as a proxy for the VD. The soil is a natural buffer for drought periods and heavy rainfall events. In the Nesterov index this was taken into account by the cumulative nature of this index. Since the VPD-based fire danger index is not cumulative, this buffering effect is taken into account by taking the monthly mean of the precipitation.
In doing so we avoid unrealistic high flammability fluctuations in time steps with isolated events of very low or very high 10 precipitation (R).
Based on this implementation in SPITFIRE, the resulting FDI was much smaller than the original F DI N I . Hence, we scaled the VPD up with a PFT-dependent scaling factor α V P Di , weighted over the corresponding FPC: α V P Di for the F DI V P D was not included in (Pechony and Shindell, 2009), but is important in order to allow different fire 15 responses for different tree and grass types. We will use the scaling factors α V P Di in the parameter optimization (Sect. 2.4).
In comparison to the NI, the F DI V P D requires more climate variables as input as it uses relative humidity and vegetation cover as additional fire-relevant variables. Vegetation cover has a direct link to fire risk by providing the number of available fuel for burning. According to many studies (e.g. Ray et al., 2005;Sedano and Randerson, 2014;Seager et al., 2015) the F DI V P D is a very accurate fire danger index with a high correlation with fire occurrence, while still being relatively easy to implement in 20 a global fire model.
The general behavior of the two indices as modelled by LPJmL in dependence of relative humidity and temperature is shown in Fig. 3. The Nesterov index shows a strong but very localized maximum for high temperatures and a small humidity. Hence a spreading fire is only possible in a very small climate range (here ca. from 25°Celsius and a relative humidity smaller than 0.5).
The VPD on the other hand shows a less pronounced maximum but a medium fire danger also for wetter and colder regions. 25 The slope of towards lower VPD values is also smaller compared to the Nesterov index. Especially in regions with temperatures colder than 20°C and a relative humidity smaller than ca. 0.6 a fire is still possible. This might increase the area in which fires can occur compared to the Nesterov index, which could be an important improvement, enabling SPITFIRE to simulate more fire in wetter and colder regions. The calculated VPD and NI values shown in Fig. 3 are based on a LPJmL-SPITFIRE run, and thus the influence of vegetation distribution on both fire danger indices.

Model input data
LPJmL4-SPITFIRE requires input data on daily air temperature, precipitation, long-wave and shortwave downward radiation, wind and specific humidity, which are taken from the NOAH Global Land Assimilation System (GLDAS, Rodell et al., 2004).
The data has a spatial resolution of 0.25 • × 0.25 • and the time step is 3 hours. We regridded and aggregated the data set to the LPJmL resolution of 0.5 • × 0.5 • and to a daily time step. We used the GLDAS 2.0 for the years 1948-1999 and the ver- Furthermore, LPJmL4-SPITFIRE is forced with gridded constant soil texture (Nachtergaele et al., 2009) and annual information on land use from Fader et al. (2010). Atmospheric CO 2 concentrations are used from Mauna Loa station (Quéré et al., 2015) and applied globally. The population density is taken from Goldewijk et al. (2011) and the lightning flashes are taken from the OTD/LIS satellite product (Christian et al., 2003).

Model optimization
To estimate parameters of LPJmL4-SPITFIRE, we aimed to calibrate model results against satellite observations of burnt area (GFED4: Giglio et al., 2013;Werf et al., 2017). However, as fire occurrence and spread impact and depend on vegetation productivity, hence fuel load, we wanted to ensure to not over-fit LPJmL4 against burnt area but to additionally achieve a realistic vegetation distribution. Therefore, we additionally included a satellite-derived data set on above-ground biomass of 10 trees (AGB, Avitabile et al., 2016) in the optimization. We combined burnt area and AGB with the corresponding model outputs within a joint cost function and applied a genetic optimization algorithm to estimate model parameters (Fig. 2). The implementation of the genetic optimization algorithm (GENetic Optimization Using Derivatives (GENOUD), Mebane and Sekhon, 2011) for LPJmL is described in Forkel et al. (2014). The used cost function is based on the Kling-Gupta efficiency (KGE), which is the Euclidean distance in a three-dimensional space of model performance measures that account for the bias, 15 ratio of variance and correlation between simulations and the observations. Gupta et al. (2009) showed that the KGE performs in an optimization setup is better than, e.g., the Nash-Sutcliffe efficiency (and hence MSE). We extended the KGE by defining it for multiple data sets d (i.e. burnt area and AGB): Several parameters of LPJmL4-SPITFIRE were included in the optimization that cover different fire processes (see Tab. 2): ignition (human ignition parameter p h , Eq. 2), fire danger (scaling factors FDI (α N Ii and α V P Di ), Eq. 10 and 13), fire spread (fire duration p t , Eq. 3), fuel bulk density (ρ b , Eq. 4), surface intensity threshold and fire effects (scorch height parameter F, Eq. 5; crown mortality parameter rCK, Eq. 7). While p t , p h and the surface intensity threshold are global parameters (for all PFTs), the others were optimized for each PFT separately. Since we focus here on tropical South America, we used tropical broadleaved evergreen (TrBE), tropical broadleaved raingreen (TrBR) and tropical herbaceous (TrH) PFTs for the optimization.
In genetic optimization algorithms, each model parameter is called an individual with a corresponding fitness, which represents 5 the cost of the model against the observations. At the beginning of the optimization process, the GENOUD algorithm creates a generation of individuals based on random sampling of parameter sets within the prescribed parameter ranges. After the calculation of the cost of all individuals of the first generation, a next generation is generated by cloning the best individuals, by mutating the genes or by crossing different individuals (Mebane and Sekhon, 2011). This results, after some generations, in a set of individuals with highest fitness, i.e. parameter sets with minimized cost. To find an optimum parameter set also used 10 the BFGS gradient search algorithm (named after the authors Broyden, 1970;Fletcher, 1970;Goldfarb, 1970;Shanno, 1970) within the GENOUD algorithm. An optimized parameter set of the BFGS algorithm is used as individual in the next generation.
We were applying the GENOUD algorithm with 20 generations and a population size of 800 individuals per generation, which corresponds to 16000 single model runs. We decided on this amount of iterations, because the cost kept almost constant in the last iterations and the parameter values did not change to the 6th digit, beyond which changes are not really relevant for model 15 applications. During the optimization we ran the model parallel for each grid cell (40 grid cells and CPU's, 3.2GHz) and had a total optimization time of ca. 24 hours.
The comparison of the two presented fire danger indices is the main objective of this study. Hence the optimization of the PFT-dependent FDI scaling factors α N Ii and α V P Di is crucial and obligatory for the VPD because of no prior values. Accordingly, we conducted two different optimization experiments using LPJmLmdi: First, using a FDI based on the VPD (VPD optim 20 hearafter) and secondly using the a FDI based on the NI (NI optim hereafter). Both resulting parameter sets were then used for LPJmL4 runs and were compared to the unoptimized original model version using the NI (NI orig hereafter) and various evaluation data sets. 25 We used burnt area from the global fire emission database (GFED4; Giglio et al., 2013;Werf et al., 2017) (2016). This data set is approximately representative for the late 2000s and therefore we compared it against simulated AGB for the years 2009-2011. We regridded all data set to a 0.5 • × 0.5 • resolution. In addition, we used maps of PFTs as derived from the ESA CCI land cover map V2.0.7 (Li et al., 2018;Forkel et al., 2014).

Evaluation metrics
To quantify the performance of the model output, we applied the Pearson Correlation between two time series, the normalized 5 mean square error (NMSE; Kelley et al., 2013) and the Willmott coefficient of agreement (W; Willmott, 1982) to describe differences between the model simulation and the reference data sets. The NMSE is calculated by: where y i is the simulated and x i the observed value in the grid cell i.x is the mean observed value. The NMSE is zero for perfect agreement between simulated and modelled results, 1.0 if the model is as good as using the observed mean as a 10 predictor and larger than 1.0 if the model performs worse than that. We chose the NMSE to represent and compare the model errors, as it has a squared error term, which puts a stronger emphasis on large deviations between simulations and observations as compared to a linear term, and due to its normalization it is comparable across different parameters. Especially for fire simulations we have a relatively large deviation between simulations and observations. The Willmott coefficient of agreement is given by: which additionally accounts for the area weight A i of the grid cell i. The Willmott coefficient is a squared index, where a value of 1 stands for perfect agreement between simulated and modelled runs and gets smaller for worse agreements with a minimum of 0. Unlike the coefficient of determination, the Willmott coefficient is additionally sensitive to biases between simulations and observations. 20 3 Results

Performance of optimized fire danger index formulations
Overall, the yearly burnt area simulated by the standard SPITFIRE model (using the original Nesterov index, NI orig ) showed poor simulation results over South America as compared to the GFED4 evaluation data set ( Fig. 4 a and    The optimized version using NI optim (Fig. 4 d), led to an overall decrease of fire, with a slight improvement of NMSE (1.09) as compared to NI orig and a worse Willmott coefficient of 0.08. While the overestimation of fire in Caatinga was reduced, all the fires across South America also have decreased significantly, which led to a general underestimation of fire by 90 % (2 million 5 ha). The optimized version, using VPD optim (Fig. 4 c), clearly improved the model performance, mainly by shifting much of the simulated burnt area from the sparsely vegetated Caatinga towards the Cerrado region (NMSE=0.82 and W=0.56). In addition, by using VPD optim , the model results also showed fire occurrence in northern South America, where fire was not at all or only minimally simulated when using NI optim or NI orig . The total burnt area was in this model version ca 20 % smaller than the evaluation data set (16 million ha).

10
The burnt area time series from 2005 to 2015 provides a more detailed view on the model performance for the fire-prone Cerrado and Caatinga region (Fig. 5)  Overall, the total amount of burnt area in the Cerrado was for all three model versions smaller than in the evaluation data set. Fire occurrence in the Caatinga was, on the other hand, largely overestimated by the NI orig and the VPD optim version.
Just in the NI optim version the burnt area of the Caatinga is in the same order of magnitude as the evaluation data set, which also led, however, to a large underestimation in the Cerrado (Tab. 1 and Fig. 4). Also the Amazonia region mostly improved by using the VPD optim version (Tab. 1, Fig. A3). The R 2 and the Willmott coefficient improved, while the NMSE increased 5 slightly. With the Nesterov index fire was strongly underestimated in the Amazon region, while the optimized VPD fixes this underestimation. The fire is only modelled (and also observed, see Fig. 4) at the edges to the Amazon, where wood density is lower and deforestation takes place. In the closed continuous forest area towards the center of the Amazon almost no fire is observed and neither modelled. The total burnt area increased from 0.7 million ha to 4.8 million ha (for VPD optim ) , which is now a bit overestimated to the observed burnt area of 3.4 million ha. Using the NI optim all error metrics as well as the total   In case of the other optimized parameters the boundaries were set smaller in order to decrease the possibility that a large error in the estimation of several parameters would lead to a better overall cost in the optimization procedure. The human ignition parameter became smaller for both optimizations, which led to a smaller amount of human ignitions (from -0.5 to -0.54 in NI optim and -0.53 in VPD optim ). The fuel bulk density increased for all three tropical PFTs in the NI optim version, while for VPD optim the fuel bulk density of the TrBE and TrH PFTs decreased and for the TrBR increased. For the NI optim 5 version, the fire duration parameter (p t ) increased, leading to a shorter fire duration (from -11.06 to -9), while the value for the VPD optim version stayed relatively similar (-11.37) to the prior value. The surface intensity threshold became slightly larger for the NI optim version than the original value (from 10 −6 to 1.03·10 −6 ). For VPD optim the parameter increased by a factor of three (3.63 · 10 −6 ). The mortality related parameters F and rCK led in the NI optim version both to a decrease in the fire-related mortality for TrBE and an increase for TrBR PFTs. The optimized parameters for VPD optim led to a decrease in the fire-related 10 mortality for both PFTs except for the TrBR rCK, which led to an increased mortality.
The relative uncertainties were for most optimized parameters very small (between 0 and 10%), hence these parameters were strongly constrained (Fig. 6). Just the fire-mortality related parameters (F and rCK) had large uncertainties for the TrBR, hence were weakly constrained. For VPD optim the uncertainty of rCK (TrBR) was 0.8 and for NI optim the uncertainty of Fwas 0.9 and for rCK 1 (TrBR).
The decrease in the model error (cost) due to the optimization process has been mainly due to improvement in the burnt area.
While for the NI optim the cost of the burnt area dataset improved by 81%, the cost of the biomass dataset improved just by 6%.
In case of the VPD optim the cost of the burnt area dataset improved by 49%, whereas the biomass dataset improved by 19% 5 (Fig. A5).

Model evaluation for South America
The modelled above-ground biomass (AGB) of trees in South America was throughout all model versions larger than the evaluation data set indicates (Fig. 7). Especially the biomass in the Amazon region is with an average of ca. 20 kgC/m 2 about one third overestimated. The drier savanna regions on the continent yielded a biomass of ca. 5-10 kgC/m 2 , which also The modelled foliage projective cover (FPC) showed for all three model versions a strong underestimation compared to the evaluation data set of the TrBE throughout the whole Amazonian region (ca. 50% compared to ca. 100% in the evaluation dateset). In the fire-prone biomes Cerrado and Caatinga, however, the TrBE PFT was sometimes overestimated (TrBE cover 15 between 0 and 40 %, Fig. 8). In the regions with less TrBE the dominant PFT was mostly TrBR (Cerrado) or TrH (Caatinga) (see Fig. A1 and A2). The VPD optim version, on the other hand, showed an slightly improved TrBE distribution (NMSE=0.41, W=0.82) but also in this case we obtained an even larger improvement, when only the fire-prone regions Cerrado or Caatinga are considered (Tab. 3). Also for the TrBR and TrH PFT distributions the optimization lead to an improved performance using the VPD optim 5 in the Caatinga and Cerrado, whereas the PFT distribution in the Amazon remained similar to the prior PFT distribution. In the NI optim version, parameter optimization only slightly reduced TrBR cover showing a worse performance compared to VPD optim . However, herbaceous cover changed only slightly in all optimization experiments ( Fig. A1 and A2).

Discussion
In summary, our results show that the implementation of a new fire danger index based on the water vapor pressure deficit 10 F DI V P D and its optimization against satellite data sets improved the simulations of fire in LPJmL4-SPITFIRE, both in terms of spatial patterns as well as temporal dynamics of burnt area. In the following, we discuss the model improvements, limitations and recommendations for future improvements of process-based global fire models within the DGVM framework.

Improvements in model performance
The VPD results showed a better model performance for fire in the spatial dimension, as well as in the temporal dimension ( Tab. 1 and 3). Compared to the Nesterov index, F DI V P D uses additional climate input as relative humidity and precipitation. In the calculation of the Nesterov index precipitation is just used as a threshold. This leads to a better accounting of the very different climatic conditions among various biomes. Furthermore, the F DI V P D includes a direct representation of the 5 vegetation density. The significance of this has been recently shown by findings of Forkel et al. (2019a) who have emphasized the importance of past plant productivity and fuel production for burnt area. This is particularly important for differentiating between fires in biomes with similar PFT distribution. For example, the vegetation density is much larger in the Cerrado, even though the Caatinga and Cerrado have a similar modelled PFT composition, which provides more fuel and therefore leads to a higher fire danger.

10
While the seasonal and interannual variability in the Caatinga has improved largely using the F DI V P D (NMSE decreased by a factor of ca. 20), the improvement in the Cerrado was relatively small (NMSE decreased by ca. 10%). This is due to the fact, that the optimization tries to obtain a compromise between the different optimized cells. As the model performance was originally much better for the Cerrado, the largest improvement could be achieved for the Caatinga. We have also chosen a large amount of cells in the Caatinga, because the model performance was here particularly bad. This leads to a large improvement 15 in the time series of the Caatinga region, while the improvement for the Cerrado was less significant. With the Nesterov index fire was strongly underestimated in the Amazonia region, while the optimized VPD increases the modelled burnt area. The fire is only present at the edges of the Amazon (both in model and observation, see Fig. 4), where tree density is lower and deforestation takes place. In the closed continuous forest area towards the center of the Amazon almost no fire is observed and also not simulated. 20 Another result of the optimizing procedure, using F DI V P D , was the improvement of the PFT distribution and the aboveground biomass of trees especially in the fire-prone biomes Caatinga and Cerrado (Fig. 8). For example, the central Amazon, where fire is a scarce event, shows almost no changes compared to the non-optimized model version. Here, it is the improvement of the vegetation model itself, and not the fire module, which can help to improve the model performance of LPJmL4-SPITFIRE. Hence, it emphasizes that we need to include further parameters in the optimization, which impact directly the PFT 25 distribution, biomass and fire to obtain a significant improvement in the spatial and temporal distribution of both, vegetation and fire. However, this study focused solely on the parameters within the SPITFIRE module. Due to the focus on fire related parameter, the cost of the burnt area dataset decreases much more than the cost of the biomass dataset (Fig. A5). Hence we only get a substantial improvement in model performance in semi-arid, fire-prone biomes, where vegetation dynamics and fire are strongly coupled.

30
During the optimization-process most of the optimized parameters were well constrained, except for the mortality-related parameters for the TrBR PFT (Fig. 6). The TrBR PFT is dominant in the fire-prone regions, where the mortality-related parameters have a large impact on vegetation dynamics. Hence, they impact multiple LPJmL routines, which are responsible for the PFT distribution and carbon cycling. This leads in turn to a less certain parameter estimation. In order to better constrain these parameters also the optimization of vegetation model parameters would be necessary to decrease the uncertainties.
The fire danger index scaling factors (α N Ii and α V P Di ) convert the quantified fire risk (NI or VPD) into the actual fire danger (FDI). Both scaling factors thus set the magnitude of the fire danger for the different PFTs. Hence they impact directly the fire spread, burnt area and the number of fires as well as indirectly fire mortality. These very important parameters vary significantly for the different PFTs. TrH has the smallest scaling factor in case of both FDIs, which leads to a lower fire danger compared 5 to the other PFTs. This indicates a prior overestimation of the fire danger of grass in tropical South America, as grasslands are generally parametrized to have a low fire resistance and moisture content and can hence burn very easily. This overestimation, compared to tree PFTs has been decreased by the optimization. In case of the VPD also the TrBR is scaled by a much smaller factor than the TrBE, which leads to a lower fire danger index. This is due to the fact, that the TrBR is dominant in dry and fire-prone regions, which experience frequent fires. Here the burnt area was often overestimated by SPITFIRE (e.g. Caatinga 10 or eastern Cerrado) and is now decreased. On the other hand, a larger FDI for the TrBE allows more fire in wetter regions at the edge between the Cerrado and the Amazon rainforests, where TrBE is more dominant. The mortality risk of TrBE for VPD optim remains close to the prior value of 1, confirming previous assumptions about its fire sensitivity. Whereas the rCK for TrBR increased to 0.48, close to the upper boundary, meaning that a mortality risk of 50% when the full crown is scorched and a 7% mortality risk when 50% of the crown is scorched, which makes the TrBR less resistant against crown damage than 15 before. Due to this changes the overestimation of biomass in the original model for the Cerrado/Caatinga region decreased (see Fig. 7).

Limitations of the optimization process
Generally, optimizing a model against burned area is challenging because 1) of the skewed statistical distribution of burned area and 2) because temporal or spatial mismatches in simulated burning can cause large model-data errors. These issues can 20 be avoided with the choice of an appropriate cost function. For example, squared-error metrics tend to underestimate the variance of burned area in comparison to, e.g., the Kling-Gupta efficiency as it has been shown in the optimization of an empirical model for burned area (see Tab. A3 in Forkel et al. (2017)). Here, the optimum parameter set for the Nesterov index-based model resulted in almost no fires across South America. Thereby the optimization algorithm tries to decrease the model error by tending towards a conservative 'no fire strategy' for all biomes. This result nicely demonstrates the need to evaluate model 25 optimization results against spatially and temporally independent data and independent variables (Keenan et al., 2011).
The Nesterov index is not able to capture fire variability within the Caatinga as well as the Cerrado at the same time. This shows that the difference in the PFT distribution between these two biomes is not adequately modelled by LPJmL or just using PFT dependent scaling factors did not sufficiently improve the model performance when using the Nesterov index. On the other hand, using the VPD fire danger index reduced the model error for burned area in both biomes, by improving the modelled 30 performance for the Caatinga and maintaining the good performance of the Cerrado region. Since improved performance of the fire model mainly had minor effect on improving FPC of the tropical PFTs, the presented optimization scheme has to go along with process-based improvements in both, in the fire and in the vegetation modules of LPJmL.
Fire largely depends on the vegetation type and their associated flammability, fire tolerance and mortality. Hence an accurately modelled vegetation distribution is crucial for a good model performance in terms of burnt area and fire effects (Forkel et al., 2019a;Rogers et al., 2015). As shown in Fig. 8, A1 and A2, the modelled PFT coverage showed an equal distribution of tropical raingreen and evergreen PFTs throughout wide parts of central-northern South America. Evaluation data shows, however, an TrBE dominance in the wet rainforest regions and a TrBR dominance in the Cerrado and Caatinga. This emphasizes the potential to improve the fire modelling further, based on an improved PFT distribution. In the tropical rainforest the TrBR 5 proportion is overestimated, which leads to problems in the optimization procedure, since TrBR has very different effects on fire spread and is more fire-tolerant (different fuel characteristics and resulting fire intensity). This leads to a lower fire-related mortality, which fits better to the drier and fire prone savanna-like regions (e.g. Cerrado). The poorly modelled PFT distribution also is responsible for the overestimation of the burnt area in the Amazon region. Because of the too large fraction of TrBE in the Cerrado/Caatinga region the scaling factor for this PFT is relatively high. This leads in turn to an overestimation in the 10 Amazon region, where the fraction of the TrBE is larger.
Since the offset is very small, the years 2000-2003 (first three years of GLDAS 2.1, before the optimization period) are enough for the model to recover from the offset and the carbon pools to return to equilibrium. To exclude the possibility that long-term trends within GLDAS 2.0 changed the modelled vegetation state significantly, we tested our optimization also just based on GLDAS 2.0 data (until 2010) and on GLDAS 2.1 data (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017) only, using the same years for model spinup, optimization 15 and evaluation. Both versions yielded similar results compared to the optimization presented in this study (results not shown).
Due to the fact that evaluation data are only available for the last 10-20 years, we are constrained to optimize the model in this relatively short time period. In South America these years were subject to an unusual high amount of severe droughts and other extreme events (Panisset et al., 2017). As a result, an optimization in this period could lead to a worse model performance in a period with less pronounced droughts. This is due to the non-linear relationship between the drought signal in the input data 20 set and the resulting modelled biosphere behavior. Nonetheless, we were able to improve the interannual variability and hence, the model performance to a great extent for the Caatinga and slightly for the Cerrado and Amazon regions ( Fig. 5 and A3). The Cerrado already had a very good modelling performance before the optimization process, which now only slightly improved.
The performance of the interannual and seasonal variability of burnt area for total South America improved substantially (Fig.   A3). The optimized SPITFIRE is now better able to simulate accurately the climate dependent seasonal and interannual vari-25 ability as well as the spatial extent of fire on natural land throughout the fire-prone woodlands of South America.
Systematic optimizations within a model-data integration setup of fire models which are embedded in a DGVM are still very rare. Previously, Rabin et al. (2018)  Fire models embedded in DGVMs should build on a FDI which is complex enough to account for various fire dynamics, while it's parameterization should be simple enough to be accurately applied on a global scale. While the VPD is more complex and takes into account more climatic input as the Nesterov index, it is still relatively easy to implement in a global fire model.
There are various other fire danger indices used for modelling purposes, as well as real fire danger assessment and fire forecast 20 purposes. For example, fire-prone countries have developed their own fire danger indices (e.g. Canada, Australia), which are suited to the unique local fire regimes and vegetation dynamics. In a global modelling approach, however, we need to find one fire danger index, which suits best for all regions of the world and has a relatively easy implementation to decrease computational cost and the number of input data sets (which might be unavailable or uncertain).
Currently, SPITFIRE does not account for fire in managed land like cropland or managed grassland. We accounted for this by 25 excluding cropland fires from the evaluation burnt area data set. We do, however, not account for the proportion of grassland, which is used for e.g. cattle ranching. Since in SPITFIRE fire is not enabled on pasture, our results show a slightly smaller burnt area throughout South America than could be expected with managed land included and hence also compared to the GFED4 evaluation data set. This effect is however small, because pasture lands cover a substantial fraction only in very few grid cells (e.g. southern Cerrado; Parente et al., 2017). Fire on managed land is generally difficult to predict in a DGVM because the 30 reason and timing of using fire depends less on climatic factors but mostly on social and political decisions which can vary between countries, regions and localities. We expect further improvement of model performance especially in regions of large land-use areas with fires on pastures included (e.g. Rabin et al., 2018;Pfeiffer et al., 2013).
We significantly improved the fire representation within LPJmL4-SPITFIRE, applied for South America, by implementing a new fire danger index and applying a model-data integration setup to optimize fire-related parameters. We improved the seasonal and interannual variability, as well as the spatial pattern of burnt area in South America. In addition, modelling of related vegetation variables, e.g., the biomass and the PFT distribution in the fire-prone Cerrado and Caatinga biomes have also 5 been improved.
Optimizing fire parameters has its limits due to error propagation of the PFT distribution and hence their fire traits influencing simulated fire spread and behavior. Furthermore, it remains a challenge to find a fire danger index that is physically interpretable and can be applied globally. In this study, the parameter-optimization by using F DI N I led to a large underestimation of fire and a generally worse model performance, when focusing on the Cerrado and Caatinga biome. However, implementing the more  that are related to model-data bias, variance ratio and correlation. The cost for burnt area for NIoptim decreased by ca. 81%, whereas the cost of the biomass only decreases by ca. 6% (a and b). For VPDoptim the cost decreased by ca. 48% for burnt area and ca. 19% for the biomass (c and d). Hence the impact of the optimization process on burnt area is much larger due to the focus on fire parameters.  Competing interests. The authors declare that they have no conflict of interest.