Drivers and distribution of forest fires in Sikkim Himalaya: A maximum entropy-based approach to spatial modelling.


 The recent episodes of forest fire in Brazil and Australia of 2019 are tragic reminders of the hazards of the forest fire. Globally incidents of forest fire events are in the rise due to human encroachment into wilderness and climate change. Sikkim with a forest cover of more than 47%, suffers seasonal instances of frequent forest fire during the dry winter months. To address this issue, a GIS-aided and MaxEnt machine learning-based forest fire prediction map has been prepared using forest fire inventory database and maps of environmental features. The study indicates that amongst the environmental features, climatic conditions and proximity to roads are the major determinants of the forest fire. Model validation criteria like ROC curve, correlation coefficient and Cohen’s Kappa show a good predictive capability (AUC = 0.95, COR = 0.78, κ = 0.78). The outcomes of this study in the form of a forest fire prediction map can aid the stakeholders of the forest in taking informed mitigation measures.


Introduction
The incidents of forest re in Sikkim Himalaya take a peak during the dry period of the year from November to March due to the accumulation of dry biomass over the forest oor. These incidents may occur by natural causes like lightning as Sikkim falls under the northeast region of India, which is considered a high lightening zone. Anthropogenic causes of forest re in Sikkim include intentional and accidental factors. Bon res by the cattle herders, burning of the forest oor to deter wild animals entering the agrarian land, logging induced decrease in forest canopy cover are the intentional causes of forest re in Sikkim. While, sparks from the uphill moving vehicles, electric transformers located in forested areas, use of traditional torch called Rankoo, throwing away of live bidi and cigarettes butts are the accidental causes of forest re (S. Sharma et al., 2014).
A forest re can be considered as a mixed blessing. Low-intensity forest re opens the canopy cover and removes dead wood, providing a niche for new plants to grow. Also, burning of the forest releases the nutrients bound with the biomass to the soil, rejuvenating forest growth. Forest re also offers new ecological niches for wildlife to proliferate. In contrast, high-intensity forest re leads to loss of soil biota, volatilization of soil nutrients, increase in soil erosion, a decline in biodiversity and forest biomass (Chandra & Bhardwaj, 2015;Parashar & Biswas, 2003).
A wide range of features have been considered for prediction of forest re. According to the review of forest re in the Indian context done by Joseph et al. (2009), topographical features like altitude, aspect, slope, Topographic Wetness Index (TWI) have been used in forest re prediction. These features are partly accountable for the type and distribution of vegetation. Also, topographic features in uence the direction and speed of forest re (Guo et al., 2016;Maingi & Henry, 2007;Minnich & Bahre, 1995). Meteorological features like average precipitation, temperature, humidity and wind speed have been used to understand forest re characteristics. In other studies lightening has been focused to predict forest re (Chen et al., 2015). These meteorological features can act as triggering and conditioning factors for the start and spread of forest re (Guo et al., 2016;Turco et al., 2013). Vegetational features like vegetation type, Normalized Vegetation Difference Index (NDVI), tree cover fraction; human induced features like proximity to road network, human habitation or Wildland-Urban Interface (WUI); and in-situ factors like soil moisture, soil texture and fuel density have also been used for forest re prediction (Gheshlaghi et al., 2020;Jaafari et al., 2018;Mhawej et al., 2015;Satir et al., 2016) A forest re or wild re prediction map has become a valuable tool for disaster management and ecological restoration. Multicriteria decision analysis such as Analytic Hierarchy Process (AHP), Analytical Network Process (ANP) and other forms of expert opinion based methods have been applied in forest re prediction mapping (Gheshlaghi et al., 2020;Goleiji et al., 2017;Ljubomir et al., 2019;Regodic et al., 2018;Yathish et al., 2019). In these methods, the model criteria and alternatives are considered as a hierarchical structure. This is followed by the ranking of the model criteria and alternatives based on a certain scale. Based on the ranking the importance or weights of the model criteria and alternatives are estimated and then used in the GIS-aided prediction mapping (Banerjee et al., 2018). However, expert opinion-based prediction mapping may suffer subjective bias. Moreover, these methods are deterministic. As a result, they may not be suitable for a phenomenon that involves uncertainty, such as a forest re (Ishizaka & Labib, 2009;Mendoza & Martins, 2006). Machine learning methods such as kernel logistic regression, support vector machine, random forest, fuzzy logic, MaxEnt, multilayer perceptron, deep learning and convolutionary neural networks have been extensively used in forest re prediction mapping. Contrary to expert opinion-based methods, machine learning methods do not suffer from subjective bias. Moreover, these methods encompass the uncertainty associated with the modelling of a phenomenon. However, machine learning may suffer issues like model over tting. These methods heavily rely on the training dataset and take time to learn. Furthermore, these methods require a large dataset of events of interest for proper training of the model. Another important limitation of machine learning method, to be speci c methods involving arti cial neural networking, is that, they achieve e ciency and accuracy at the cost of interpretability of the model (Nami et al. 2018;Tien Bui, Le, and Hoang 2018;Ghorbanzadeh, Kamran, and Blaschke 2019;Tehrany et al. 2019;Tien Bui, Hoang, and Samui 2019;Zhang, Wang, and Liu 2019).
Maximum entropy or MaxEnt is a popular machine learning method widely being used in species distribution and earth hazard modelling (Feng & Hong, 2009;Harte, 2011;Pourghasemi & Rossi, 2018). Unlike most machine learning methods such as logistic regression, support vector machine, random forest, k-nearest neighbour and arti cial neural network, that uses presence-absence instances dataset for training, the MaxEnt uses presence-background instances dataset for training. MaxEnt is based on the principle, that the probability distribution that maximizes entropy for the current state of knowledge subject to the constraints of the features is the best t model for the phenomenon under consideration (De Martino & De Martino, 2018). It is popular primarily because it considers 'minimum assumption' while selecting a probability distribution (Warton, 2013). Moreover, this method considers more realistic presencebackground dataset, in the sense that in nature hardly any absence data is available. On the other hand, MaxEnt needs a large presence-dataset to perform reliable prediction. Also, a study has suggested that MaxEnt is equivalent to the Generalized Linear Model (GLM) when it comes to Point Process Models (PPMs) such as forest re events (Fithian & Hastie, 2013). MaxEnt has been used in several forest re prediction mappings. Studies indicate that MaxEnt has performed equally well in comparison to other machine learning methods in predicting a forest re. (Arpaci et al., 2014;Fernández-Manso & Quintano, 2020;Fonseca et al., 2016;Kim et al., 2015;Lim et al., 2017;Massada et al., 2013;Peters et al., 2013) In this study, MaxEnt has been applied to prepare a forest re prediction map of Sikkim Himalaya using MODIS and Ground data-based forest re inventory. As features, meteorological, topological, ecological, in-situ and human-induced data have been used to train the MaxEnt model. Model validation criteria have been used to evaluate the model. The study indicates that MaxEnt is a reliable machine learning method in predicting areas prone to forest re events in Sikkim Himalaya.

