Spatial correlates of forest and land fires in Indonesia

Keywords: biomass burning, climate change, fire hotspot, haze, south-east Asia.

Biomass burning, the combustion of organic matter generally in the form of vegetation or the partially decomposed remains of plants that form peat, is a natural, if rare, phenomenon in tropical ecosystems in the absence of human activity. In humid, highly biodiverse rainforests, natural fires are characterised by long return intervals with fairly low, generally reversible impacts on forest stands (Goldammer 2016). Recent decades have, however, seen an increase in the frequency and magnitude of major pantropical fires associated with human activity such as land-use changes, with some of the most severe biomass burnings occurring in peatlands in south-east Asia (Page and Hooijer 2016). Increasing pressure from humans coupled with climate change are likely to drive increases in biomass burning in the future and may lead to alternative stable states in which fires and fire-tolerant vegetation are the norm (Brando et al. 2019).

Peatlands are important terrestrial carbon sinks storing up to 46% of global soil carbon (Page and Hooijer 2016). Degraded peatlands are highly susceptible to burning and can shift from being carbon sinks to sources through oxidation and fire-related emissions; degraded peatlands in the tropics account for ~75% of annual greenhouse gas emissions from modified peatlands globally (Leifeld and Menichetti 2018). Fire management, especially in peatlands, is thus increasingly viewed as a priority for global climate change mitigation.

Indonesia accounts for ~80% of the total peatland area in south-east Asia (Evers et al. 2017; Fig. 1) and >75% of total pantropical wetland (peatland and mangrove combined) climate change mitigation potential (Griscom et al. 2020). Since the late 1980s, large tracts of peat swamp forests in Indonesia have been rapidly converted through state-supported agribusiness schemes into oil palm and pulpwood plantations (Miettinen et al. 2012a; Goldstein 2016). Fire is commonly used as a cost-effective means of preparing land for plantations and agriculture (Simorangkir 2007), and has been resorted to in disputes over land and resources (Dennis et al. 2005).

Fig. 1.  Peatland distribution (highlighted in black) in Sumatra, Kalimantan and Papua, Indonesia. Grey-shaded provinces represent the areas of study. Abbreviated province names are: AC, Aceh; SU, North Sumatra; RI, Riau; SB, West Sumatra; JA, Jambi; BE, Bengkulu; SS, South Sumatra; BB, Bangka Belitung; LA, Lampung; KR, Riau Islands; KS, South Kalimantan; KB, West Kalimantan; KT, Central Kalimantan; KI, East Kalimantan; PB, West Papua; PA, Papua.

Carbon released during biomass burning contributes to atmospheric loadings of greenhouse gases (Parker et al. 2016). Particulate emissions from peatland fires form haze pollution, leading to adverse health impacts (Koplitz et al. 2016; Tan-Soo and Pattanayak 2019) and heavy economic losses (World Bank 2016) in Indonesia and in neighbouring countries. Notable haze events in 1997–98 and 2015 coincided with El Niño, when lower-than-normal rainfall levels generated prolonged drought conditions conducive for burning (Reid et al. 2012; Fanin and Van Der Werf 2017). Haze events were also documented in 2013 and 2014, however, when there were no prolonged droughts in the region. This suggests a decoupling of fires from climate cycles and a possible reduction in overall fire-resistance of the environment (Gaveau et al. 2014).

Statistical modelling of biomass fire distributions as a function of spatially explicit variables related to land use contributes to the understanding of the drivers involved (Stolle et al. 2003). Yet many such studies in Indonesia tend to be restricted to individual provinces and to particular land-cover types (e.g. Stolle et al. 2003; Cattau et al. 2016; Sumarga 2017). Given that biomass fires are widely distributed in Indonesia, and that many regulations regarding fire and land use are established at the national level, the absence of a large-scale investigation of the factors influencing fires that spans multiple, biomass-fire-prone provinces impedes development of targeted management strategies.

