A hydroclimatic model for the distribution of fire on Earth

The distribution of fire on Earth has been monitored from space for several decades, yet the geography of global fire regimes has proven difficult to reproduce from interactions of climate, vegetation, terrain and land use by empirical and process-based fire models. Here, we propose a simple, yet robust, model for global fire potential based on fundamental biophysical constraints controlling fire activity in all biomes. In our ‘top-down’ approach we ignored the dynamics of 10 individual fires and focus on capturing hydroclimatic constraints on the production and (seasonal) desiccation of fuels to predict the potential mean annual fractional burned area, here estimated by the 0.99 percentile of the observed mean annual fractional burned area (F . ). We show that 80% of the global variation in F . can be explained from a combination of mean annual precipitation and potential evapotranspiration. The proposed hydroclimatic model reproduced observed fire activity levels equally well across all biomes and provided the first objective underpinning for the dichotomy of global fire regimes in 15 two domains characterised by either fuel production limitations on fire or fuel dryness limitations on fire. A sharp transition between the two climate-fire domains was found to occur at a mean annual aridity index of 1.9 (1.94±0.02). Our model provides a simple but comprehensive basis for predicting fire potential under current and future climates, as well as an overarching framework for estimating effects of human activity via ignition regimes and manipulation of vegetation. 20

complexity to global fire patterns by, amongst others, changing vegetation community composition, structure and flammability, active use of fire to clear land or reduce fire hazard, and fire suppression Archibald et al., 2012;Parisien et al., 2016;Bowman et al., 2011). The partial success of current models in predicting global fire patterns suggests that their fire modules fail to capture some aspects of the biophysics that control fire activity across different environments.
Incomplete understanding of biophysical drivers and constraints that underlie current global fire patterns creates uncertainty 35 in model predictions of how fire regimes, fire-prone ecosystems and related biogeochemical cycles may respond to rising atmospheric [CO2] and climate change (Harris et al., 2016).
Here, we go back to the fundamental biophysics of fire to formulate a simple, yet robust, model for the prediction of potential mean annual levels of fire activity across all global biomes. For this study we define fire activity in terms of the mean annual 40 fractional burned area, , of the landscape or land area unit, ranging from 0 in long unburnt landscapes to 1 in landscapes that burn completely on an annual basis (Giglio et al., 2013). By predicting we ignore the dynamics of individual fires, but this is appropriate in our view given the fact that most fires are orders of magnitude smaller (e.g., Archibald et al., 2013;Malamud et al., 2005) than the typical grid size (e.g. 0.5° x 0.5°) of LSMs and DGVMs.

45
Building on Bradstock's (2010) 'four-switch' concept we assume that four fundamental conditions need to be met for a landscape fire to occur: i) there must be enough plant biomass (i.e. fuel) to carry a fire, ii) the extant fuel must be dry enough to be ignitable, iii) weather conditions need to be favourable (i.e. hot, dry and windy) for a fire to spread, and iv) there must be an ignition. Bradstock (2010) conceptualized these conditions as four 'switches' in a series circuit that need to be 'on' for a fire to occur. While any fire will require alignment of all four switches at some point in time, in the context of modelling 50 global fire patterns we emphasize that the four switches operate at disparate rates, with an associated hierarchy of conditional constraints on fire: production of plant biomass and build-up of fuel loads occurs over months to years, fuels dry out over weeks to months, and fire weather varies over time scales of hours to days, while ignitions are instantaneous events. Therefore, we hypothesize that the mean annual fractional burned area ( ) can be predicted from long-term fuel production and fuel drying rates (i.e. switches 1-2), while information on fire weather and ignitions (i.e. switches 3-4) is only required when 55 modelling the specific attributes of individual fires (e.g. size and burn pattern). The hierarchical organization of the four switches further implies that the fractional burned area predicted from long-term fuel production and fuel drying rates alone represents an upper limit of fire activity that is only reached when fire weather and ignition limitations are minimal.
To test the hypothesis that the upper quantiles of can be predicted from long-term fuel production and drying rates, we 60 analysed global burned area data together with indices of fuel productivity and fuel dryness, both calculated from the climatic water balance (Stephenson, 1998). Building on the methods developed by a previous study of Australian fire regimes , we propose a new global model that predicts the upper limit (i.e. the 0.99 quantile) of the mean annual fractional burned area, . , from two basic hydroclimatic variables: mean annual precipitation, , and potential evapotranspiration, .
The model also predicts the relative importance of either fuel productivity or fuel dryness constraints on . , which we 65 hypothesized to vary with the global distribution of land cover types and corresponding fuel types. Consistent with the variable constraints hypothesis , in dryland environments with grassy fuels we expected . to be primarily limited by fuel productivity constraints, while in more mesic environments dominated by woody vegetation and litter fuels . was expected to be limited primarily by fuel dryness constraints. Finally, we explored whether the global distribution of contemporary fire regime classes or 'pyromes' as classified by Archibald et al. (2013) on the basis of four fire 70 regime metrics (i.e. fire return interval, maximum fire intensity, length of fire season, maximum fire area.) was associated with the relative importance of fuel productivity or fuel dryness constraints on . .

