Modeling Burned Areas in Indonesia: The FLAM Approach

Abstract

Large-scale wildfires affect millions of hectares of land in Indonesia annually and produce severe smoke haze pollution and carbon emissions, with negative impacts on climate change, health, the economy and biodiversity. In this study, we apply a mechanistic fire model to estimate burned area in Indonesia for the first time. We use the Wildfire Climate Impacts and Adaptation Model (FLAM) that operates with a daily time step on the grid cell of 0.25 arc degrees, the same spatio-temporal resolution as in the Global Fire Emissions Database v4 (GFED). GFED data accumulated from 2000–2009 are used for calibrating spatially-explicit suppression efficiency in FLAM. Very low suppression levels are found in peatland of Kalimantan and Sumatra, where individual fires can burn for very long periods of time despite extensive rains and fire-fighting attempts. For 2010–2016, we validate FLAM estimated burned area temporally and spatially using annual GFED observations. From the validation for burned areas aggregated over Indonesia, we obtain Pearson’s correlation coefficient separately for wildfires and peat fires, which equals 0.988 in both cases. Spatial correlation analysis shows that in areas where around 70% is burned, the correlation coefficients are above 0.6, and in those where 30% is burned, above 0.9.

1. Introduction

Fires affect terrestrial ecosystems and have profound consequences on global climate, air quality and vegetation [1,2]. Fire regimes are determined by climate variability, vegetation and direct human influence. Seasonal weather is the most variable factor and has the greatest effect on regional burned area [3,4]. Climate change lengthens fire seasons, with an associated trend towards increasing extent, frequency and severity of fires, as well as the global occurrence of extreme (mega) fire events [5] in many regions of the Earth (e.g., in the Mediterranean [6] and boreal [7] regions). At the same time, human activity is the primary source of ignition in global forest, savanna and agricultural regions. Indonesia is one of the countries that has experienced particularly large-scale peatland and forest fire disasters [8]. According to the World Bank [9], 2.6 million hectares of Indonesian land burned between June and October 2015, costing Indonesia USD 16.1 billion.

Most fires in Indonesia are induced by human activity [9]. Forest and land clearing for plantations and agriculture have been cited by various researchers as the major causes of these fires [10]. Indonesia contains 47% of the global area of tropical peatland, located mostly in peat swamp forest [11]. Large peatland areas have been converted to agriculture, initially under the Indonesian Government’s transmigration program, and later for commercial oil palm and paper pulp tree plantations, especially in Sumatra and Kalimantan [12,13]. Draining and conversion of peatland, driven largely by palm oil production, contributes to the intensity of haze from fire [9]. Peatlands account for about one-third of the area burned, but are responsible for the vast majority of the haze and CO${}_2$ emissions from the fires [9].