To fill this gap, the present paper addresses the following research questions. First, what are the major factors influencing the presence of biomass fires in Indonesia? Second, how have the impacts of these factors changed between El Niño and La Niña years? Our investigation draws on the conceptual framework developed by Stolle et al. (2003), which divides the factors influencing biomass burning into predisposing conditions (biophysical conditions conducive to spreading and sustaining fires) and ignition (actions that have the potential to initiate burning, such as land clearing and arson). For the latter, proxies were used to denote human activities, such as population density and accessibility.

Indonesia is located between 15°S and 8°N and 90°E and 150°E and surrounded by the Pacific and Indian oceans. Our study focuses primarily on provinces on the island of Sumatra, Indonesian Borneo (Kalimantan) and western New Guinea (Papua; Table S1, available as Supplementary material to this paper). All provinces included in this study are peatland-rich and have been the foci of major biomass fires in the recent past.

The country can be divided into three broad climatic regions on the basis of the timing of the main wet seasons: (1) southern Sumatra to Papua, with above-average rainfall from November to March; (2) northern Sumatra, including Riau, and north-western Kalimantan, characterised by heaviest rainfall during October to November and March to May, with intervening drier periods; and (3) Maluku and the northern regions of Sulawesi, with June to July the wettest period (Aldrian and Dwi Susanto 2003; Fanin and Van Der Werf 2017). Seasonal precipitation patterns are also affected by the El Niño–Southern Oscillation (ENSO), which exerts strong influence on fire activities (Fanin and Van Der Werf 2017) from June to November (Aldrian and Dwi Susanto 2003; Spessa et al. 2015), and by interactions between ENSO and the Indian Ocean Dipole (IOD) (Pan et al. 2018). In general, anomalously low rainfall levels are recorded across Indonesia during the northern hemisphere summer and autumn months, i.e. during El Niño development stages, with the opposite case in La Niña years. Notable differences do occur, particularly in the spatial pattern of impacts and the temporal occurrence and frequency of extremes associated with the two phenomena (Supari et al. 2018).

The study utilised data from four years, 2002, 2005, 2011 and 2015, based on a range of climate conditions across the study region and the availability of data in Sumatra, Kalimantan and Papua. The years 2002 and 2015 were climatically dry years associated with the development of moderate to very strong El Niño signals, whereas 2005 and 2011 corresponded to weak and moderate La Niña events respectively and were thus relatively wet (Hendon 2003; Null 2018).

Fire hotspot data were obtained from NASA FIRMS’ (Fire Information for Resource Management System) MCD14ML collection 6 (Giglio et al. 2018). Each fire hotspot represents a 1-km pixel in which a fire or multiple fires have been detected by algorithms used to process data from MODIS (Moderate-Resolution Imaging Spectroradiometer) sensors. The product comes with confidence estimates to indicate the quality of detection. To minimise the chances of false fire hotspot detections, we used data with ≥30% confidence (i.e. nominal to high confidence of true fire hotspot detection; Giglio et al. 2018). A more stringent detection confidence threshold of 80%, which is often reported in Indonesia’s forest and fire monitoring system (Ministry of Environment and Forestry 2015), was also used for a second set of models. Non-fire-related thermal anomalies from volcanic activities were removed from a 5-km buffer around volcanic eruption sites. We did not distinguish between fire ignition and spread (Vayda 2006), as peatland fires can continue to burn belowground, making estimates of the duration of individual fire events difficult.

A mixture of spatial predictor variables was used (Table S2). We used south-east Asia land-cover maps with classes ranging from intact mangrove and montane forests to anthropogenic land-covers, such as mosaics of farmland, regenerating vegetation and urban areas (Miettinen et al. 2012b, 2016). Land-cover classes were reclassified to avoid an overabundance of categories and therefore an insufficient number of data points per class (Sumarga 2017) and to remove any ambiguity over different categories of plantations (Table S3). Additionally, aboveground live woody vegetation biomass (Avitabile et al. 2016) was used as a proxy of potential fuel (Hoscilo et al. 2011).

