Modelling Global Burned Area and Fire Regime Printer-friendly Version Interactive Discussion Geoscientific Model Development Modelling Fires in the Terrestrial Carbon Balance by Incorporating Spitfire into the Global Vegetation Model Orchidee – Part 1: Simulating Historical Global Burned Area and Fi

Discussions This discussion paper is/has been under review for the journal Geoscientific Model Development (GMD). Please refer to the corresponding final paper in GMD if available. Abstract Fire is an important global ecological process that determines the distribution of biomes, with consequences for carbon, water, and energy budgets. The modelling of fire is critical for understanding its role in both historical and future changes in terrestrial ecosystems and the climate system. This study incorporates the process-based 5 prognostic fire module SPITFIRE into the global vegetation model ORCHIDEE, which was then used to simulate the historical burned area and the fire regime for the 20th century. For 2001–2006, the simulated global spatial extent of fire occurrence agrees well with that given by the satellite-derived burned area datasets (L3JRC, GLOB-CARBON, GFED3.1) and captures 78–92 % of global total burned area depending 10 on which dataset is used for comparison. The simulated global annual burned area is 329 Mha yr −1 , which falls within the range of 287–384 Mha yr −1 given by the three global observation datasets and is close to the 344 Mha yr −1 given by GFED3.1 data when crop fires are excluded. The simulated long-term trends of burned area agree best with the observation data in regions where fire is mainly driven by the climate 15 variation, such as boreal Russia (1920–2009), and the US state of Alaska and Canada (1950–2009). At the global scale, the simulated decadal fire trend over the 20th century is in moderate agreement with the historical reconstruction, possibly because of the uncertainties of past estimates, and because land-use change fires and fire suppression are not explicitly included in the model. Over the globe, the size of large fires (the 95th 20 quantile fire size) is systematically underestimated by the model compared with the fire patch data as reconstructed from MODIS 500 m burned area data. Two case studies of fire size distribution in boreal North America and southern Africa indicate that both the number and the size of big fires are underestimated, which could be related with too low fire spread rate (in the case of static vegetation) and fire duration time. Future 25 efforts should be directed towards building consistent spatial observation datasets for key parameters of the model in order to constrain the model error at each key step of the fire modelling.