Study area
Sikkim is a small eastern Himalayan state of India neighboured by Tibet in the North, Nepal in the West, Bhutan in the East and the state of West Bengal in the south. It extends from 27 • 00′ 46′′ N to 28 • 07′ 48′′ N latitude and 88 • 00′ 58′′ E to 88 • 55′ 25′′ E longitude. The elevation of Sikkim varies from 280 m in the South to 8586 m in the North, crowned by the world's third-highest mountain peak, Mt. Khangchendzonga (Shukla et al., 2018). Sikkim, apart from having four seasons of winter, summer, spring, autumn, has a monsoon season lasting from June to September. It has a sub-tropical climate in the south and tundra climate in the north. The two main rivers of Sikkim include the Teesta River and its tributary, the Rangeet (ENVIS Sikkim, 2019) (Figure 1a Sikkim is endowed with various vegetation ecotypes based on the elevation and climatic conditions, like Himalayan subtropical broadleaf forests in the lower elevations, Eastern Himalayan broadleaf forests in the temperate zone above the elevation of 1500 metres, Eastern Himalayan subalpine conifer forests from 3500 to 5000 metres and Eastern Himalayan alpine shrub and meadows in the higher elevations (O'Neill, 2019;O'Neill et al., 2020). The dry winter season from December to March in Sikkim characterised by windy weather and dry forest biomass, create the right conditions for a forest re. It is more common in the deciduous Sal forest ecosystem followed by the temperate oak and sub-alpine conifer forests. Erratic rainfall, climate change and conversion of forest land into other land uses have increased the vulnerability of the forests in Sikkim, leading to a growing trend of forest re incidents ( Figure 2) (Banerjee et al., 2020;R. Sharma et al., 2012).