An averaged fire weather index (FWI) for the whole year and dry season from June to November (JJASON) was taken from the Global Fire Weather Database (GFWED; NASA GISS 2019). Reduced and enhanced rainfall during JJASON are strongly correlated with the progression of El Niño and La Niña respectively (Hendon 2003; Dowdy et al. 2016), and act as a strong predictor of fire activity (Fig. 2) (Spessa et al. 2015; Sze et al. 2019). FWI combines multiple meteorological variables that influence the likelihood of fire ignition and spread, such as temperature, relative humidity, wind speed and precipitation (Field et al. 2015). Given the levels of uncertainties associated with different estimates of precipitation (Heyer et al. 2018), we averaged the Modern-Era Retrospective analysis for Research and Applications, Version 2 (MERRA-2), rain gauge-corrected MERRA-2 and Climate Prediction Centre gauge-based precipitation data. These data were chosen based on availability for all 4 years of analysis. Two separate components of the FWI system – the drought code (DC) and fine fuel moisture code (FFMC) – were also acquired (De Groot et al. 2007; Field et al. 2015).

Fig. 2.  Averaged Fire Weather Index (FWI) based on three gauge-based precipitation data and fire hotspots in the study region from 2000 to 2015. Years with stronger El Niño signals (2002 and 2015) had higher FWI values, corresponding to a higher number of fire hotspots. Years with stronger La Niña signals (2005 and 2011) had lower FWI values and relatively lower number of fire hotspots.

Topographic variables, such as slope and elevation (Jarvis et al. 2008), were used as proxies for the level of difficulty in cultivating crops in an area (Sumarga 2017) and ease of facilitating fire spread on steeper land (Sze et al. 2019). Peatland depth intervals were averaged to obtain mean depth values of 25, 75, 150, 250, 300, 600 and 1000 cm (Wahyunto et al. 2003, 2004, 2006).

Rivers facilitate the access of people to otherwise remote, inaccessible locations and are also important sources of peatland moisture (Sumarga 2017). We included large river networks, expressed as polylines with an upstream catchment area of ~8 km² (Lehner et al. 2008; Lehner 2013). As each fire hotspot denotes the centre of a 1-km pixel with detection of ≥1 fire, the precise locations of the fires are not known; therefore, a distance-based measure was not used to avoid introducing uncertainty. Instead, we calculated river density per pixel using the line density function in ArcGIS (ESRI 2014), generating a raster measuring the number of river segments per pixel.

Most of the data for logging and pulpwood concessions (Greenpeace 2014) were dated according to the year of permit issuance. We assumed that no plantation can be established before receiving permits and removed concessions with no permit issuance date. Thus, our data exclude illegal concessions and those with incomplete information. By comparison, the oil palm concessions dataset (Greenpeace 2014) lacked dated information for most of the permits. Given that the average life cycle of an oil palm plantation is 25–30 years and that 70% of industrial plantations in the region were established since 2000 (Miettinen et al. 2012a), we assumed that the concessions were in existence throughout the modelled years. All concessions were expressed as presence–absence values.

We also investigated the inclusion of land-cover change as a possible driver of biomass fires, but encountered difficulties in arriving at a single variable, owing to the large number of different land-covers and uses and discrepancies among various map data available. We were able, however, to incorporate primary forest cover changes from 2002 to 2012 (Margono et al. 2014) as a proxy of changes in land-cover as a result of human activities (Table S4).

Finally, population density from the Global Population World Grid ver. 4.0 (Center for International Earth Science Information Network (CIESIN) 2018) and accessibility (Nelson 2008; Weiss et al. 2018) were included to reflect the intensity of land use and ignition potential in an area.

All data layers were resampled to a 5-km grid with nearest-neighbour interpolation for categorical variables and bilinear interpolation for continuous variables. A resolution of 5 km was chosen as a trade-off between accuracy and computing power. Continuous variables were standardised (deduct the mean, divide by the standard deviation) to account for different scales in measurement units across variables. In an initial set of analyses, pixels with fire hotspots (hereafter referred to as ‘fire pixels’) were assigned a value of 1 to denote presence. The same number of pseudo-absence pixels as presence pixels was randomly selected to avoid the problem of unbalanced sampling (Sumarga 2017).