Modelling approach 75
In this study our aim was to predict global patterns of potential mean annual fractional burned area ( . ) as a function of basic biophysical constraints on the production and dryness of fuel material. Following Boer et al. (2016), we assumed that both the production and drying of fuel material are essentially functions of the local water and energy budgets available for the production and desiccation of plant biomass. The long-term climatic water balance, calculated from mean annual and , captures these interactions of biologically available water and energy (Stephenson, 1998). When calculated over long time 80 scales (>>years) and broad spatial scales (>>km 2 ), changes in the soil water store and lateral water inputs can be assumed to be negligible so that the climatic water balance is reduced to: 0, where is runoff/drainage losses and is actual evapotranspiration. is a reliable predictor of continental patterns of annual primary productivity at annual timescales (Rosenzweig, 1968;Yang et al., 2013) and therefore a reasonable proxy for fuel production rates (Meentemeyer et al., 1982;Matthews, 1997). Similarly, the potential for drying of fuel material can be assumed to be proportional to the atmospheric 85 moisture demand (∝ ) that cannot be met by available water (∝ ), which is the climatic water deficit as defined by Stephenson (1998), . For long timescales and large land areas, can be estimated from and using the semiempirical Budyko curve (Budyko, 1958;Zhang et al., 2010;Williams et al., 2012): where is the global maximum of . , here set at 1. The shape of . and . is set by two parameters (Yin et al., 2003): (or ) is the value at which . increases most strongly with (or ), and (or ) the value at which . becomes irresponsive to further increase of (or ). Since the two predictor variables, and , are not independent, the fitted , , . response surface was mapped to (orthogonal) axes of and before interpreting the shape of the response surface 100 in terms of biophysical constrains on global fire activity.

Burned area
Global annual burned area data at 0.25° x 0.25° spatial resolution for the period July 1995 -June 2016 were obtained from the 105 GFED4 database (Giglio et al., 2013). The mean annual fractional burned area, , was calculated by summing the burned areas within each grid cell for the entire observation period and dividing by the area of the grid cell and the duration of the observation period.

Landcover and fuel types
This study focused on areas of (semi-)natural vegetation. The corresponding vegetation mask was constructed from the Land Cover Type product (MCD12Q1) of the Moderate Resolution Imaging Spectroradiometer (MODIS) by first resampling the 120 data layer to the 0.25° grid of the GFED4 data base, using nearest neighbour interpolation, and reclassifying grid cells of excluded land cover types (i.e. water, cropland, urban and built-up, cropland/natural vegetation mosaics, snow and ice, barren or sparsely vegetated) to missing values.

Model fitting and validation
Data analyses focused on the modelling of global . as a function of mean annual and (Eq. 2-4) and on the interpretation of the model in terms of its consistency with current understanding of global fire patterns. Our complete data set consisted of estimates of mean annual , , , and for 193,476 grid cells covering the selected land cover types of the global land area, except Antarctica, at 0.25° x 0.25° spatial resolution. A randomly selected sample of 50% of the grid cells was used for model 130 fitting, while the other 50% of the data was set aside for model validation. We used R (R Core Team, 2018) for all data analyses, in particular the 'raster' (Hijmans, 2016) and 'quantreg' packages (Koenker, 2017).
Non-linear quantile regression was used to fit Eq (4) to the 0.99 quantile of as a function of mean annual and . To minimize bias in the model towards the most common global climates (e.g. desert or boreal climates), we used a simple 135 bootstrap procedure of two steps: i) the global , space was divided into 100 mm by 100 mm bins and all climate bins with a minimum of 100 grid cells identified (n=145), ii) a random sample (with replacement) of 100 grid cells was drawn from these bins, and iii) Eq (4) was fitted to the sample data. This procedure was run 1000 times to generate 1000 response surfaces from which a mean . response surface was calculated (Geyer, 2011). The means and confidence intervals of the fitted model coefficients are listed in Table 1. 140 Model predictions were validated against the 50% of the data that was not used for model fitting. To do so the validation data 100 mm , bins. The relationship between observed and predicted . was evaluated by linear regression analysis (Supplementary Material, S1). Deviations between observed and predicted . were quantified using the mean difference (MD) and root mean squared difference (RMSD): where , are predicted and observed values of . and n=191. To measure agreement between the spatial patterns of observed and predicted . , we used the normalised mean error (NME) and normalised mean squared error (NMSE), as proposed by Kelley et al. (2013): where is the predicted value of . at grid cell i, the corresponding observed value, and the mean of all observed values. By normalising by the spatial variability of the observations, NME and NMSE provide a measure of the spatial error of the model, with NME or NMSE close to 0 indicating perfect agreement between observed and predicted patterns, and both metrics approaching unity when agreement is similar to that of a model that predicts a spatially uniform value equal to the mean of all observations (Kelley et al., 2013). 160