prognostic fire module SPITFIRE into the global vegetation model ORCHIDEE, which was then used to simulate the historical burned area and the fire regime for the 20th century. For 2001, the simulated global spatial extent of fire occurrence agrees well with that given by the satellite-derived burned area datasets (L3JRC, GLOB-CARBON, GFED3.1) and captures 78-92 % of global total burned area depending 10 on which dataset is used for comparison. The simulated global annual burned area is 329 Mha yr −1 , which falls within the range of 287-384 Mha yr −1 given by the three global observation datasets and is close to the 344 Mha yr −1 given by GFED3.1 data when crop fires are excluded. The simulated long-term trends of burned area agree best with the observation data in regions where fire is mainly driven by the climate 15 variation, such as boreal Russia (1920Russia ( -2009, and the US state of Alaska andCanada (1950-2009). At the global scale, the simulated decadal fire trend over the 20th century is in moderate agreement with the historical reconstruction, possibly because of the uncertainties of past estimates, and because land-use change fires and fire suppression are not explicitly included in the model. Over the globe, the size of large fires (the 95th 20 quantile fire size) is systematically underestimated by the model compared with the fire patch data as reconstructed from MODIS 500 m burned area data. Two case studies of fire size distribution in boreal North America and southern Africa indicate that both the number and the size of big fires are underestimated, which could be related with too low fire spread rate (in the case of static vegetation) and fire duration time. Future

Introduction
Fire is an important process in the Earth system, that existed long before the largescale appropriation of natural ecosystems by humans (Bowman et al., 2009;Daniau et al., 2013). Fires have multiple biophysical and ecological consequences, and they are also an important source of atmospheric trace gases and aerosol particles (Lang-5 mann et al., 2009;van der Werf et al., 2010). By damaging some plant types and concurrently promoting others, fires play an important role in shaping vegetation structure and function (Bond et al., 2005;Pausas and Keeley, 2009). Fire changes the surface albedo, aerodynamic roughness, and the sensible and latent heat fluxes; fireinduced ecosystem change may thus influence the surface energy budget, and in turn Introduction

Model description
The processes and equations of the fire model SPITFIRE, as described by Thonicke et al. (2010), were implemented in the vegetation model ORCHIDEE (Krinner et al., 2005). The SPITFIRE model operates on a daily time step, consistent with the STO-5 MATE sub-module in ORCHIDEE, which simulates vegetation carbon cycle processes (photosynthates allocation, litterfall, litter and soil carbon decomposition). The major processes within SPITFIRE are briefly described below with applicable minor modifications (see Thonicke et al., 2010 for more detailed information).
1. Potential ignition includes ignitions by both lightning and human activity. Re-10 motely sensed lightning flash counts (Cecil et al., 2012) were obtained from the High Resolution Monthly Climatology of lightning flashes by the Lightning Imaging Sensor-Optical Transient Detector (LIS/OTD) (http://gcmd.nasa.gov/records/ GCMD_lohrmc.html). The LIS/OTD dataset provides annual mean flash rate over the period of 1995-2000 on a 0.5 • grid with a monthly time step. This annual data 15 was used each year throughout the simulation. Following Prentice et al. (2011), the proportion of lightning flashes that reach the ground with sufficient energy to ignite a fire is taken as 0.03. This value differs from that in the original SPITFIRE model as implemented in LPJ-DGVM (Thonicke et al., 2010); there a cloud-toground (CG) flashing ratio of 0.2 was used, followed by a further ignition efficiency 20 of 0.04.
To estimate potential ignitions by humans, the original Eqs. (3) and (4) in Thonicke et al. (2010) were modified for the purpose of unit adjustment, as below: √ PD × a(ND)/10 000 (1) 25 where I H is the daily ignition number (1 day −1 km −2 ), PD is population density (individuals km −2 ). The parameter a(ND) ( propensity of people to produce ignition events; it is a spatially explicit parameter in the Thonicke et al. (2010) LPJ-DGVM implementation, but following Venevsky et al. (2002), a(ND) is assumed to be a global constant of 0.22 in our simulation.
The spatial distribution of a(ND) used by Thonicke et al. (2010) is given in Supplement Fig. S1. The a(ND) values for many regions are within 50 % of the uni-5 form value used here, although northern Australia is an exception. Test simulations showed that using a uniform a(ND) value rather than the original spatial dataset has little influence on the simulated regional burned area (see Supplement Fig. S2).
2. Fire numbers are derived by scaling the potential ignitions (which include both 10 human and lightning ignitions) using the fire danger index (FDI), which is derived by comparing the simulated daily fuel moisture to a Plant Functional Type (PFT) dependent moisture of extinction. All fires with a fireline intensity lower than 50 kW m −1 are assumed unable to propagate and therefore are suppressed as stated by Thonicke et al. (2010). 15 3. Mean fire size is calculated by assuming an elliptical shape of fire, with the major axis length being the product of fire spread rate and fire duration time. Fire spread rate is obtained using the Rothermel equation (Rothermel, 1972;Wilson, 1982 whose designation in terms of hours describes the order of magnitude of time required for fuel to lose (or gain) 63 % of the difference between its current moisture content and the equilibrium moisture content under defined atmospheric conditions (Thonicke et al., 2010). Fuel combustion completeness is simulated as a function of daily fuel moisture, with a smaller fraction of fuel being consumed at 10 higher fuel moisture. Crown fuel consumption is related to the fraction of crown that is scorched by the fire flame. All PFT-dependent parameters follow Table 1 in Thonicke et al. (2010).

Further modifications made to the SPITFIRE equations
Fires in dry climate regions are limited by the availability of fuel on the ground 15 (Krawchuk and Moritz, 2010;Prentice et al., 2011;van der Werf et al., 2008b). This constraint is implicitly included in the SPITFIRE equations because fires are limited by a threshold of fireline intensity of 50 kW m −1 . However in model testing, we found that this threshold is too low to limit fires in low productivity regions (with modelled annual Net Primary Productivity or NPP of 0-400 g C m −2 yr −1 , corresponding to an 20 annual precipitation of 0-400 mm), with too much burned area being simulated for arid and semi-arid regions (see Supplement Fig. S3). Following Arora and Boer (2005), we therefore introduced a new factor that limits the ignition efficiency depending on the availability of ground fuel. Ignition efficiency varies linearly between zero when ground fuel is lower than 200 g C m −2 , to unity when ground fuel is above 1000 g C m −2 . Here, Introduction The equations for surface fuel combustion completeness given by Thonicke et al. (2010) follow Peterson and Ryan (1986), which allow combustion completeness to decrease with increasing fuel wetness and level out when fuel wetness drops below a threshold (Fig. 1). During model testing, we found that because fuel wetness frequently approaches zero, the simulated fuel combustion completeness is much higher 5 than the field experiment values reported by van Leeuwen et al. (2014), who compiled the combustion completeness measurements for fires in different biomes across the globe. We therefore modified the maximum combustion completeness for fuel classes of 100 h and 1000 h to be the same as the mean combustion completeness compiled by van Leeuwen et al. (2014). This biome-dependent maximum combustion completeness is 0.48 for tropical broadleaf evergreen and seasonal forests, 0.45 for temperate forest, 0.41 for boreal forest, 0.85 for grassland, and 0.35 for agricultural land.

Input dataset and the simulation protocol
The six-hourly climate fields used to drive the model were from the CRU-NCEP dataset (http://dods.extra.cea.fr/store/p529viov/cruncep/V4_1901_2012/readme.htm). 15 Population density for the 20th century was retrieved from the History Database of the Global Environment (HYDE) as compiled by the Netherlands Environmental Assessment Agency (http://themasites.pbl.nl/tridion/en/themasites/hyde/download/ index-2.html). From 1850 till 2005, the HYDE gridded data are available for the beginning of each decade and for 2005. Annual data were linearly interpolated within each 20 decade, and further re-gridded to 0.5 • resolution. For 2006-2009, population density was set constant at the 2005 value.
A three-step simulation protocol was used. For the first two steps, the atmospheric CO 2 was fixed at the pre-industrial level (285 ppm) and the climate forcing data of 1901-1930 were repeated in loop. The first step was a without-fire spin-up from bare 25 ground lasting for 200 years (including a 3000-year run of soil-only processes to speed up the equilibrium of mineral soil carbon). The second step was a fire-disturbed spin-up lasting for 150 years, with fire being switched on to account for fire disturbances in the 2385 Introduction Agricultural fires are not simulated in ORCHIDEE, even though the model has two PFTs that approximately represent C3 and C4 crops (but without realistic speciesspecific phenology). Magi et al. (2012) show that agricultural fire seasons differ significantly from those of natural fires, warranting a special treatment of agricultural fires 20 in global fire modelling (Li et al., 2012). Agricultural fires make up a rather small proportion of the total global amount in terms of both burned area (4.3 % according to the GFED3.1 dataset) and carbon emissions (less than 2 % according to GFED3.1). Further, deforestation fires are not explicitly simulated. Evidences show that deforestation fires occur during the "time window" when climate is dry enough to allow complete 25 burning (van der Werf et al., 2008a). The simulated fire danger index should partly capture this fire climate window, and that the model is able to account for some "deforestation" fire activity in tropical and subtropical forests, but not all, because our land cover map is static.

Datasets used to evaluate model performance
Several datasets were prepared and used to compare simulated burned area and fire regime with global Earth observations and government fire agency survey data.

Spatially gridded burned area data
Satellite-derived burned area data 5 The GFED3.1 fire carbon emission dataset provides monthly burned area data at 0.5 • resolution for 1997-2009 with global coverage . The GFED3.1 data were mainly generated using MODIS imagery with additional images from TRMM IRS and ERS ATSR. The fire carbon emissions data were also provided which are model simulation results by applying the CASA model .
L3JRC dataset provides daily global burned area data at 1 km resolution for April 2000 to March 2007; these data were generated from the 1 km SPOT VEGE-TATION satellite imagery (Tansey et al., 2008). This dataset was assembled at 0.5 • resolution at monthly time step for use in the present study.
GLOBCARBON burned area dataset was produced from a combination of SPOT 15 VEGETATION and ERS2-ATSR2/ENVISAT AATSR data as one of the four land products of the ESA GLOBCARBON initiative (Plummer et al., 2007). Global burned area data were provided at monthly resolution with four different spatial resolutions (1 km/10 km/0.25 Historical burned area reconstruction for the 20th century 20 To evaluate the simulated burned area for the 20th century, historical burned area data were used. These data, which cover the period 1900-2000, were compiled by Mouillot and Field (2005) at 1 • resolution and monthly time step (hereafter referred to as the Mouillot data). The data were generated by first synthesizing the burned area information from published data at national or regional scale for the periods of the 1980s Introduction

Tables Figures
Back Close

Full Screen / Esc
Printer-friendly Version

Interactive Discussion
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | or 1990s, further interpolated spatially at the global scale at 1 • resolution using the available satellite-derived active fire distribution. Then national fire statistics, historical land-use practices and other fire-relevant quantitative information (such as tree ring reconstruction) were used to build the historical fire temporal trend to interpolate historical burned area.

5
When comparing the burned area given by GFED3.1 and the Mouillot data for their overlapping period of 1997-2000, the global total burned area given by Mouillot and Field (2005) is 52 % higher than GFED3.1 with significant regional discrepancies. As the satellite-derived data is considered to be more reliable than the national or regional statistical data at a large spatial scale, a bias correction was performed on the Mouillot 10 data. We calculated the ratio of burned area by the Mouillot data to GFED3.1 for 1997-2009 for each region (Fig. 2b) and also the globe. This ratio was then applied to correct for each decade the burned area data in the Mouillot dataset.
Note that the regional breakdowns of the globe by GFED3.1, and Mouillot and Field (2005) are different (Fig. 2). For comparison of burned area for the 20th century, the 15 regional breakdown by Mouillot and Field (2005) was adopted as it is based on maximum temporal stability (error consistency), a highly important factor when comparing long-term data.

Fire patch data
Alaskan and Canadian fire management agencies have maintained historical fire mon-20 itoring for a relatively long time (dating back to 1950s). Their observations provide robust datasets for model validation of fire trend and variability. The historical fire information for the US state of Alaska was retrieved from the Alaska Interagency Coordination Center (AICC, http://afsmaps.blm.gov/imf_firehistory/imf.jsp?site=firehistory). The fire information for Canada was from the Canadian Wildland Fire Information System 25 (http://cwfis.cfs.nrcan.gc.ca/ha/nfdb). These datasets contain information on fire location, fire size (burned area), fire date and fire cause. Note, in these two datasets fires with all sizes are included.  Archibald et al. (2010) classified the MCD45A1 500 m MODIS fire burned area data into individual fire patches for the southern African region. This fire patch information includes location, patch size (with minimum fire size of 0.25 km 2 ), and the start and end dates of burning. The fire patch data for boreal North America and southern Africa are used to evaluate the model simulated fire size distribution (Fig. 3). Global observation datasets for fire season length, the 95th quantile fire size and 95th quantile fire intensity are provided by Archibald et al. (2013). The fire season length was quantified as the number of months required to reach 80 % of the total average annual 10 burned area by using GFED3.1 data. The MCD45A1 burned area product at 500 m resolution was used to derive the individual fires by applying a flood-fill algorithm, and the 95th quantile fire size in each grid cell was extracted to represent the size of large fires. The fire intensity is approximated by fire radiative power (FRP), which is derived from the satellite middle-IR wavelength measurements sensed over actively burning 15 fires. The MCD14ML active fire product was used to calculate the 95th quantile FRP as a proxy for the fireline intensity of large fires. The fire season length simulated by the model was also quantified as the number of months required to reach 80 % of the simulated mean annual burned area for the period of 1997-2009. As a homogeneous set of fires with a single fire size is simulated 20 by the model over a given grid cell on a daily time step, there is no corresponding fire patch information available on each grid cell. Thus the simulated daily mean fire size and fireline intensity over a given grid cell during the period of 1997-2009 were pooled together to create a "pool of fire patch", which was further used to extract the 95th quantile of fire size and fireline intensity. As the GFED3.1 data are most widely used by the fire modelling community, the model results are evaluated against GFED3.1 data for the period of 1997-2009. Three aspects were examined: the mean annual value of burned area, the similarity in interannual variability and the seasonality. The evaluation has been done for each GFED3.1 5 region (Fig. 2a).
For the model error in terms of mean annual burned area (BA), we use the relative difference: where BA model is the simulated burned area averaged over 1997-2009, and BA GFED is GFED3.1 mean annual burned area for the same period. The similarity in interannual variability (S interannual ) is estimated by the correlation coefficient of the two linearly detrended monthly burned area time series by model simulation and GFED3.1 data. Finally, the seasonality similarity (S season ) is given by: where frac_model i and frac_GFED i are the fraction of burned area for the i th month relative to the annual burned area. Fires in grassland-dominated systems are well captured by the model, including steppe fires in central and east Asia, savanna fires in northern Africa, northern Australia 15 and central to east South America. Two regions could be identified where model simulation is different from all the three observation datasets. One is the woodland savanna (miombo) in southern Africa, where burned area is underestimated by the model (simulated annual burned fraction is ∼ 4 %, but 14-24 % is observed). The other is western and central continental US (dominated by C3 and C4 grass in the land-cover map used 20 by ORCHIDEE) where fires are overestimated by the model (simulated annual fraction is ∼ 6 %, but 1-2 % is observed). For the fires in high-latitude (> 45 • N) boreal forest and tundra, the magnitude of burned fraction by ORCHIDEE falls between that found in the GFED3.1 and in the L3JRC/GLOBCARBON datasets (Fig. 5).

Results
The simulation pixels are divided into five classes according to their simulation qual- 25 ity, as shown in Fig. 4e. Table 1 shows the mean burned area for each category and dataset. The grid cells with burning collocated for both the model and observation data (labelled as ORC-good, ORC-max, ORC-min in Fig. 4)  can reproduce the major spatial extent of burning. However, discrepancy still remains, in that 50 % of the modelled burned area is classified as ORC-max (i.e., overestimation of burned fraction by the model), whereas observation datasets have half of the burned area labelled as ORC-min (i.e., underestimation of burned fraction by the model). Shifting from the single value of a(ND) (0.22 over the globe in Eq. 1) to a spatial dataset did 5 not improve the result (Supplement Table S1). Figure 5 also illustrates the lower burned area found in L3JRC and GLOBCARBON in comparison with GFED3.1 for the Southern Hemisphere and subtropical Northern Hemisphere, in contrast to the higher burned area in the middle-to-high latitude region in the Northern Hemisphere. The simulated latitudinal distribution of burned area generally falls within the minimum-maximum range of the three observation products.
Exceptions are the regions of ∼ 5-15 • S and 30-40 • N, corresponding to the underestimated burning in southern African savanna and the overestimate in western and central US, discussed above.
The interannual (12-months smoothed global total burned area) time series of 15 burned area by the GFED3.1, model simulation and the GLOBCARBON data are quite similar, but the J3JRC data has a different temporal phase (Fig. 6). The model result agrees best with the GFED3.1 variability, except that the model fails to reproduce the 1998 peak burning. Shifting from the single a(ND) value (0.22 over the globe) to the spatial dataset (Supplement Fig. S1) improves the agreement of the simulated global 20 total burned area with the GFED3.1 data. This indicates that the model is able to reproduce the temporal trend of burning on the global scale as observed by the satellite data.

Model evaluation against GFED3.1 burned area data
The simulated burned area for 1997-2009 is evaluated against the GFED3.1 data for 25 each region (shown in Fig. 2a)  presented in Table 2. The comparison of annual burned area at regional scale is shown in the Supplement Fig. S4. The correlation coefficient for the linearly detrended global total annual burned area between the ORCHIDEE simulation and GFED3.1 is 0.57, indicating that the model captures the interannual variability of burned area rather well, because errors are compensated among regions. On regional scale, the similarity of interannual variability and 15 seasonality between the model and GFED3.1 is generally good, with S interannual values bigger than 0.5 and S season bigger than 0.6 for most of the regions (Table 2).

Fire and precipitation relationship
The model captures well the empirical relationship between burned area and precipitation found in tropical and subtropical regions ( Fig. 7; see also Prentice et al., 2011;20 van der Werf et al., 2008b). In low precipitation regions (< 400 mm yr −1 ), the climate is favourable for fire but burning is limited by the available fuel. In contrast, regions with higher precipitation (> 2000 mm yr −1 ) always support a relatively large fuel amount but

Peak fire month and fire season length
The spatial distributions of the peak months for burned area by model simulation and as given by GFED3.1 are compared in Fig. 8. The peak fire month is defined as the month 10 of the year with the maximum monthly burned area. The global spatial pattern of simulated peak fire month is in general agreement with that given by GFED3.1 data. The difference between simulated and observed peak month is quantified by the following index, after Prentice et al. (2011): 15 where θ j is the angle between vectors representing the simulated and observed peak fire month (with January to December resembling one to twelve on a clock), n is the total number of grid cells and A j is the burned area by GFED3.1 data. According to Eq. (4), the value of D 2 is zero when simulated peak month is perfectly in phase with 20 the observation, 0.5 if the timing is off by 3 months in either direction, and one (the maximum) if the timing is off by 6 months. The model simulation gives D 2 = 0.51, indicating that the simulated and observed peak fire month differ by an average of three months. At regional scale, the simulated peak fire months for most regions are within one month of those given by the GFED3.1 data, with the exceptions of the Middle East and South-Introduction  Fig. S6 for more details of the comparison of modelled and observed seasonal pattern of burning for different GFED3.1 regions. Figure 9 compares the model results with the GFED3.1-derived fire season length from Archibald et al. (2013). The spatial pattern of the fire season length by model simulation agrees well with that given by Archibald et al. (2013), with the fire season length 5 lasting 1-3 months in the boreal region, and 4-7 months in semi-arid grasslands and savannas. The fire season length in the southern fringe of northern African savanna, and in the savannas of southern Africa and South America seems to be underestimated by the model by 2-4 months.

Long-term trends of burned area during the 20th century
Over the 20th century, the historical trend of modelled global total burned area generally follows the Mouillot reconstruction data (Mouillot and Field, 2005) as corrected by GFED3.1 data (see Sect. 2.4.1), with a gradually increasing burned area after the 1930s until the 1990s, after which global burned area began to decrease (Fig. 10). Regionally, simulated trends in burned area agree best with the Mouillot data in boreal 15 Russia. The simulated burned area also agrees well with fire agency survey data for the US state of Alaska and Canada (Supplement Fig. S7). This reflects the model capability to capture fire trends driven by climate variation relatively well. This agreement also suggests that the CRU-NCEP climate data can be used for simulating trends in fires despite the uncertainties. The model fails, however, to capture burned area trends 20 for regions where fires resulting from changed land use likely played a bigger role in the earlier 20th century according to Mouillot data (Mouillot and Field, 2005), e.g., Australia and New Zealand, USA and southern South America, or in regions where the implementation of modern fire prevention has drastically reduced the burned area, as occurred in the 1960s in boreal North America.

Fire size distribution in boreal North America and southern Africa
Above a certain threshold, the relationship between the number of fires and fire patch size conforms to the power-law distribution (Malamud et al., 1998(Malamud et al., , 2005 and can be described as: where N is the number of fires, S is the individual fire patch size, β is the slope on a log-log plot. In this section we compare the simulated fire size distribution against observational data over two regions: boreal North America and southern Africa (see map in Fig. 3).

10
The simulated daily fire number and mean fire size for the same period as the observation data were used to generate the modelled fire patch data. For both simulated and observed fires, all fires within the test region were pooled together. Fires were binned according to the fire size in an equal logarithmic distance manner (with the minimummaximum size range being divided into 100 bins). The mean annual number of fires on 15 an area basis for each bin was calculated.
The minimum fire size for southern Africa given by Archibald et al. (2010) is 25 ha (which is the resolution of a 500 m MODIS burned area pixel) and this value was arbitrarily used as the fire size threshold in the regression analysis for both regions. The mean value of each size range is used as the explaining variable. Both fire number 20 and fire size data are log transformed, and OLS regression is used to estimate the parameter values (slope beta, intercept alpha).
For both boreal North America and the southern Africa, the slopes obtained from the model simulation are higher than the ones derived from fire management agency data (boreal North America) or satellite-derived data (southern Africa) (Fig. 11). This 25 indicates that the simulated fire distribution is skewed towards small fires and that the number of big fires is underestimated. However, the model successfully reproduced the power-law (heavy-tailed) distribution pattern of fires. Introduction We further calculate the cumulative fraction of the total burned area by fires below a given quantile of fire size (the minimum size, every tenth quantile from 10th to 90th quantile, and the maximum size) (Fig. 12). According to the observation data, in boreal North America, the total burned area is mainly dominated by a few big fires, with the top 10 % of fires (90th quantile to the maximum size) accounting for 99.8 % of the 5 total burned area. By contrast, the same group of fires (i.e., the highest 10 % big fires) account for only ∼ 60 % of the simulated total burned area, with the remaining being accounted for by many small fires. For southern Africa, the satellite observation data show the top 30 % of fires (70th quantile to the maximum size) make up 90 % of the total burned area, whereas the model simulates 90 % of the total burned area coming 10 from half of all fires (median to the maximum size). Figure 13 compares the simulated 95th quantile fire size with the global observations. According to the observation data, the biggest fires (500-10 000 km 2 ) are grassland- The FRP quantifies the energy released per unit of time by an active fire based on the pixel area of the satellite, while the fireline intensity is the rate of heat transfer per 25 unit length of the fireline (kW m −1 ) (Byram 1959). Thus the FRP can only function as a temporally and spatially continuous proxy for the fireline intensity over the globe, as the length of flame front is needed to convert the former into the latter ( Smith and Wooster, 2005). So the comparison here between the simulated 95th quantile fireline intensity with the 95th FRP is mainly for qualitative purposes. As shown in Fig. 14, at the global scale, the model is able to reproduce the spatial pattern of the fire intensity revealed by the FRP data, although with rather different absolute values due to the nature of the two different datasets. Fires with high intensity 5 occur prominently in two different ecosystems over the globe, one is the grasslanddominated system in central and eastern Asia, inland Australia, and the southern Africa and America. The other is the boreal forest system in the US state of Alaska and northern Canada, and eastern Russia. The model captures this pattern relatively well, despite the fact that the highest fire intensity is simulated to occur in Arctic tundra.

General model performance
The SPITFIRE module was first presented by Thonicke et al. (2010). It was later adapted in LPX by Prentice et al. (2011), notably with the removal of ignition sources created by humans. In the current study, the SPITFIRE module was fully integrated 15 into the global process-based vegetation model ORCHIDEE for the first time. Simulated burned area, fire season, fire size distribution and fire intensity were evaluated in a comprehensive way against observation data for a short recent period (1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009), and against reconstructed burned area for the 20th century.  Prentice et al., 2011;van der Werf et al., 2008b), indicating that precipitation, as an 5 important driver for fire, has been successfully captured by the model. For the boreal region where climate plays a dominant role in determining the long-term trend of burned area, the simulated burned area generally agrees well with the historical reconstruction data (boreal Russia for 1920-2009) and the government fire agency data (US state of Alaska andCanada for 1950-2009), indicating that the model is capturing the climate 10 as a driver of fire. However, on the global scale because the fire trend is determined by multiple factors including climate, land-use practice and fire suppression (Mouillot and Field, 2005), the simulated burned area trend only agrees moderately well with the reconstruction data for the 20th century. 15 Fire is a complex, regional process that involves dimensions of vegetation, climate and human activity (Bowman et al., 2009). The uncertainties in all these factors will contribute to the overall uncertainty in simulated fire activities. Because SPITFIRE simulates fire occurrence through a complex chain that includes multiple processes from potential ignition to fire climate, to fire spread rate and tree mortality, identifying the 20 contributions of each modelling step to the ultimate error in the simulated fire regime is problematic. A complete error analysis involving all the model parameters is beyond the scope of this study, but the following sections are intended to serve as preliminary investigation of model errors.

Ignition sources and fire numbers
On the global scale, due to the limitation of fire by fuel load on the dry regions, the modelled annual burned area is more closely related to the total fire numbers rather than the fire danger index (Suplement Fig. S8). This might lead to speculation that the potential ignition source is the first identified source of error for the simulated burned area. 5 However, the shift from a constant uniform value of a(ND) = 0.22 in the human ignition equation to the spatial dataset (Supplement Fig. S1) does not significantly improve the regional model-observation agreement (Supplement Fig. S2). The other possible error in ignition is that potential lightning ignitions are not suppressed by humans in densely populated areas, which might cause lightning-ignited fires being overestimated. We 10 have tested this possibility by applying a population density dependent human suppression of lightning-ignited fires following Li et al. (2012), and the result showed that part of the overestimation of burned fraction in western US and central South America could be reduced ( Fig. 4 and Supplement Fig. S9), but the burned area in Africa was further underestimated. 15 Published fire models (Kloster et al., 2012;Li et al., 2012;Pechony and Shindell, 2009;Venevsky et al., 2002) generally include ignitions from both lightning and human sources, together with explicit or implicit human suppression of fires. However, one common challenge is to properly calibrate ignition parameters. One option is to use active fire counts (as in case of Li et al., 2012;Pechony and Shindell, 2009), but 20 fire counts are not exactly real fire numbers, because a single widespread fire could be seen as many fire counts and the burned area per hotspot could vary by an order of magnitude depending on vegetation composition (Hantson et al., 2013). Little information is provided in the literature on the comparison of simulated and observed fire numbers. 25 We compared the simulated number of fires with the observed fire patch number for Canada and southern Africa (Supplement Fig. S10). Simulated fire numbers are higher than observation in both regions, with an overestimation of ∼ 70 % in Canada and ∼ 140 % in southern Africa. However it should be noted that the observation data are also subject to error, as very small fires might be missed in the Canadian fire agency survey data. The fire patch data given by Archibald et al. (2010) were reclassified from MODIS 500 m burned area data and as a result the fires smaller than 0.25 km 2 threshold will be missed. To make the model practically useful for burned area prediction at 5 the regional scale, in principle one could use the ratio of observed to the simulated burned area to correct the potential ignition numbers when there is a lack of observed fire patch data (ideally separated by different sources).

Fire size
Our simulated burned area for Canada agrees well with the fire agency survey data 10 (Supplement Fig. S7), but simulated burned area for southern Africa is 60 % lower than GFED3.1 data (Table 2). Given the overestimation of the number of fires, the mean fire size is underestimated by the model. The fire size distribution analysis reveals that the model mainly underestimates the number and size of the large fires in these two regions. Over the globe, despite the fact that the model correctly identifies some of the 15 regions where large fires occur (mainly with grassland fraction higher than ∼ 70 % by the land cover used in the model), the large fire size is grossly underestimated in the model -by up to two orders of magnitude (Fig. 13). Within the model, the fire size is determined by the fire duration time and fire spread rate. Fire spread rate is determined using the Rothermel equation (Rothermel, 1972;20 Wilson, 1982) and depends on multiple factors including wind speed, fuel bulk density, net fuel load and fuel packing ratio, etc. To gain some insight into the model's behaviour, we tried to relate the 95th quantile fire size to some parameters (grassland fraction, fuel bulk density). As shown in Supplement Fig. S11, the size of large fires is exponentially dependent on the fire spread rate. The fire spread rate is very sensitive to the fuel bulk 25 density and grass fraction beyond some threshold (e.g., fire spread rate surges when grass coverage exceeds ∼ 70 %), with the fuel bulk density being inversely dependent on the grass fraction. Thus the simulated fire size could be sensitive to the land-cover 2401 Introduction

Tables Figures
Back Close

Full Screen / Esc
Printer-friendly Version

Interactive Discussion
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | map (grass fraction) used in the simulation. As a static land cover map is used in our simulation, the grassland fraction is not allowed to vary in line with the fire disturbance and thus a full fire-climate-vegetation feedback is limited, this could help to explain the underestimation of fire size.  Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | a 0.5 • grid cell, the model simulates (on a daily basis and when fire does occur on a specific day) a given number of homogeneous fires with a single fire size, while in reality a big fire might grow from a single or multiple ignitions (in the case of separate fires merging with each other) over a period of days or even weeks (Keeley et al., 2009). The daily time step of the model does not allow fires to span more than one day.

5
To fully represent the big fire process in reality, improvements need to be made to the model to allow fire to span multiple days when the climate is favourable. More mechanistic process modelling is needed, for example, by incorporating the information on fuel, weather and topography to predict the fire extinction. 10 In the current simulation the dynamic vegetation module of ORCHIDEE was switched off and a static land-cover map was used. Tree mortality was affected by fire-induced tree damage, but tree coverage within a given grid cell was static and not allowed to vary with fire occurrence. A test simulation has been done for Southern Hemisphere Africa following the same simulation protocol as in Sect. 2.3 but with the dynamic veg- In dynamic vegetation mode, the simulated burned area suddenly begins to increase 20 around 1965, in response to increased fire danger index. The increase in fire activity further increases the grass coverage and reduces tree coverage, causing a positive feedback to finally induce a peak of burned area around 1975, after which the burned area decreases. In static vegetation mode, the simulated burned area shows a similar peak in response to the peak in the fire danger index, however, with a much smaller 25 peak of burned area than that simulated in dynamic vegetation mode because of the lack of fire-vegetation-climate feedback. The simulated burned area by both simulations is still lower than that given by the GFED3. ing the competitiveness of trees vs. grass might also play a role in the error of fire modelling.

Summary
The preliminary investigation of modelling error reveals that fire size is generally underestimated by the model over the entire globe and the ignition error is playing an im-10 portant role in determining the ultimate simulated burned area. On the regional scale, ignition numbers (fire numbers) are either overestimated to compensate the fire size underestimation to cause a moderate or overestimated burned area, or are not enough that the simulated burned area is underestimated as well. The underestimation of fire size is most likely due to the constrained maximum fire duration time, although in some 15 regions the fire spread rate is also underestimated relative to measured data.

Future model improvement directions and needed datasets
Currently many efforts in global fire modelling are directed at reproducing the temporal and spatial pattern of burned areas (Kloster et al., 2010;Li et al., 2012;Prentice et al., 2011;Thonicke et al., 2010). Total burned area is determined by ignition frequency and 20 fire size, which itself is controlled by fire spread rate (fire intensity) and fire duration. More work is needed to investigate if a model can reproduce the mechanisms that drive burned area: i.e. the rate of spread, fire duration, fire size, ignition frequency, and fireline intensity. Comparing observed and simulated fire regimes, which combine information on fire timing (fire season), size, numbers and intensity (Gill and  2008) will help to reveal gaps in this understanding. The present study is a step in this direction, bringing new in-depth model evaluation.
In summary, the fire processes in the SPITFIRE model are complex enough to include many aspects of wildland and human-caused fire processes in nature. However, little is known about the parameter sensitivities and their contribution to model error. 5 The simulated intermediate model parameters (e.g., fire numbers, fire size, fire spread rate, fire intensity) are poorly constrained by the observation data. As a result, error compensation could be prevalent in the model and a wider application of the model is impeded.
To advance model development, global measurement datasets of the key fire pa-10 rameters, including fire spread rate, fuel bulk density, wind speed, fire intensity, fire duration, etc., should be established and used to calibrate fire models. On the other hand, the complexity of fire model parameters and the regional nature of fire processes make it unlikely that these parameters could be calibrated in a parameter-by-parameter and site-by-site way, but some more advanced techniques such as data assimilation or 15 model-data fusion could be helpful. Finally, we should consider incorporating some more mechanistic fire processes into the model, such as crown fire spread, impacts of land fragmentation, multi-day fire duration, topography effects, and the mechanistic process of fire extinction.