We dropped the urban and water classes from the land-cover data because these classes are unlikely to be associated with biomass fires. Predictor variables were tested for multicollinearity and dropped if the general variance inflation factor was greater than 4 (Zuur et al. 2009). Through this process, we dropped the annual FWI data and separated the FWI subcomponents into different analyses. For the analyses using 80% hotspot detection confidence, the mangrove class was combined with the forest class owing to an insufficient number of fire pixels.

Generalised linear mixed-effects models (GLMMs) with a binomial error structure, a logit-transformation and provinces or subregion as random effect(s) were fitted to the data in the first set of analyses. Regional differences in fire occurrences were expected but were not the main interest of our models and therefore included as random effects. A model was fitted for each year to capture years with El Niño and La Niña events. The models followed the general equation:

where: Y_ij = 1 if pixel j in province i has fire and 0 otherwise; p = probability of a fire pixel; α = intercept; β_n = coefficients of predictor variables; x_n = predictor variables; a_i = random intercept, which is assumed to be normally distributed with mean = 0 and variance = .

Model results were expressed as the odds ratio (OR) between the probability of a pixel being a fire pixel against the probability that it is not (Stolle et al. 2003). The OR was calculated by taking the exponent of the coefficient estimate (Peng et al. 2002), where OR > 1 indicated a positive effect and 0 < OR <1 a reduced effect of the predictor variable on the response variable. For categorical variables, OR was the ratio between the odds of a class compared with a baseline class: concession variables compared with the absence of the concessions; mangrove, mosaic, open area, water and plantation covers compared with the forest cover; and forest degradation status with cleared forest (Stolle et al. 2003). We then conducted model simplification through backwards stepwise selection to remove non-significant predictor variables (P > 0.05) based on chi-square tests from an initial global model that contained all variables (Crawley 2013). This resulted in different models reported for each year.

In another set of analyses, the response variable was fire hotspot count per pixel. Count data are often used as a proxy for the extent of actual fires in which the number of fire hotspots is correlated with the average size of area burned (Tansey et al. 2008). We used zero-inflated negative binomial (ZINB) regression to account for overdispersion and excess zero counts arising from underlying processes independent of those generating count values. Where the zero-inflation component was insignificant, we used a negative binomial model (NB) instead. Two NB calculation methods are available for modelling linear and quadratic mean-variance relationships. We selected the appropriate model based on the greatest reduced dispersion ratio and Akaike Information Criterion (AIC) (Bolker et al. 2012). Model coefficients of >0 indicated a positive correlation and <0 a negative correlation between the predictor variable and response variable.

All data were processed in ArcGIS 10.3 (ESRI 2014). Analyses were performed in Rstudio using the lme4 (Bates et al. 2015) and glmmTMB (Brooks et al. 2017) packages.

The fit of the GLMMs was partially assessed by the R² values using the MuMIn package (Bartoń 2018). The marginal R² refers to the percentage of variance explained by the fixed effects, whereas the conditional R² refers to the variance explained by the overall model. The package further reports the theoretical and delta R² for binomial distributions. The former assumes a theoretical variance of π²/3 for all observation-level data, whereas the latter uses the observation-level variance (Nakagawa et al. 2017). Additionally, we calculated the Brier Score and Area under the Receiver Operating Characteristic Curve (AUC) using the ROCR package (Sing et al. 2005). The AUC measures model accuracy where an AUC = 0.5 amounts to random model prediction and AUC = 1 indicates an accurate model (Elith et al. 2011). The Brier Score measures the models’ predictive accuracy along a range of 0–1, where a score of 0 refers to high predictive accuracy and vice versa (Brier 1950). For the ZINB models, the mean squared error (MSE), root-mean-square error (RMSE) and mean absolute error (MAE) were calculated using the Metrics package (Hamner and Frasco 2018).

Fire hotspots were distributed unequally across the study area, with most concentrated in a few provinces, such as Central Kalimantan and South Sumatra. Although the percentage of fires in the two provinces in Papua constituted only a small fraction of the total fire hotspots detected, the number of hotspots in 2015 (14 793 hotspots) greatly exceeded those in 2002 (5465 hotspots). Fires also tended to cluster on peatlands compared with sites on mineral soils (Table 1), especially in areas where there are plantations and other forms of human activities (Table S5).