Identification of climate-fire domains
Following Boer et al. (2016) two climate-fire domains were distinguished on the , , . response surface depending on whether the direction of the . gradient was more parallel to the local gradient or gradient, indicating predominance of fuel productivity or fuel dryness limitation on fractional burned area, respectively. The boundary between the two climate-fire 165 domains was identified analytically using a gradient analysis of the . response surface relative to gradients of and in , space (Supplementary Material, S3).
As in Boer et al. (2016) we refer to these two domains as productivity-limited (PL) fire and dryness-limited (DL) fire domains and analysed whether the affinity to either domain was related to the vegetation type being dominated by grasses/herbaceous or woody plants. To this end, areas of homogeneous land cover type were identified on the GlobeLand30 map (Chen et al., 170 2015). The GlobeLand30 product  is based on 30 m resolution Landsat imagery and classifies land cover types according to the dominant plant life form (e.g. forest, shrubland and grassland), which can be more readily related to distinct fuel types than biome classifications that often include classes of mixed life forms (e.g. woody savanna in MCD12Q1). pyrome map (i.e. Fig. 2 in Archibald et al., 2013) to the , , . response surface via their corresponding mean annual and , ii) dividing the global , space in 100 mm x 100 mm bins and identifying the pyrome class with the highest relative frequency in all , bins with a minimum of 50 data points (N=185), iii) a qualitative evaluation of the consistency of the predicted . and relative importance of fuel production and fuel dryness constraints on fire with the defining characteristics of the five pyrome classes . 185

Burned area
According to the GFED4 data base (Giglio et al., 2013), global burned areas amounted to 2.3% of the terrestrial land area per The 0.99 quantile of the mean annual fractional burned area, . , was found to be a highly predictable function of the climatic water balance terms and , and therefore of mean annual precipitation ( ) and potential evapotranspiration ( ) (Fig. 1a-c).
Linear regression analysis of observed versus predicted . for validation sites showed that the hydroclimatic model (Eq. 4,195  Africa and Australia, where tropical wet-dry climates combine high levels of fuel production during the wet season with intense drought during the dry season, producing values as high as 0.8-1.0 ( Fig. 2a-b). In woody ecosystems outside of the tropical savannas and semiarid grasslands . seldom exceeded 0.5, which is consistent with the fact that most predominantly woody vegetation communities cannot survive such high levels of fire activity over long periods (Bond and Keeley, 2005

Climate-fire domains
The fitted , , . response surface consists of two distinct domains (Fig. 1d) characterized by different climate constraints on fire . The first domain (green zone in Fig. 1d) is characterized by . increasing with mean annual actual 210 evapotranspiration ( ) but not with variation in climatic water deficit ( ), consistent with potential fire activity levels being primarily limited by fuel production (hereafter: PL-type fire). In the second domain (orange zone in Fig. 1d), . increases strongly with but varies little with increasing , consistent with potential fire activity levels being primarily limited by the capacity of the atmosphere to dry-out fuel material to ignitable levels (DL-type fire).
A visual interpretation of the fitted , , . response surface suggests that the domain shift from PL-type fire to DL-type 215 fire occurs at some threshold aridity index (i.e. the ratio of mean annual potential evapotranspiration and precipitation: ). (1-4)  Using gridded mean annual climate data, the domains of PL-and DL-type fire were mapped to geographical space (Fig. 2c). 220

The equation for the boundary between the two domains can be derived analytically from Equations
The global pattern of PL-and DL-type fire is similar to the global distribution of climate aridity and biome types, as expected given that the boundary between PL-and DL-type fire corresponds to an aridity index of ~1.9. Dryland environments at mid latitudes on all continents were classified in the domain of PL-type fire, while wet and cold environments at all latitudes were classified in the domain of DL-type fire. The domain classification indicates whether the primary limiting factor on mean annual fire activity levels was fuel production or fuel dryness, which can be expected to correlate strongly with the dominant 225 vegetation lifeform and fuel type (litter versus grass). Using the GLC30 global lifeform mapping  we found that areas of forest, wetland and tundra were predominantly classified in the domain of DL-type fire (i.e. they are typically too wet to burn for much of the time), whereas grasslands, shrublands and barren lands were predominantly classified in the domain of PL-type fire (i.e. fuels are typically sparse and discontinuous for much of the time) (Fig. 3). The shrubland class is an interesting case: these ecosystems occur most frequently in the domain of PL-type fire even though they are dominated by 230 woody vegetation and woody plant litter forms an important fuel component. However, classification as PL-type fire makes sense because (semiarid and arid) shrublands can often only support large fires when herbaceous vegetation fills in the gaps and connects the fuel array after above-average rainfall (e.g., O'Donnell et al., 2011;Prior et al., 2017).
As expected, the five pyromes distinguished by Archibald et al. (2013) were found to occupy distinct domains in , climate space ( Fig. 4): i) the 'RIL' pyrome, characterized by rare, intense, large fires was restricted to low environments in the 235 climate domain of PL-type fire (green grid cells, Fig. 4) with a median predicted . of 0.05, consistent with its geographical distribution in the more arid zones of boreal forests and temperate coniferous forests, plus areas of Mediterranean vegetation and xeric vegetation; ii) the 'FIL' pyrome, characterized by frequent, intense, large fires (yellow grid cells, Fig. 4 in the climate domain of PL-type fire across a broad range of mean annual combined with high and with a median predicted . of 0.55, which is typical for the tropical savannas in Australia and Africa where this pyrome prevails; iii) the 240 'FCS' pyrome, characterized by frequent, cool (low-intensity), small fires (orange grid cells, Fig. 4), was found across both climate domains in the region that combines high , high , and very high predicted . (median=0.63), corresponding mainly to tropical grasslands and shrublands as well as tropical dry broadleaf forests; iv) the 'RCS' pyrome, characterized by rare, cool, small fires (pink grid cells, Fig. 4)