Data sources
The active re data from the year 2000-2019 of Moderate Resolution Imaging Spectroradiometer (MODIS) was accessed from the data archive at the Fire Information for Resource Management System (FIRMS) site. The MODIS dataset was combined with the forest re inventory prepared from the GPS-tagged dataset of the Forest and Environment Department, Government of Sikkim. The re incidents dataset thus prepared was intersected with the forest fraction raster data (Shimada et al., 2014) to exclude re incidents beyond the forest cover of the study area. This generated a re dataset of 754 events.
The environmental feature raster maps or simply features used in this study included precipitation (avgPrep), ambient air temperature (avgTemp) and wind speed (avgWind) averaged over the dry period of Sikkim prepared from the monthly data accessed from Worldclim 2 (Fick & Hijmans, 2017). For this study, amongst several data resolutions, the 30 seconds resolution average monthly climate data for 1970-2000 was taken from Worldclim 2. Also, features like aspect, plan curvature (PlanCurv), pro le curvature (pro leCurv), slope and TWI were derived from the Digital Elevation Model (DEM) (Jarvis, A., H.I. Reuter, A., Nelson, E. Guevara. 2008). The DEM used in this study was the product of Shuttle Radar Topography Mission (SRTM) of 90m resolution derived from CGIAR/SRTM90_V4 image with a data collection timeframe from 2000-02-11 to 2000-02-22. The NDVI of the study area was prepared from the Advanced Spaceborne Thermal Emission and Re ection Radiometer data (ASTER Mount Gariwang image, 2018). The NDVI data of 250m resolution was prepared from MODIS/006/MOD13Q1 image, computed from atmospherically corrected bi-directional surface re ectance that has been masked for water, clouds, heavy aerosols, and cloud shadows. The timeframe of NDVI data collection was from 2000-02-18 to 2020-07-02. Features like percent tree cover (Sexton et al., 2013) and population density (CIESIN 2018) were also used in the modelling. The percent tree cover (treeCov) was the product of 30m resolution image of

Maximum entropy -the MaxEnt Model
The MaxEnt algorithm does prediction by minimizing the relative entropy between the probability density of the presence-only instances of the target variable from that of the instances of background landscape data (Elith et al., 2011). For instances, for a landscape, L the algorithm uses forest re occurrence, (y = 1) over a vector of environmental features, z. The MaxEnt algorithm attempts to minimize the distance of the probability densities of features in case of forest re occurrences, from the probability densities of features of the background or the null model, over the landscape L. The minimization of the distance function is achieved by maximizing a penalized likelihood model subject to the model constraints given by the probability densities of features of the landscape (Steven J. Phillips & Dudík, 2008).
Unlike conventional machine learning methods like logistic regression or random forest which uses the presence-absence data, MaxEnt uses presencebackground data for prediction of the forest re. This makes MaxEnt prone to sample selection bias, a condition where some areas of the landscape may be over-sampled than other areas. Also, MaxEnt is prone to over tting the predictive model to the presence-only data (Devisscher et al., 2016). However, recent MaxEnt software uses transformation algorithms like linear, quadratic, product, categorical, threshold and hinge for the standardization of the prediction features and regularization of the model to prevent model over tting (Elith et al., 2011).
Forest res, like most natural phenomena, rarely have data on their absence. This makes MaxEnt more appropriate for forest re studies as it requires presence-only data for their predictions (Arnold et al., 2014;Steven J. Phillips & Elith, 2013).
Data processing and preparation of forest re prediction map Initially, all feature maps were projected from geographic projection system of GCS-WGS-1984 to plane projection system of WGS-1984-UTM-Zone-45N, which is suitable for the study area. Next, Euclidean distance raster maps were prepared from the polyline vector maps of the road network and waterbodies network and point vector map of human habitations. These raster maps were prepared to measure the proximity of re events from roads, waterbodies and human habitations. Topographic feature maps of aspect, slope and TWI were prepared from DEM. Thereafter, all the feature maps were changed to have the same cell size and same extent. Next, the feature maps were normalised, such that the pixel values of the maps were in the range from zero to one (Chang, 2017). The normalized maps were exported in GeoTiff format as they are readily readable by R-programming language. Furthermore, the presence-only dataset of forest re events was stored as a CSV le.
The forest re prediction map was prepared in RStudio environment using R packages named 'raster' (Hijmans, 2020), 'rgdal' (Bivand et al., 2019) and 'dismo' . The rst two packages were mainly used for raster images related spatial operations while dismo was used for bridging between R and MaxEnt software. The MaxEnt software used in this study is a java program-based package (S.J. Phillips et al., n.d.). During the preparation of the prediction map, all the feature maps were stacked with matching extents and feature attributes were extracted from the features stack using the re event coordinates. The re dataset was divided into ve-folds for crossvalidation. This was followed by the preparation of background dataset by selecting 1000 random points from the extent of the study area. Similar to the preparation of re event dataset, the background dataset was populated with the feature attributes and divided into ve-fold datasets for crossvalidation. In this process, repetitively any one sub-dataset out of the ve sub-dataset is used for testing while the rest are used for training the MaxEnt algorithm. An average of the errors generated from the repeated tests help in the tuning of the model parameters for better performance of MaxEnt.
The MaxEnt method-based prediction was applied to the entire extent of the study area, considering every cell as an instant. The model output was exported as a GeoTiff raster and classi ed into ve qualitative categories applying natural breaks method of ArcGIS framework. All Model validation criteria like Receiver Operating Characteristic (ROC) curve (Peres & Cancelliere, 2014), correlation coe cient and Cohen's kappa (Vakhshoori & Zare, 2018) were used to evaluate the model performance. ROC curve is a popular diagnostic and visualization tool that plots sensitivity of the model against speci city of the model.
The correlation coe cient compares how well the predicted values match with the observed values. Cohen's Kappa statistic measures the agreement between two categorical scales, such as binary outcomes of predicted and observed events of forest re. The kappa is the ratio of the deviation of the predicted value from the observed to one less predicted value (Sim & Wright, 2005). In all the three cases, a value close to one is satisfactory for model validation.
Moreover, importance of the features of the model and sensitivity analysis were performed. The methodology of the study is illustrated below (Figure 3).