Table 1.  Fire hotspot density between non-peat and peat soils

Fire hotspots were 2–3 times more numerous in the 2 dry years 2002 (n = 91 183) and 2015 (n = 143 111), compared with 2005 (n = 60 821) and 2011 (n = 50 196) with wetter weathers. The sensitivity to meteorological effects varied by provinces across the study region. For example, Riau, on the island of Sumatra, experienced a high number of fire hotspots from January to March, irrespective of ENSO strengths. This suggests either precipitation patterns in Riau are not closely regulated by ENSO, or a decoupling of burning from rainfall levels in the province. In terms of temporal distribution, fire hotspots generally started appearing in July and peaked in August to September. A smaller peak in fire hotspots during February was attributed to several provinces in northern Sumatra and north-west Kalimantan, which have two dry seasons. By comparison, fire hotspots in the two provinces of Papua lagged by a month, peaking in September to October (Fig. S1).

The first set of analyses had 31 206, 23 856, 24 376 and 34 424 pixels for the years 2002, 2005, 2011 and 2015 respectively. Provinces in south-west Sumatra, West and Central Kalimantan, and south Papua evinced higher probabilities of fire (Fig. 3). Burning tended to occur in low-lying areas that were fairly accessible, but also supported low population density (Fig. 4). Areas with mixed-production systems (i.e. mosaic and open areas) where forest cover had been replaced by agriculture or shrub vegetation were most prone to fires (OR 1.41–2.86 compared with forested areas). Plantation and regrowth land-cover were also significant predictors of fires (OR 1.70–2.16), although the effects depended on the type of concessions. Of the three types of concessions considered, odds of fire increased within pulpwood and oil palm and decreased within logging. Where peatlands were converted to agriculture, susceptibility to burning increased with peat depth (OR 1.05–1.4 in mosaic, open-area and plantation covers on peatlands). Primary intact forests were less likely to burn compared with non-primary forests and degraded primary forests. Relatively high FWI values were associated with greater probability of fire pixels that characterised the 2 dry years (2002 and 2015) included in the analysis.

Fig. 3.  Predicted fire probabilities (scaled from 0 to 1) for sampled pixels in the study area based on the generalised linear mixed-effects models (GLMMs) with 30% fire hotspot detection confidence. Higher fire probability (dark red) was consistently found in several provinces such as South Sumatra, West and Central Kalimantan, and the southern tip of Papua in the 4 study years. A colour version of this figure is available online.

Fig. 4.  Odds ratios (ORs) of predictor variables extracted from the generalised linear mixed-effects models (GLMMs) for the 4 years of analysis. Predictor variables with OR > 1 have a positive effect on increasing fire occurrences and those with OR < 1 indicate a negative effect. The predictor variables showed similar ORs for all 4 years.

Approximately 40–42% of the variance was explained by the models, 29–33% of which was explained by the predictor variables. Model fits were acceptable with AUCs of 0.83–0.86 and Brier Scores of 0.15–0.16 (Table S6).

Setting 80% as the detection confidence removed ~50% of pixels from our existing data, resulting in sample sizes of 16 320, 10 442, 10 390 and 18 762 for respectively 2002, 2005, 2011 and 2015. The overall effects of the predictor variables in terms of direction and magnitude remained consistent to the models with ≥30% fire hotspot detection confidence and are shown in Fig. S2 and S3 rather than described here. A notable difference is the omission of logging and oil palm concessions from the 2015 model, oil palm concessions and biomass from the 2011 model, and pulpwood concessions from the 2005 model because they were not significant.