Fires have a large disruptive impact in Indonesia, affecting millions of hectares of land, producing severe smoke haze pollution. There are multiple negative effects of wildfires in Indonesia on climate change [14], health [15], the economy [8,9] and biodiversity (http://fires.globalforestwatch.org/about/). Indonesia’s biomass burning emissions are comparable to the annual fossil fuel CO${}_2$ emissions of countries such as Japan and India [16,17].

The overall fire situation in Indonesia is described in various papers [18,19]. In this study, we apply a mechanistic wildfire model to estimate burned area in Indonesia. Fire models are developed for hazard assessment, as well as forest and fuel management planning, including prescribed burning and suppression efforts. The paper led by Herawati [20] reviewed available fire models (e.g., the Risco de Queimada/fire risk model (RisQue), Coupled Atmosphere-Wildland Fire-Environment (CAWFE) and FireFamilyPlus), which could be potentially used in Indonesia. However, we have not come across a mechanistic model that specifically focuses on modeling burned areas in the forest and peatland of Indonesia in the literature. This motivated us to apply the Wildfire Climate Impacts and Adaptation Model (FLAM) to modeling burned areas in Indonesia, as described in this paper.

There have been studies from different territories where historical climate data in combination with Moderate Resolution Imaging Spectroradiometer (MODIS) satellite data were used in fire mapping [21,22]. However, they did not deal with explicit modeling of burned areas. At the same time, area burned is the most commonly-used indicator when fire trends are being examined [23], and the assessment of burned area is an important step in estimating the cost of fire [9], as well as the air quality impacts [24]. Previously in studies led by Khabarov [25] and Krasovskii [26], FLAM (formerly known as the Standalone Fire Model (SFM)) was applied to modeling burned areas in Europe. In this study, we enhance the model by estimating several wildfire process parameters to better fit regional conditions in Indonesia. For the purposes of this paper, we separate peat fires from wildfires.

Methodologically, our paper contributes to the development of the original algorithm by Arora and Boer [27]. In a recent study, Arora and Melton [28] linked fire extinguishing probability with population density. We demonstrate that this simplified approach can be improved by an explicit calibration of suppression efficiency.

In this paper, we test the following hypotheses: (1) the observed burned area in Indonesia reported by the Global Fire Emission Dataset (GFED) [29] (which is based on satellite observations) can be modeled independently by FLAM; (2) the approach integrating the Canadian Fire Weather Index (FWI) [30,31,32] into a mechanistic fire model provides good accuracy in assessment of burned area; (3) the spatial suppression efficiency calibrated using observed burned area improves the original fire algorithm for area burned by wildfires and peat fires.

2. Materials and Methods

FLAM is used to reproduce historical and to project future burned area, as well as to assess climate change impacts and adaptation options (http://www.iiasa.ac.at/flam). FLAM is able to capture impacts of climate, human activity and fuel availability on burned areas. It uses a process-based fire parametrization algorithm that was developed to link a fire model with the Dynamic Global Vegetation Models (DGVMs) in the literature [27,28,33,34,35,36]. FLAM is a standalone fire model based on the Arora and Boer [27] algorithm, being decoupled from DGVM. This allows for faster running times and application of the calibration approach, described below. The key features implemented in FLAM include fuel moisture computation based on the Fine Fuel Moisture Code (FFMC) of the Canadian FWI [30] and a procedure to calibrate regional fire suppression efficiency as proposed in previous studies [25,26].

The FLAM flowchart is depicted in Figure 1. The model operates with a daily time-step at 0.25-arc degree spatial resolution. All inputs were re-sampled to fit the model resolution, using the resample and aggregate functions of the raster package [37] in the R software; we used the country shape of Indonesia from the wrld_simpl dataset available in R.

To model daily fire occurrences, FLAM computes three conditional probabilities of fire:

• probability of ignition depending on natural and human sources
• probability of ignition subject to fuel availability
• probability of ignition depending on weather

2.1. Ignition Probability

There are two causes of ignition: (i) natural and (ii) human activity. Natural ignition is represented in the model by lightning frequency. Our approach to modeling human ignition and suppression consists of two steps. Firstly, we use a standard procedure, where population density is transformed into probability functions of human ignition and simultaneously suppression [35], i.e., “unsuppressed ignition”. Secondly, we apply a calibration procedure that identifies pixels with higher or lower suppression efficiency [26], i.e., probability of extinguishing a fire within one day, based on observed burned areas. In the case when the calibrated probability in a pixel equals one, the fire is suppressed “immediately”. Thus, observations help us to additionally identify the human activity impact independently from population density.

In this study, population density was taken from the dynamic maps of the Gridded Population of the World, v4 [38]. We used maps for the years 2000, 2005, 2010 and 2015, inside the modeling interval 2000–2016. Original data were averaged using the mean function excluding empty data cells. The corresponding human ignition probability ($P_h$) is implemented in FLAM using the equation that relates fire occurrence to population density $\rho$:

$$P_h=min\left(1,{\left(\frac{\rho }{{\rho }_{up}}\right)}^{0.43}\right),$$

(1)

where a threshold value (${\rho }_{up}$) was set to 300 inhabitants/km${}^2$ [35]. The suppression probability ($F_{supp}$) depends on the population density in the following manner. More likely, fire suppression is provided in densely-populated areas, where typically high property values are at risk and more resources and infrastructures are available for suppression [35,36]. The probability function is parameterized as follows [39]:

$$F_{supp}=1-(0.1+exp(-0.025\rho )).$$

(2)

Equation (2) assumes that fire suppression increases for low values of $\rho$, and then declines, reflecting the fact that in more densely-populated areas, 90% of the fires are almost immediately suppressed [35].

Lightning represents a natural ignition source in the model. Here, we use monthly climatology averaged from 1995–2014 at 0.5-degree resolution from LIS/OTD Gridded Lightning Climatology Data Collection Version 2.3.2015 [40]. According to a recent study [41], excluding the inter-annual variability of monthly lightning forcing has only a minor affect on simulated burned area. The LIS/OTD 0.5 Degree High Resolution Monthly Climatology (HRMC) contains gridded data of total lightning flash rates obtained from two lightning detection sensors: the spaceborne Optical Transient Detector (OTD) on Orbview-1 and the Lightning Imaging Sensor (LIS) onboard the Tropical Rainfall Measuring Mission (TRMM) satellite. The long LIS (equatorward of about 38 degree) record makes the merged climatology most robust in the tropics and subtropics. Data were resampled to 0.25-arc degree resolution.

The algorithm used in FLAM is based on the assumption that cloud-to-ground lightning frequency ($L_f$ [flashes/km${}^2$/month]) varies from a specified lower value of essentially no lightning ($L_{f,low}=0.02$ flashes/km${}^2$/month) to an upper value close to the maximum observed ($L_{f,up}=0.85$ flashes/km${}^2$/month) [27]. Therefore, the lightning frequency can be transformed to a value from 0–1 according to the equation:

$${\beta }_L=max\left[0,min(\frac{L_f-L_{f,low}}{L_{f,up}-L_{f,low}})\right].$$

The probability of lightning ignition ($P_l$) is modeled using the function:

$$P_l=\frac{{\beta }_L}{{\beta }_L+exp(1.5-6{\beta }_L)}.$$

(3)

Using Equations (1)–(3), the daily probability of unsuppressed ignition is calculated as follows:

$$P_i=(P_l+(1-P_l)P_h)(1-F_{supp}).$$

(4)

2.2. Fuel Availability and Weather

In this study, we used the tropical biomass map of [42], providing the Aboveground Biomass (AGB) density of vegetation in units of Mg/ha with the cell size of 0.008 arc degrees, corresponding to around one kilometer. The value for each pixel of the map was converted to gC/m${}^2$ by multiplying the original values by 50 [43] and aggregated using the R’s mean function under the exclusion of no available data cells.

To model fuel available for burning, we need information about litter and Coarse Woody Debris (CWD). We used the data from the Global Forest Resources Assessment 2015 (FRA 2015) of the Food and Agriculture Organization of the United Nations to obtain the shares of these components in the above-ground biomass (AGB). As these values were not available for Indonesia, we used the values reported for neighboring tropical countries. The dead wood information was available for Philippines, where the ratio of dead wood to AGB was $0.1358$. The litter fraction in AGB was also not available for Indonesia; therefore, an average of the neighboring countries (Laos, Malaysia, Myanmar, Philippines) was calculated, resulting in a litter-to-biomass ratio of $0.0328$. Biomass available for burning (B), which is the sum of litter and CWD, was transformed to the ignition probability conditional on fuel availability ($P_b$) according to the equation:

$$P_b=max\left[0,min\left(1,\frac{B-B_l}{B_u-B_l}\right)\right],$$

(5)

where $B_l$ and $B_u$ are constants set to 200 gC × m${}^{-2}$ and 1000 gC × m${}^{-2}$, respectively [36].

Climate is one of the main drivers of fire dynamics in FLAM. Climate change effects on wildfire risk are discussed in many studies [17,25,44,45]. A well-established Fire Weather Index (FWI), initially developed for Canada by Van Wagner and Pickett [30], is used to assess fire danger in many regions of the world, outside of Canada. The applications of FWI include European countries [46], as well as Indonesia and Malaysia [47]. Daily Fine Fuel Moisture Code (FFMC) values from the Global Fire WEather Database (GWFED) [48] were used as FLAM inputs in the period 2000–2016. GFWED is comprised of eight different sets of FWI calculations, all using temperature, relative humidity, wind speed and snow depth estimates from the NASA Modern Era Retrospective Analysis for Research and Applications Version 2 (MERRA-2) [49]. We used the version based on the observation-corrected precipitations [50]. The data are available from 1981 globally at a resolution of $0.5\times 2/3$ arc degrees. We used the data till 2016 for our modeling, as there was no complete data for 2017 in the dataset. The FFMC values were resampled to 0.25-arc degree resolution and converted to moisture content values, used in FLAM. Here, we applied the transformation derived from the FWI equation [30]:

$$m=\frac{0.01(59.5\times 250-147.2\times V_{FFMC})}{V_{FFMC}+59.5},$$

where $V_{FFMC}$ is the FFMC value and m is the fuel moisture content (ratio between 0 and 1). The daily probability ($P_m$) of ignition is calculated as:

$$P_m=1-tanh{(1.75\times m/m_e)}^2,$$

(6)

where $m_e=0.35$ is the moisture of extinction (ratio between 0 and 1) [33,35,51]. Note that $m=0.35$ corresponds to $V_{FFMC}\approx 70$.

Based on Equations (4)–(6), the total daily probability of fire is calculated by multiplication:

$$P=P_b\times P_m\times P_i.$$

(7)

2.3. Fire Spread Rate and Expected Burned Area

Fire spread rate is limited by the moisture content and wind speed. The daily wind speed data were taken from the MERRA-2 dataset and resampled to 0.25-arc degree resolution. Fire spread rate in (km/h) in the downwind direction is calculated using the equation [35]:

$$u_p=g\left(w\right)\times h\left(m\right)\times u_{\max },$$

where w is wind speed (km/h) and $u_{\max }=0.45$ km/h is the maximum spread rate. Here, function g describes the dependence on the wind speed:

$$g\left(w\right)=0.1exp(-\frac{w^2}{2500}),$$

while function h, the dependence on moisture content:

$$h\left(m\right)={(1-{\beta }_m)}^2,\mathrm{where}{\beta }_m=min(1,\frac{m}{m_e}).$$

The total area burned is calculated as follows:

$$a=0.36\pi \frac{u_p^2}{l_b},$$

(8)

where $l_b$ is the length to breath ratio of an ellipse approximating burned area:

$$l_b=1+10(1-exp(-0.017w_s)).$$

Assuming that a particular fire is put out on a given day with probability q, the expected area burned within one day can be calculated as follows:

$$A_{burn}=P\times a(1-q)(2-q)q^{-2},$$

(9)

where $q\in (0,1]$ is the probability of extinguishing a fire in one day and P is the total probability of fire (Equation (7)).

2.4. Calibration Procedure Implemented in FLAM

In FLAM as in the original fire algorithm [27], the efficiency of fire suppression is defined as the probability q of extinguishing a fire within one given day (see Figure 1). Potential area burned within one day (Equation (9)) and the cumulative burned area over any time period of N days for a grid pixel are calculated as follows:

$$A\left(q\right)=\frac{\alpha (1-q)(2-q)}{q^2},\alpha =\sum _{j=1}^N{P}_j{a}_j,$$

(10)

where coefficient $\alpha$ reflects the availability of fuel, ignition sources and weather conditions, but is not a function of q [27,35], $P_j$ is probability (7) for a particular day j and $a_j$ is calculated for the same day j using Equation (8). In our calibration procedure, we find a value of the variable $q=q_c$ such that $A(q_c)=A_{obs}$, where $A_{obs}$ is the observed cumulative burned area in a grid cell over a given time period [25,26]. Based on a non-calibrated model run with an arbitrary value of $q=q_0$ ($0<q_0<1$) delivering accumulated burned area $A(q_0)$ for a time period for a given grid cell, the calibrated value $q_c$ is defined by the following equation:

$$q_c=\frac{-3+\sqrt{8\beta +1}}{2(\beta -1)};$$

(11)

here:

$$\beta =\frac{A_{obs}}{A(q_0)}(1-q_0)(2-q_0)/q_0^2.$$

Substituting the value $A(q_0)$ (10) into Equation (11), we get that parameter $\beta$ equals $A_{obs}/\alpha$, and therefore, the calibrated value of suppression efficiency $q_c$ does not depend on the arbitrary selected value $q_0$. Nevertheless, model runs with a determinate (place-holding) value for $q_0$ (e.g., 0.5) are necessary to obtain the value of $\alpha$. We apply this calibration procedure at the pixel level, forcing the model to fit the total accumulated burned area (as reported in GFED datasets) over the calibration period, which is long enough relative to the model’s operating daily time step. Let us note that calibrated values are sometimes considerably lower than the $0.5$ threshold, discussed by Arora and Melton [28], corresponding to the average duration of fire in the absence of human influence as 1 day. The expected fire duration is calculated according to the equation [28]:

$$\overline{\tau }=\frac{1-q}{q}.$$

(12)

Below, we will show that in peatland, the values of $q_c$ could be below 0.2, meaning that in these areas, average fire duration, $\overline{\tau }(q_c)$, could be much longer than 1 day, i.e., up to 173 days depending on the location (see the modeling results in Section 3.2).

The calibrated values of $q_c$ are then used in FLAM for calculating daily burned areas, i.e., (9) is modified as:

$$A_{burn}=P\times a(1-q_c)(2-q_c)q_c^{-2}.$$

Below, we describe three special cases arising during the calibration procedure.

2.4.1. Case I: Probability of Suppression Equals One

In the case when the accumulated burned area in a grid cell modeled by FLAM during the calibration period is above zero and, at the same time, the accumulated burned area reported in GFED in the same grid cell is zero, the value of parameter $q_c=1$. This means that although in this pixel, the total probability of fire (Equation (7)) is positive and results in modeled burned area, the calibration procedure identified that in this pixel, fire is not possible, i.e., the fire was never observed during the calibration period. In terms of modeling, this means that the fire is suppressed immediately after ignition. This fact could have different interpretations: firstly, this could be a correcting factor to the probability of human ignition discussed in Section 2.1; secondly, this could also be a proxy of a change in land cover, zeroing the probability conditional on biomass, which in turn is also a result of human activity, decreasing fire occurrence [52].

2.4.2. Case II: No Burned Area in Both GFED and FLAM

In the case when in a grid cell, both GFED and FLAM report zero burned area, the calibration procedure is not able to set a value of $q_c$ for this grid cell. To avoid indifference in choosing $q_c$, we keep the value of $q_c=0.5$ in those pixels in order to validate FLAM in the validation period. Note that this case is rare for Indonesia.

2.4.3. Case III: Not Feasible for Calibration

During the calibration procedure, there could be pixels where GFED reports burned area, while FLAM does not. In this case, the calibration procedure is not feasible mathematically. We identify those pixels during the calibration procedure, explain the reasons for their occurrence and exclude them from analysis as described in Section 3.2.

2.5. GFED Data Used in the Study

By applying historical climate data, we compare FLAM results with observed burned areas from the most recent Global Fire Emission Database GFEDv4 [29] for the period 2000–2016. The starting year 2000 is chosen due to the MODIS era starting in GFED in 2000, and the end year is determined by the availability of GFED data. In the study, the time period is divided into two intervals: the period 2000–2009 is used for FLAM calibration, while the period 2010–2016 is used for FLAM validation. In the calibration period, we use GFED data accumulated from 2000–2009. It is important to note that calibration and validation are two independent processes: for the period 2010–2016, FLAM is not calibrated against data or information from GFED; rather, GFED is used only for validation purposes, i.e., the comparison of burned area and location with the output of FLAM.

Currently, two versions of GFED are available: the one that includes small fires and another one that does not include them [53]; these are compared for Indonesia in Figure A1. One can see that small fires considerably increased the reported burned area, while the shape of the inter-annual variability stayed similar to the version without small fires. In this study, we use the version without small fires for the following reasons. Firstly, it includes the year 2016 (Version 4.1s is till 2015). Secondly, it contains information from the University of Maryland (UMD) land cover classification scheme applied to burned area and provides ratios of area burned in each of those classes, as well as the ratios of area burned in peatland. Note that GFED4.1s contains a split only for emissions. In our study, we model wildfires, and therefore, we choose burned area in the following 10 classes (out of 18): evergreen needleleaf forest, evergreen broadleaf forest, deciduous needleleaf forest, deciduous broadleaf forest, mixed forests, closed shrublands, open shrubland, woody savannas, savannas, grasslands [54]. Additionally, we use information on uncertainty intervals provided in GFED4, referred to as GFED in this paper.

A schematic illustration of how burned area is presented in GFED is shown in Figure 2a. Each grid cell contains information about total burned area. Additionally, we use the information of the fraction of area of wildfire (10 UMD land classes as mentioned above) and the fraction of area allocated to peat fire. Those fractions are independent, meaning that area burned by wildfire and peat fire can intersect in some areas, as indicated in the figure. Figure 2b shows actual data extracted for Indonesia from GFED separately for the total area burned, in wildland and peatland. The sum of areas burned in peatland and wildland is higher than the total burned area.

In our modeling, we use a similar approach. FLAM is calibrated separately for wildfire and peat fire, i.e., there are different suppression efficiencies. This approach allows us to distinguish between smoldering fires in peatland and flaming fires in forest.

2.6. Peatland Ratio Map

To calculate the fraction of area burned in peatland, we use a map produced by the Indonesian Ministry of Agriculture [55]. This map provides information on total share of peat in a grid cell in addition to burned peat within total burned area available in GFED. Original map indicated values of 1 for peat and 0 for no peat in each pixel at 30-m resolution. Based on these data, we calculated the ratio of peatland in each pixel at the 0.25-arc degree resolution as shown in Figure 3. The GFED dataset contains information on peat fires for Kalimantan and Sumatra. Therefore, in our analysis, we use the peatland map only for these islands. Nevertheless, we show the entire map with Papua greyed out.

3. Results

In this section we present the results of the FLAM model applied to Indonesia in their temporal and spatial dimensions. This is comprised of results for the model calibration and validation.

3.1. Temporal Validation of FLAM Modeling Results

Figure 4 shows the fit of annual burned area during wildfires over the entirety of Indonesia. In the calibration period, observed burned area in GFED (black line), accumulated from 2000–2009, was used for FLAM calibration (see Section 2.4). In the period 2010–2016, annual GFED data were used for validation of the burned area estimated by FLAM (red line). The figure demonstrates that FLAM fit actual observations quite well. The red line follows the black line and usually is in the range of indicative uncertainty intervals of GFED. Here, we used burned area uncertainty (in ha) reported in GFED [29].

We compared burned areas modeled by FLAM and reported by GFED using two indicators: Pearson’s correlation coefficient (r) and Mean Absolute Error (MAE). The corresponding values of MAE and Pearson’s r are provided in

Table 1. Mean Absolute Error (MAE) and Pearson’s r for annual country level burned areas reported in GFED and modeled in FLAM for the calibration period 2000–2009 and the validation period 2010–2016.

Period	Calibration Period (2000–2009)		Validation Period (2010–2016)
Type of Fire	Peat Fire	Wildfire	Peat Fire	Wildfire
Pearson’s r	0.857	0.865	0.988	0.988
MAE (thousands of ha)	169.25	257.93	35.56	106.71
MAE relative to average burned area	47%	43%	11%	19%

separately for calibration and validation periods. For the validation period, we obtained very high values of Pearson’s $r=0.988$, both for wildfires and peat fires. MAE values were 106.71 and 35.56 thousands hectares, respectively, for wild- and peat fires (19% and 11% relative to the average annual burned area). Lower values of Pearson’s r, 0.857 and 0.865, respectively, were obtained in the calibration period, which can be explained by the extreme fires in 2006 reported by GFED, but not entirely captured by FLAM (see Figure 4 and Figure 5). This could also explain the higher values of MAE in the calibration period compared to the validation period. The ability to capture large fires in 2015 by FLAM in the validation period was high. Here, no information from GFED was used, and the model was running only with updated input data in the period 2010–2016.

In order to compare modeled areas burned by peat fires with the data reported by GFED, we used the peatland mask as described in Section 2.6. The corresponding dynamics of burned area in peatland of Sumatra and Kalimantan is shown in Figure 5. Accumulated burned areas for Figure 4 and Figure 5 can be found in Table 4 below.

3.2. Analysis of the Calibrated Suppression Efficiency Maps

Here, we present detailed results of the calibration procedure (see Section 2.4) that was applied in the validation period. The suppression efficiency maps were obtained using GFED observations, accumulated over 2000–2009, separately for wildfires and peat fires. These maps were used by FLAM for estimation of the spatial and temporal dynamics of burned area in the calibration and validation periods.

3.2.1. Spatial Distribution of Suppression Efficiency for Wildfire

During the calibration procedure, 6% of area in GFED burned in pixels that could not be calibrated in FLAM due to the absence of climate data (see Section 2.4.3). The problem was that the country mask used in the climatic data of GWFED did not coincide with the country mask used by GFED (see Figure A2). For this reason, we excluded those pixels (where climate data were not available) from the analysis. We removed from our consideration about 1% of pixels marked as having any burned area in GFED (totaling about 1% of all burned area), where FLAM was not able to reproduce any fire due to various reasons (e.g., no fuel, no ignition source).

Figure 6 shows how calibrated suppression efficiency varies across Indonesia. Areas with high values of parameter $q_c$ are colored in green. The dark green color is attached to pixels with $q_c=1$, where fire occurrence was excluded by the calibration procedure (see Section 2.4.1). Statistical analysis of ratios of pixels with various thresholds for suppression efficiency, as well as the ratio of accumulated area burned in those pixels are provided in

Table 2. Calibrated suppression efficiency ($q_c$) in the wildfire and its connection to area burned.

Suppression Efficiency Threshold	${\mathit{q}}_{\mathit{c}}<0.1$	${\mathit{q}}_{\mathit{c}}<0.2$	${\mathit{q}}_{\mathit{c}}<0.3$	${\mathit{q}}_{\mathit{c}}<0.4$	${\mathit{q}}_{\mathit{c}}<0.5$	$0.5\le {\mathit{q}}_{\mathit{c}}<1$	${\mathit{q}}_{\mathit{c}}=1$
Ratio of pixels within threshold	0.02	0.07	0.14	0.24	0.35	0.42	0.23
Ratio of aggregated burned area
FLAM = GFED	0.35	0.63	0.8	0.9	0.95	0.05	0
Original fire algorithm ($q=0.5$)	0.01	0.02	0.07	0.13	0.19	0.35	0.46

. We also present the ratio of accumulated burned area in the same pixels for the original fire algorithm without the application of the calibration procedure and with the fixed parameter $q=0.5$ for the entirety of Indonesia. Note that 0.5 was also the smallest value in the study by Arora and Melton [28]. Actually, there was a substantial amount of pixels in diapason $q_c\in [0.5,1)$, i.e., 42% of all pixels. However, only 5% of the total burned area was attributed to these pixels (both in GFED and FLAM) versus 35% in the original algorithm. Note that as FLAM was calibrated using accumulated GFED, their accumulated areas coincided in the calibration period in each pixel, while there was annual variability in burned areas driven by FLAM inputs. The original algorithm models 46% of area burned in pixels, where $q_c=1$, meaning that in fact, according to GFED, there were no fires in those pixels (0% of area burned in FLAM). This illustrates a considerable improvement in spatial accuracy by application of the calibration procedure realized in this study.

Note that according to Table 2, 95% of area burned in pixels where $q_c<0.5$, and 63% of area burned in pixels where $q_c<0.2$ (orange and red areas in Figure 6). This analysis shows that considerable area burned in a small amount of pixels with a very low suppression efficiency. For instance, 35% of area burned in 2% of pixels within threshold $q_c<0.2$. Figure 3 and Figure 6 show that pixels with low suppression efficiency were mostly located in peatland areas. The pixel with the lowest value of $q_c=0.021$ is located in Borneo swamp forest. According to Equation (12), this value corresponds to the average fire duration of 47 days.

3.2.2. Spatial Distribution of Suppression Efficiency for Peat Fire

During the calibration procedure, we found out that 20% of area burned in peatland in GFED was located in pixels where there were no peat fires according to FLAM. This can be explained by the difference in peatland masks used in GFED and FLAM. It seems that GFED allocated peatland to a wider class of land than actual peatland. Our understanding from the publication describing GFED [54] is that the Terrestrial Ecosystems of the World map [56] was used as the source for peatland data. It seems that all types of swamp forests according to the map were considered as peatland. For this reason, we removed those pixels, which were outside of the peatland mask (see Section 2.6) in our analysis, from the GFED data (see Figure A3). We removed from our consideration about 1% of pixels marked as having any burned area in GFED (totaling about 1% of all burned area), where FLAM was not able to reproduce any fire due to various reasons (e.g., no fuel, no ignition source).

Spatial distribution of suppression efficiency in peatland is shown in Figure 7. Again, we outline the pixels, where the probability of extinguishing a fire equals one. The statistics of pixels with various thresholds for $q_c$ in peatland, as well as information for corresponding areas burned in FLAM and the original algorithm are given in

Table 3. Calibrated suppression efficiency ($q_c$) of peat fire and its connection to area burned.

Suppression Efficiency Threshold	${\mathit{q}}_{\mathit{c}}<0.05$	${\mathit{q}}_{\mathit{c}}<0.1$	${\mathit{q}}_{\mathit{c}}<0.15$	${\mathit{q}}_{\mathit{c}}<0.2$	$0.2\le {\mathit{q}}_{\mathit{c}}<1$	${\mathit{q}}_{\mathit{c}}=1$
Ratio of pixels within threshold	0.1	0.28	0.38	0.46	0.27	0.28
Ratio of aggregated burned area
FLAM = GFED	0.36	0.77	0.91	0.94	0.06	0
Original fire algorithm ($q=0.5$)	0.02	0.12	0.25	0.31	0.49	0.2

. Interestingly, 94% of area burned in pixels with very low suppression efficiency, i.e., $q_c<0.2$. The minimum value of $q_c$ was 0.00574, which according to Equation (12), corresponds to the average duration of 173 days. This reflects the actual situation in Southeast Asia, where once established, peat fires may burn for months and are very difficult to control; in particular, when they occur in remote, off-road locations where they are difficult to extinguish using conventional fire-fighting techniques [57].

Identifying pixels with very low suppression efficiency, FLAM is able to capture the phenomena of smoldering fires, the most difficult types of fires to extinguish [58]. Comparing Figure 7 for peat fire and Figure 6 for wildfire, one can see that the areas with low suppression efficiency (red and orange colored areas) often coincide. This observation may indicate possible transition from flaming fires (in forest) to smoldering fires (in peatland). The opposite direction, from smoldering to flaming fires, could also happen in rare cases [58].

3.2.3. Bimodal Density Functions of Suppression Efficiency

To summarize the previous sub-subsection, we have produced histograms and smoothing density functions in R for calibrated suppression efficiency maps. They are presented in Figure 8 for (a) wildfire and (b) peat fire. Both distributions have a bimodal shape, which shows that calibration procedure corrects the areas estimated by the original algorithm in both directions. It removes the areas where the fires could not occur according to observations. Where the observed areas are underestimated, suppression efficiency is adjusted to lower values.

Although mean values of the density functions in Figure 8 (0.64 for wildfire and 0.42 for peat fire) are close to 0.5 (the minimum value in [28]), our analysis above shows that most of the area burned in pixels with very low suppression efficiency. Lower suppression efficiency corresponds to longer fire duration. According to our calibration procedure (see Table 2 and Table 3) and Equation (12), for wildfires, 63% of area burned in pixels where the average duration of fire was more than four days, and in peatland, 77% of area burned in pixels with an average fire duration of more than nine days.

3.3. FLAM vs. Original Fire Algorithm

Finally, we show how the calibration of suppression efficiency improves the agreement between modeled and observed total burned areas, as well as the ratio of peat fire to wildfire in FLAM, as compared to the original fire algorithm.

Table 4. Aggregated burned areas in FLAM, GFED and the original fire algorithm in mln.ha.

Period	Calibration Period (2000–2009)			Validation Period (2010–2016)
Type of Fire	Peat Fire	Wildfire	Ratio	Peat Fire	Wildfire	Ratio
GFED	3.59	5.99	0.59	2.27	3.92	0.58
Original fire algorithm	0.12	3.25	0.04	0.08	1.86	0.04
FLAM	3.59	5.99	0.59	2.27	3.67	0.62

shows the corresponding burned areas for two periods for GFED, FLAM and the original fire algorithm ($q=0.5$ uniformly). In the calibration period, accumulated GFED and FLAM areas coincided due to the calibration procedure (see Section 2.4). Note, that inter-annual dynamics of burned area was not forced to coincide and was driven by FLAM inputs (mostly weather) in the calibration and validation periods (see Figure 4 and Figure 5). The ratio of peat fire to wildfire was 59% in the calibration period. The original algorithm underestimated burned area by wildfire reported in GFED by around 2.74 million ha and very significantly underestimated area burned in peatland (only 0.12 vs. 3.59 mln.ha). This led to an unrealistic ratio of 4% of peat to wildfire areas. This is another illustration of how much the calibration procedure proposed in FLAM improves model performance compared to using a fixed value of $q=0.5$.

Moreover, Table 4 shows that in the validation period, where areas estimated by FLAM and observed by GFED did not necessary coincide, FLAM kept good track of GFED observations. In particular, total areas burned in peatland were almost the same and in wildland were very close. This demonstrates a good temporal fit as reported in previous Section 3.1 for the validation period. Table 4 indicates that the original fire algorithm considerably underestimated burned area in the validation period similar to the calibration period.

3.4. Spatial Validation of Estimated Burned Area

Here, we present the results of the spatial validation of the estimated burned area for the period 2010–2016. We calculated Pearson’s correlation coefficients r based on annual dynamics of burned area (in GFED and FLAM) in each grid cell separately for wildfires and fires in peatland.

3.4.1. Spatial Correlation for Peat Fire

Results in

Table 5. Results of the spatial correlation analysis of estimated area burned by peat fire.

Pearson’s Correlation Coefficient	$\mathit{r}>0$	$\mathit{r}>0.5$	$\mathit{r}>0.6$	$\mathit{r}>0.7$	$\mathit{r}>0.8$	$\mathit{r}>0.9$
Ratio of pixels in correlation map	0.9	0.64	0.57	0.49	0.38	0.29
Ratio of aggregated burned area (FLAM)	0.98	0.84	0.75	0.65	0.48	0.28
Ratio of aggregated burned area (GFED)	0.97	0.82	0.75	0.64	0.5	0.3

show the ratio of pixels with various correlation thresholds, as well as the ratio of aggregated burned area with respect to burned area in all pixels in FLAM and GFED, respectively. One can see that positive values of Pearson’s r were achieved in 90% of pixels, and at the same time, 98% of area was burned in these pixels during the entire time interval (2010–2016). In pixels where 75% of total area was burned in GFED, the correlation is above 0.6. Figure 9 shows the distribution of Pearson’s r across peatland area, and Figure 10 shows the mean annual burned area in each pixel. Comparing those figures, we see that higher correlation coefficients (green-colored pixels in Figure 9) often corresponded to the area where there was high burned area values (purple pixels in Figure 10). At the same time, pixels with a low correlation (orange and red in Figure 9) often coincided with pixels where burned area was relatively small.

3.4.2. Spatial Correlation for Wildfire

Similar to burned area in peatland, we present correlation results for wildfire. As follows from

Table 6. Results of spatial correlation analysis of estimated area burned by wildfires.

Pearson’s Correlation Coefficient	$\mathit{r}>0$	$\mathit{r}>0.5$	$\mathit{r}>0.6$	$\mathit{r}>0.7$	$\mathit{r}>0.8$	$\mathit{r}>0.9$
Ratio of pixels in correlation map	0.83	0.62	0.58	0.53	0.47	0.4
Ratio of aggregated burned area (FLAM)	0.94	0.76	0.67	0.57	0.43	0.24
Ratio of aggregated burned area (GFED)	0.94	0.78	0.7	0.59	0.49	0.28

, in pixels where 76% of area burned, the Pearson’s correlation coefficient was above 0.5. Remarkably in 40% of pixels, the correlation was above 0.9. We present the spatial distribution of Pearson’s r in wildfires in Figure 11 and mean annual burned area in GFED in Figure 12. They show that a high correlation was achieved in areas with large fires, i.e., in Sumatra and Kalimantan, while lower correlation coefficients were often in the areas with relatively small fires.

4. Discussion

In this study, we applied FLAM to the assessment of burned areas in Indonesia. The model was calibrated in the 10-year period (2000–2009) and validated over the seven-year period (2010–2016) using GFED observed burned area. In order to check how results would change with respect to different calibration/validation periods, we considered an alternative split of the historical period into two intervals: seven years for calibration (2000–2006) and 10 years for validation (2007–2016). The corresponding modeling results are given in

Table A1. Mean absolute error (MAE) and Pearson’s r for annual country level burned areas reported in GFED and modeled in FLAM for the calibration period 2000–2006 and the validation period 2007–2016.

Period	Calibration Period (2000–2006)		Validation Period (2007–2016)
Type of Fire	Peat Fire	Wildfire	Peat Fire	Wildfire
Pearson’s r	0.838	0.85	0.914	0.956
MAE (thousands of ha)	155.7	257.51	87.65	154.9
MAE relative to average burned area	46%	43%	38%	36%

. The results show that because of a shorter calibration period, the values of the Pearson’s r are lower compared to those in Table 1. Nevertheless, the correlation with GFED is still very good in the extended validation period (2007–2016): 0.914 for peat fire and 0.956 for wildfire. The corresponding MAE values are 87.65 and 154.9 thousands ha (38% and 36% with respect to average annual burned area). This increase in MAE compared to that of the 10-year calibration period is due to the fact that years 2007–2009, where FLAM is challenged in reproducing burned areas, are now included in the validation period. Experiments with shorter calibration periods provide a decrease in Pearson’s r in validation periods as the model requires sufficient number of years with “peaks” (2002, 2004, 2006) for calibration. In our view, the 10-year calibration period is sufficient for calibrating FLAM.

Calibration of the spatial suppression efficiency performed in this study provides an important development in the methodology of the application of the original Arora and Boer [27] algorithm. Our findings show that the distribution of suppression efficiency has a bimodal density function. A recent publication [28] applying the same core algorithm by Arora and Boer [27] for fire modeling states that in the absence of any human influence, the probability of extinguishing a fire within one day equals 0.5 and corresponds to an average fire duration of one day; the extinguishing probability increases with population density. We show that this global value is unrealistic for tropical forests using the example of Indonesia. The results of model calibration using historical GFED data [29] show that in Indonesia, 95% of area burned by wildfire is located in pixels where the extinguishing probability is below 0.5. It could be as low as 0.0057 in the peatland of the Borneo swamp forest, where the Mega Rice project [59] took place, corresponding to an average fire duration of 173 days. This reflects the actual situation in Southeast Asia, where once established, peat fires may burn for months and are very difficult to extinguish [57]. We show that the application of the algorithm without calibration considerably underestimated burned area, particularly in peatland. Our conclusion is that the approach proposed by Arora and Melton [28] may help to fit a decreasing trend in global burned area, but at the price of considerably underestimating regional burned area, e.g., in Indonesia, which produces a great deal of global CO${}_2$ emissions. Underestimation of burned area by the original fire algorithm with $q=0.5$ could also be the reason why the majority of global fire models underestimate the global burned area, as reported in the study led by Andela [52].

Validation results on a pixel scale show that in pixels where around 70% of area is burned, the correlation coefficients are above 0.6 and in those where around 30% of area is burned, above 0.9, both for peat and wildfire. We show figures for the spatial distribution of correlation coefficients and compare them with the maps of mean annual burned area. The comparison of maps shows that high correlation is often achieved in the pixels with a large burned area, mostly located in the peatland of Sumatra and Kalimantan. At the same time, lower correlation coefficients are often in the pixels with small burned area. Capturing pixels with small fires could be a potential direction for spatial improvement of FLAM modeling, where the GFED version with small fires can be used [53].

We would like to emphasize that the simplified approach used for modeling peat fires produced good results. It is also interesting to mention that areas with large wildfires and peat fires often coincide. Our analysis of the spatial variability of the suppression efficiency for wildfires and peat fires can potentially allow one to find the areas where smoldering and flaming fires are tightly related [58]. However, further analysis has to be carried out in this direction.

FLAM was able to correctly estimate extreme burned area in 2015, which caused a lot of damage in Indonesia [9]. While the calibration procedure uses only accumulated burned area from GFED, the good inter-annual variability is due to adequate representation of daily weather dynamics in FLAM. An interesting observation is that FLAM was not able to capture 2006 in the calibration period. This could indicate an uncertainty in the GFED dataset, bias in climate data or an increase of ignition potential, e.g., due to human activity. Finally, we explored the fire dynamics at a finer temporal scale and showed a comparison of monthly dynamics between GFED and FLAM areas burned by wildfires in Figure A4. Preliminary analysis based on the figure shows that FLAM is not always capturing the highest burned area in September. Potentially, seasonal calibration of the suppression efficiency could lead to the improvement of FLAM monthly estimates, and consequently better spatial and temporal results. A promising direction for future research would be the assessment of greenhouse gases emissions based on burned areas modeled by FLAM.

Current peat fires and wildfires can affect future fires in two main areas. On the non-anthropogenic side, current fires will highly affect soil moisture and the water content of biomass. This will increase the likelihood of fire occurrence within a relatively short period until adequate rainfall compensates the loss of soil moisture and water content. From the anthropogenic perspective, physical alteration of the environment such as canal drainage in peat areas or land clearance for crop production can have an effect. If no restoration activities are carried out, peat fire is likely to reoccur.

5. Conclusions

In this study, we applied FLAM to estimate burned area in Indonesia and validated these results independently of the model calibration for the first time. For burned areas aggregated over Indonesia, we obtained Pearson’s correlation coefficient separately for wildfires and peat fires, which turns out to be equal to 0.988 in both cases. Spatial correlation analysis showed that in pixels where around 70% of area is burned, the correlation coefficients are above 0.6, and in those, where 30% of area is burned, above 0.9. We showed that FLAM is capable of capturing not only flaming, but also smoldering peat fires, which are crucially important for Indonesia.

The suggested approach for the spatial calibration of the probability of extinguishing a fire improves the results of the original Arora and Boer fire algorithm. We showed that the calibrated parameter has a bimodal distribution over Indonesia; in areas where the probability is found to be close to one, the fire occurrence could be diminished due to human activity; at the same time, in areas where this parameter is found to be close to zero, human activity, on the contrary, leads to extensive fire, for example due to slash and burn agriculture or peatland drainage. Very low suppression levels are found in the peatland of Kalimantan and Sumatra, where fire can burn for very long periods of time despite extensive rains, weather changes and fire-fighting attempts.

The results presented in this paper are a prerequisite for estimating future burned area under projected climate change and for modeling adaptation options. Modeling burned area in Indonesia and an assessment of corresponding impacts, e.g., CO${}_2$ emissions, can be potentially used for climate/economic damage assessments. This could also help to adequately represent the risks of wildland fires in REDD+ programs [60,61]. Modeling fire ignition and suppression can also help the implementation of fire prevention policies and provide important information for building adequate and cost-efficient fire response infrastructure.