Results
Correlation analysis indicated that feature like popDen was strongly correlated with avgPrep and avgTemp, primarily because most of the human population was in the south of Sikkim. Also, avgTemp had a strong correlation with avgPrep, as the subtropical Sikkim gets the bulk of rainfall (Table 1). Multicollinearity analysis indicated that Variance In ation Factor (VIF) of the independent feature variables against avgTemp was within acceptable limits (VIF < 5), except for popDen. Hence, multicollinear variable like population density was dropped from the model (Table 2).
Starting with the meteorological features, the re events were more common in moderate to warmer areas (11 -24 o C). Like avgTemp, the re events were more common in areas with moderate to higher avgPrep (35 -55mm). In contrast, re events were common in areas with low avgWind (1.4 -1.7 ms -1 ). Moving onto topographic features, bulk of the re events were observed in the atter slope (5 -7 degree), lower TWI (4 -6 value) and moderate aspect (81 -217 degree). The PlanCurv and pro leCurv had nominal in uence in the prediction of forest res. In case of PlanCurv, convex curvature of 0.765 had some in uence on forest re occurrences. Concave curvature of 0.488 of pro leCurv had some in uence in predicting forest res. Looking at the ecological features, re events were common to areas with moderate to high NDVI value (0.5 -0.7) having moderate treeCov (31-55% of the area). In-site features indicated that areas with moderate soilCarbon and low soilWater were more common to have forest res. Furthermore, re events were skewed towards areas close to the waterbodies (0 -800m), human habitations (0 -3000m) and roadways (0 -800m) (Figure 5, supplement).
In terms of contribution of the features towards prediction of forest re, proxRoads explained 43% of events, followed by avgTemp that explained almost 17% of re events. Environmental features, namely, proxRoads, avgTemp, treeCov, avgPrep, proxPlace and avgWind together explained 85% of the forest re events.
TWI had no contribution in the prediction value. Hence, TWI was dropped from the MaxEnt model ( Figure 6).
The MaxEnt method-based forest re prediction map of 30.7m resolution showed a probability range from 0 to 1, indicating no chances of forest re occurrences to very high chances of forest re occurrences (Figure 7). The prediction map was further categorized into very low, low, medium, high and very high chances of forest re incidents for the sake of convenience (Figure 8). Out of the total study area, 351 km 2 was under very high probability of forest re, while 451 km 2 was under high probability of forest re. Thereby, these two categories together constitute 11% of the total study area.
Model validation criteria like ROC curve showed an Area Under the Curve (AUC) of 0.957 (Figure 9a). The correlation coe cient was estimated to be 0.809 and Cohen's Kappa was estimated to be 0.78 (Figure 9b).