The results of the ZINB models were largely consistent with the results of the GLMMs. There were 52 216, 52 216, 51 947 and 52 719 pixels for the models of the years 2002, 2005, 2011 and 2015 respectively. The presence of pulpwood concessions, disturbed land-covers and degraded peatlands was positively correlated with the occurrence of higher numbers of fire hotspots (Fig. 5). For intact peatlands (not developed into plantations or modified in other ways), peat depth was negatively correlated with the number of hotspots. Oil palm concessions were associated with an increase in fires in 2002 (β_OilPalm = 0.11, 95% confidence interval = [0.07, 0.15]) and 2005 (β_OilPalm = 0.44 [0.38, 0.51]), but not in 2011 (β_OilPalm = −0.05 [−0.1, 0]) and 2015 (β_OilPalm = −0.06 [−0.1, −0.02]). In contrast, lower number of fires were associated with higher river density, biomass, population density and elevation. Primary intact forests showed a marked decline in fire count compared with primary degraded or non-primary forests. Model evaluation metrics are reported in Table S7.

Fig. 5.  Non-transformed coefficients of the negative binomial and logit components for the zero-inflated negative binomial (ZINB) models using fire hotspot count as the response variable. Predictor variables with coefficient > 0 were correlated with higher number of fire hotspots. The zero-inflated intercept represents the odds of encountering a structural zero.

Several provinces such as West Kalimantan, Central Kalimantan and South Sumatra had relatively large numbers of fires (Fig. S4) due to unaccounted for subregional and province-level factors. Lampung, Aceh and West Papua had comparatively low numbers of fires across all years.

Land-covers characterised by mixed production systems and plantations, in addition to the degree of access, exhibited strong influence over the occurrence of fires. Our models indicate that fires occur in response to anthropogenic land use and cover and are thus in agreement with the findings of Vetrita and Cochrane (2020) based on MODIS satellite data for Sumatra and Kalimantan over the last two decades. However, our results were novel in distinguishing two different levels of detection confidence. We found that predictor variables generally had the same effects on fires across the years, irrespective of the occurrence of drought conditions, and whether fire hotspots were detected using the 30 or 80% detection confidence. This finding has important implications for evaluating the performance of fire interventions, usually reported as the reduction in hotspot count. Only comparing the absolute fire hotspot numbers pre- and post-intervention overlooks how different detection confidence thresholds, oftentimes underreported, can affect hotspot counts, and subsequently, the fire pixels. Clearer reporting of the detection confidence, along with comparing the relative effects of predictor variables on fire occurrence across larger geographical regions, can yield a clearer picture of intervention effectiveness.

Although Papua and West Papua provinces only account for a small percentage of total fire hotspots, the increase in number of fire hotspots between 2002 and 2015 may reflect a rapid clearing of primary forests (~0.2 Mha loss from 2002 to 2012) that rendered the landscape increasingly flammable (Margono et al. 2014). These results highlight the importance of not only prioritising areas of high fire probability for fire management based on specific land use and predisposing conditions, but also serve to focus preventative interventions on regions such as Papua with rising fire occurrences in more recent years.

Our results corroborated previous work identifying the importance of anthropogenic drivers in biomass fires in Indonesia. However, we found that high population density lowered the number and probability of fire pixels and hotspots, a result supported by Sze et al. (2019) and Lilleskov et al. (2019) in their more finely resolved, province-specific analyses. The relationship is reflective of fire-setting activities being initiated on land that is fairly remote from settlements and therefore more likely to be less closely monitored and regulated. Where population densities are high, such as in and close to urban areas, fire occurrences are low owing to the shortage of combustible biomass and higher levels of fire suppression (Syphard et al. 2009; Price and Bradstock 2014; Cattau et al. 2016). Under such conditions, identifying perpetrators is difficult (Thung 2018). Interventions targeted at building local capacity and institutions to manage uncontrolled fires can be important in reducing the overall fire risk in these areas. Examples include training and financing firefighters, such as local firefighting communities established by the Ministry of Environment and Forestry (Ni’mah et al. 2018), and providing incentives for alternative zero-burning practices.