Discussion 255
This study has shown that climate sets strong and highly predictable constraints on the global distribution of fire on Earth. In particular, climate constrains the amounts and timing of plant available water and energy and thereby determines the probability that the two most basic conditions for fire are met, namely the production and desiccation of plant biomass and derived fuels. We showed that the strength of those two fundamental climate constraints on fire, and global variation therein, are captured well by the mean annual climatic water balance, which provides a simple yet biophysically sound basis for the 260 prediction of potential fractional burned area from just two readily available climate variables: mean annual precipitation and potential evapotranspiration. The proposed hydroclimatic model was validated against independent burned area data and found to explain 80% of the global variation in potential fractional burned area ( . ) with a slight tendency to overpredict . , but an NME of 0.45 indicating good agreement between predicted and observed spatial patterns of . . A direct comparison of model performance with existing global fire models is difficult since most existing models predict rather than . and 265 systematic evaluation of their performance is ongoing as part of the fire modelling intercomparison project (FIREMIP) (Hantson et al., 2016;Rabin et al., 2017;Teckentrup et al., 2019;Forkel et al., 2019) for their model predictions of mean annual burned area against the GFED4s data set (Randerson et al., 2012;Giglio et al., 270 2013), indicating a lack of agreement between the spatial pattern of predicted and observed burned area.
Overprediction of . by our model was somewhat more pronounced in environments of the DL-type fire than in environments of PL-type fire (see Supplementary Material, S1), but model residuals had similar distributions for a wide range of MODIS landcover types, supporting the notion that the model captures the main climatic constraints on global fire activity levels. Existing fire models tend to struggle predicting the fire regimes of temperate and boreal forest regions, that burn much 275 more infrequently than tropical savannas , but we did not observe that model fit suffered in any particular biome because we employed a mechanistic, hydroclimatic approach that captured the fundamental biophysics underlying global fire regimes.
The bias in performance of existing global fire models is likely due to a limited capacity to simulate fuel drying dynamics in forest and woodland environments. Whereas the seasonally wet-dry climate of tropical savannas makes annual production and 280 subsequent desiccation of fuels highly predictable, fire activity in forests and woodlands is primarily constrained by the moisture content of the (inherently abundant) fuels, which varies at (sub)daily to monthly timescales (Nolan et al., 2016;Caccamo et al., 2012;Boer et al., 2017). The fuel moisture models within several of the existing global fire models (Rabin et al., 2017) are driven by some implementation of the Nesterov index: ∑ , , where is the number of days since the last rainfall exceeding 3 mm d -1 , is the daily 3 pm air temperature (°C) and , the daily 3pm dew point 285 temperature (°C). As shown by Schunk et al.(2017) the Nesterov Index was a poor predictor of the moisture content of 1-hour and 10-hour fuels from four main forest species in Germany. With fuel moisture variation modelled by the Nesterov index global fire models lack the critical capacity to predict when, or how frequently, forests and woodlands switch from a moist/nonflammable state to a dry/flammable state and are therefore unlikely to reproduce observed spatiotemporal variation in burned area in those biomes. In contrast, our hydroclimatic model does not predict fuel moisture content; instead, the mean annual 290 climatic water deficit ( ) is used as an estimate of the probability of extant fuels drying out to ignitable levels during some fraction of the year. With being a measure of the atmospheric moisture demand that cannot be met by the soil water store (Stephenson, 1998), captures the basic biophysics involved in the desiccation of fuels and accounts for the fact that sparse fuels (low ) require less energy ( ) to dry out than dense/heavy fuels.