Discussion
In this study an attempt has been made to prepare the forest re prediction map of Sikkim Himalaya using MaxEnt machine learning method. The study indicated that estimation of probability of forest re by MaxEnt was satisfactory as per the model validation criteria.
It was observed that out of the 15 predictor features, six of them, namely proxRoad, avgTemp, treeCov, avgPrep, proxPlace, avgWind, were able to explain 85% of forest re events. Due to the di cult terrain, much of the human activities in Sikkim are con ned around the roadways. This also includes the human habitations. The limited contribution of topographic and ecological features in explaining the forest res further reinforces the disproportionate in uence of human activities and meteorological factors in the occurrences of forest re. Together, these two anthropogenic features explain 47.5% of forest res. This indicates with a fair amount of con dence, that forest res in Sikkim Himalaya primarily occur due to human activities and they are facilitated by climatic conditions. This observation was similar to a previous study done in the Amazonian forest of Bolivia (Devisscher et al., 2016). Also, studies conducted in the Huron-Manistee National Forest, Michigan, USA suggested that development activities such as road networks were the major determinants of forest re occurrences (Massada et al., 2013). Other studies on forest re prediction modelling also suggest that proximity to roads and human settlements are the major determinants of forest res (Arpaci et al., 2014;Jaafari et al., 2018;Maingi & Henry, 2007). Areas with low to moderate tree cover were more susceptible to forest re. Also, areas with moderate temperature and precipitation with very low wind speed were more prone to forest res. Again, these areas are prevalent in the valleys of southern Sikkim Himalaya which are dominated by agrarian land, human habitations and road networks. This observation further indicates that the ecological and meteorological conditions are the conditioning factors of forest re in Sikkim.
In contrast to the mainland India where forest re is common during the hot dry summer (Joseph et al., 2009), majority of forest re in Sikkim Himalaya occurs during the cold dry period from November to March. As observed in this study, bulk of the forest re events were in the southern part of Sikkim. This is mainly due to the logging activities there. Also, higher vehicular tra c explained by greater road network in southern Sikkim makes the dry vegetation vulnerable to re due to engine spark and cigarette or bidi butts. The limited number of forest re events in the high altitude of Sikkim is primarily due to lightening. Moreover, the high contrast of warmer climate in southern Sikkim as compared to very cold climatic conditions of northern Sikkim makes the former more vulnerable to forest re (R. Sharma et al., 2012). This study also indicated that forested areas close to human habitations are at a higher risk of forest re. An aspect from East to South-West direction had more contribution towards forest re. Aspect in uences soil moisture, solar radiation, vegetation composition and density (Estes et al., 2017). Also, forest patches of valley areas that receive moderate rainfall, have moderate temperature and low wind speed were prone to forest re. The higher values of model validation criteria suggested that the model prediction was satisfactory.
Being fundamentally distinct from other machine learning methods, MaxEnt uses presence-only dataset to train itself (Elith et al., 2011). However, like many studies have shown earlier, the present study also indicated that this distinction of MaxEnt does not limit its capability in generating reliable hazard prediction maps (Arpaci et al., 2014;Fernández-Manso & Quintano, 2020;Fonseca et al., 2016;Kim et al., 2015;Lim et al., 2017;Massada et al., 2013;Peters et al., 2013). In this study a limited set of features have been considered for forest re prediction. This was primarily due to availability of reliable data. However, other features like lightening activities in the North Sikkim, dry fuel biomass and vegetation type could be considered to improve the model.
The forest re prediction map of Sikkim Himalaya can be considered as a decision support tool for stakeholders of forest resources. The forest managers such as forest rangers and forest dependent communities can mitigate forest re by allocating their re control resources to areas more prone to forest re.
Population of forest re prone areas can be educated about the impacts of their activities on the occurrence of forest re. Targeted law enforcement against irresponsible activities like illegal logging, negligent smoking and bon re, slash and burn farming and tra c management can be achieved from the forest re prediction map.

Conclusion
Applications of remote sensing imageries, machine learning and geospatial analysis can mitigate forest re by identifying areas that were relatively more prone to forest re. MaxEnt-based forest re prediction map of Sikkim Himalaya indicated that anthropogenic features, mainly road network, tree cover fraction and meteorological features were mainly accountable for forest re incidents. Although a forest re can be considered as an opportunity for the forest to rejuvenate, increase in frequency and extent of forest re can only lead to damage to the forest health. The prediction map can be used as a decision support tool by the stakeholders to mitigate the occurrences forest re. The applications of MaxEnt can be extended to other forms of earth hazards like landslide, ood and drought predictions. The MaxEnt model can be further improved by expanding the feature set, followed by factor analysis to identify the most relevant explanatory features of forest re. A comparative analysis of the MaxEnt along with other machine learning methods can be performed to assess the e ciency and e cacy of MaxEnt. The outcomes of this study can be internalized into forest management policies by applying geographically targeted resource allocation and law enforcement towards forest re mitigation.