The effects of anthropogenic activities on combustion increase with peat depth. In contrast, organic-rich soils supporting forest and mangrove vegetation had lower fire risks due to the presence of relatively high watertables. The results support regulations such as the Indonesian Government’s moratorium restricting development on deep peatlands. Where peatlands have been degraded, programs to restore watertable levels and regenerate native vegetation can be key to reducing fires and transboundary haze. In 2016, the Indonesian Government launched a National Peatland Restoration Agency (Badan Restorasi Gambut, BRG) with the goal of restoring 2.6 Mha of peatlands, two-thirds of which are in corporate concession areas (Astuti 2020) within a 5-year-long period. To improve fire prevention, BRG established a system for monitoring peatland watertable levels across Indonesia in anticipation of drought conditions (Arumingtyas 2019). Launch of the BRG constitutes a first, ambitious step in combating peatland degradation, although further involvement of plantation owners and smallholders in providing effective support for efforts to manage fires on peatlands in the long term is crucial.

We expected lower probability and number of fires within logging concessions, as the initial stages of logging operations require little to no land preparation involving fire use. Using similar datasets, Abood et al. (2015) showed that forest cover is higher within logging concessions compared with pulpwood and oil palm concessions. Leftover biomass litter from logging to clear land for monoculture concessions can provide fuel for subsequent burning, increasing the risk of repeated fires at later stages of development (Hoscilo et al. 2011).

Oil palm concessions were not a significant predictor of fires in 2015. Fires are frequently employed during the initial stages of land preparation for planting oil palm, and fire risks decline as the crop matures (Gutiérrez-Vélez et al. 2014). A closer examination of the dataset revealed that dated permits were mostly granted in the late 1990s–early 2000s, corresponding to the period of growth in privately owned oil palm plantations (Goldstein 2016). By 2015, many of these concessions – if established – would have reached production maturity. However, the development of palm plantations is a multistage and multi-actor process, with establishment often triggering an influx of migrants to surrounding areas, especially when smallholders are subcontracted by plantation companies to grow oil palm on land they have access to. This pattern of land use is consistent with data showing rises in smallholder participation in oil palm production (Sheil et al. 2009) and aggregation of smallholder agriculture around plantations (Sloan et al. 2017). Our models did not factor in distance-based predictors for concessions, which can fail to detect the interactions between concessions and surrounding land use. Longer-term studies accounting for distance-related measurements from concessions could reveal more insights into the agents and drivers involved in burning.

FWI was also a major predictor of fires in our study, especially for 2015 with stronger El Niño signals. Its subcomponents, the FFMC and DC, are indicative of probability of ignition and potential to burn respectively (De Groot et al. 2007), and are affected strongly by El Niño type and IOD phase meteorological conditions (Pan et al. 2018). We found all these indices to have greater effect on fire probability during the years with stronger El Niño signals (Figs S5 and S6). However, our models used a fairly coarse averaged-value approach and therefore only captured parts of the seasonal and geographical variations in meteorological patterns and fires. The literature suggests that relationships between meteorological factors and fires may be non-linear and based on threshold values (Spessa et al. 2015; Field et al. 2016). For example, 200–305 mm of rainfall accumulation is a crucial threshold for determining fire susceptibility in southern Sumatra and Kalimantan (Fanin and Van Der Werf 2017). El Niño episodes are projected to become more frequent and intense in the region in the future (Cai et al. 2014), thus increasing the frequency of meteorological conditions conducive to major biomass fires. Early warning systems that enable effective fire management interventions, based on accurate predictions of the sensitivity of peatland biomass to burning, will therefore become an increasing part of attempts to mitigate climate change and protect important forest and peatland ecosystems and their services, including carbon storage.

The current study has some limitations related to data and model choice. Owing to data paucity, we excluded crucial predictors that are known to be strong predictors of fires on peatlands, such as canal infrastructure and land-use zonings (Hoscilo et al. 2011), potentially reducing the R² as a result. With regards to the land-cover data, improvements were made in the 2015 mapping approach to better classify land-cover classes. As such, the 2015 land-cover map was not directly comparable with the previous years, resulting in only 4 years selected for modelling, which may bias against years with less ENSO influence. Resampling was conducted owing to different spatial resolutions in the variables, which are a potential source of error. Occasionally spatial errors can be as high as 50% of the pixel in the case of categorical nearest neighbour resampling, such as the reduction in land-cover heterogeneity, especially when resampling from finer to coarser resolutions (Christman and Rogan 2012). For continuous surfaces, resampling may result in loss of accuracy that can affect model results, such as generating the opposite estimation directions (Dixon and Earls 2009). However, resampling was necessary to ensure that datasets were comparable for analysis.