295
Further evidence for being a reasonable indicator of the mean annual probability of forests and woodlands being in a dry fuel state can be derived from previous studies on climate-fire relations in the southwest Unites States (Williams et al., 2015;Abatzoglou et al., 2017) and studies on dead fuel moisture (Resco de Dios et al. 2015) and fire activity in temperate forests of SE Australia (Nolan et al., 2016) and Portugal (Boer et al., 2017) that showed that cumulative burned area in these different forest regions responds strongly non-linearly to predicted fine dead fuel moisture content dropping below thresholds 300 identified at 14-18% and 10-12%, respectively. Using the Resco de Dios et al. (2015) fuel moisture model with gridded global vapour pressure deficit data to predict daily fine dead fuel moisture content for global forests and woodlands, we found that https://doi.org/10.5194/bg-2019-441 Preprint. Discussion started: 25 November 2019 c Author(s) 2019. CC BY 4.0 License. mean annual is strongly and linearly related (adj. R 2 : 0.76, p<2e-16) to the mean annual frequency of predicted daily fine dead fuel moisture content dropping below 10% (Supplementary Material, S4).
The hydroclimatic model allowed us to objectively distinguish climatic domains for a predominance of either fuel productivity 305 (PL-type) or fuel dryness (DL-type) constraints on mean annual fractional burned area, and showed their geographical distribution to be consistent with global patterns of herbaceous vs. woody vegetation types and corresponding fuel types, in accordance with the variable constraints hypothesis . Further, we demonstrated (see Supplementary Material, S3) that the boundary between the domains of PL-and DL-type fire is well approximated by an aridity index of 1.9 (1.94±0.02), providing the first objective identification of where in climate space (Fig. 1d), and in 310 geographical space (Fig. 2c), fire regimes switch from fuel load limitations to fuel moisture limitations. We showed that the hydroclimatic model and associated classification in two contrasting climate-fire domains is consistent with the hydroclimatic distribution of pyromes , indicating that key aspects of global fire regimes vary in a predictable way with global gradients in mean annual precipitation and potential evapotranspiration.
By focusing exclusively on the roles of fuel production and fuel dryness constraints, our hydroclimatic model was designed to 315 predict the potential or maximum mean annual fractional burned area, which provides a useful reference for bottom-up processbased modelling approaches used in many DGVMs (Hantson et al., 2016). Further work should be extended to predicting annual, rather than potential mean annual, fractional burned area by modelling the deviations between predicted . and observed annual fractional burned area as a function of fire weather and ignition constraints on fire activity, thus completing the formalisation and implementation of the four-switch concept (Bradstock, 2010). Other future work could also examine the 320 drivers, such as fire management and other human activities, of geographical variation between and . .

Conclusion
At long time scales the global distribution of fire on Earth is highly predictable from fundamental biophysical constraints on the production and seasonal desiccation of plant biomass (i.e. fuel), which in turn are proportional to mean annual precipitation 325 and potential evapotranspiration. The sharp transition of global fire regimes from domains of fuel production limitations on fire (PL-type) to fuel dryness limitations on fire (DL-type) can be identified from the mean annual aridity index being above or below a threshold value of 1.9. Our model provides a simple but comprehensive basis for predicting fire potential under current and future climates, as well as an overarching framework for estimating effects of human activity via ignition regimes and manipulation of vegetation. In these respects, it offers a significant advance on existing global fire models and therefore a 330 basis for improving predictions from coupled global vegetation models.