A further limitation is that only a cross-section of fire events was modelled, whereas many of the processes and conditions related to the predictor variables have delayed effects. For example, fire impacts may lag behind land-use changes because of the time taken for fuel loads to accumulate beyond a threshold conducive for fires to start and subsequently spread. Our models were also unable to test for the effects of climate patterns on fire occurrence, given that we only conducted analyses for 4 years with El Niño and La Niña signals. Incorporating temporal lags and time series analyses into future models may reveal more nuanced impacts of the predictor variables.

Land-use and land-cover factors were persistent determinants of fire hotspots across Indonesia over the 4 years of study. This highlights the importance of preserving forests and peatlands in relatively pristine states irrespective of climate conditions. Although ENSO and IOD events can alter the susceptibility to fires, the prominence of human-related drivers indicates the importance of land interventions and appropriate policies in reducing landscape flammability. Given the spread of fires into areas where there are deep accumulations of peat, interventions to restore and protect peatlands are likely to assist recovery and reduce the risk of major biomass burning events. Current moratoria on modifying peatlands for plantation agriculture are useful in halting concession expansion but are unlikely to reduce flammability in already-degraded sites or the use of fire in particular settings. Regulations influencing how peatlands are used should be accompanied by projects aimed at restoring and rehabilitating degraded ecosystems and promoting alternative livelihoods less reliant on drainage and fire. In this context, BRG holds promise for improving future resilience to fires by rewetting and revegetating peatlands, so long as the program can be sustained over the medium to long term. Additional assistance to smallholders in the forms of technical expertise and financial assistance can provide alternative means of land preparation based on zero burning. In addition to contributing to the conservation of biodiversity, reducing the risk of biomass fires will be a critical part of ensuring the long-term security of important carbon stores, in the form of tropical forests and peatlands, and thus mitigating climate change through management of the global carbon budget.

The fire hotspot data were obtained from NASA FIRMS (https://earthdata.nasa.gov/earth-observation-data/near-real-time/firms/active-fire-data), aboveground woody biomass data from GEOCARBON (https://www.bgc-jena.mpg.de/geodb/projects/Home.php), south-east Asia land cover data from Centre for Remote Imaging, Sensing and Processing (https://ormt-crisp.nus.edu.sg/ormt/Home/Disclaimer), digital elevation data from Consultative Group for International Agricultural Research (http://srtm.csi.cgiar.org/srtmdata/), Hydrological data and maps based on SHuttle Elevation Derivatives at multiple Scales (HydroSHEDS) river data from USGS (https://hydrosheds.cr.usgs.gov/datadownload.php?reqdata=15rivs), fire weather index data from the global fire weather database (https://data.giss.nasa.gov/impacts/gfwed/), accessibility data from the Malaria Atlas Project (https://malariaatlas.org/research-project/accessibility_to_cities/) and Joint Research Centre of the European Commission (https://forobs.jrc.ec.europa.eu/products/gam/download.php), concession data from Greenpeace (https://www.greenpeace.org/archive-indonesia/Global/seasia/Indonesia/Code/Forest-Map/en/data.html), population data from Socioeconomic Data and Applications Center (https://sedac.ciesin.columbia.edu/data/set/gpw-v4-population-density-adjusted-to-2015-unwpp-country-totals-rev11) and forest cover loss data from Global Land Analysis and Discovery (https://glad.umd.edu/dataset/primary-forest-cover-loss-indonesia-2000-2012). Peatland distribution data can be purchased from Wetlands International (http://www.wetlands.or.id/publications _maps.php) at a low cost. Most data accessed 25 November 2019, fire weather index data accessed 3 March 2020.

The authors declare no conflicts of interest.