Forecasting wildfires in major forest types of India

Introduction

Numerous studies have revealed that forest fires have become increasingly severe (Reilly et al., 2017; Singleton et al., 2019). Several incidences of catastrophic wildfires have been reported in the last 4–5 years for example in California, USA (Du et al., 2019) and Australia (Bowman et al., 2020). In the Indian context, an unprecedented wildfire occurred in the state of Uttarakhand in 2016 (Sati and Juyal, 2016), which burnt 2,166 km² within and outside the forest area (Jha et al., 2016). As the forests are habitats for the tribal/natural forest dwellers and the storehouse of floral and faunal biodiversity, it is important to safeguard them from a devastating wildfire. It is difficult and costly to contain the wildfire after it has spread significantly, thus mid- to long-duration forecasts (monthly and yearly) of the potential forest fires are important to develop well-informed mitigation strategies.

Some of the previous studies were carried out to forecast the number of forest fires using the ARIMA (Autoregressive Integrated Moving Average) model (Slavia et al., 2019; Kadir et al., 2020); estimate the probability of fire on a defined day using Probit model (Albertson et al., 2009); model the number of forest fires using longitudinal negative bionomial (NB) and zero-inflated negative bionomial (ZNIB) mixed models (Viedma et al., 2018); predict fire activity for 1–5 days in the future using MODIS satellite data active fire count and ERA-interim reanalysis-based weather data using Poisson’s regression method (Graff et al., 2020); and map the fire risk using machine learning and satellite vegetation index time series (Michael et al., 2021) and an analytical hierarchical process (Sivrikaya and Kucuk, 2021).

The impact of different drivers on the occurrence of forest fires has also been studied in numerous studies (Westerling et al., 2002; Rodrigues et al., 2014; Biswas et al., 2015; Kale et al., 2017; Viedma et al., 2018). Extreme climate events, i.e., ENSO (El Niño Southern Oscillation), NAO (North Atlantic Oscillation), and AMO (Atlantic Multi-decadal Oscillation) are found to be associated with forest fires (Siegert and Hoffmann, 2000; Siegert et al., 2001; Patra et al., 2005; Kidzberger et al., 2007; N’Datchoh et al., 2015; Devischer et al., 2016). Kale et al. (2017) found an interrelation of NINO 3 index with a number of forest fires in India. The drivers such as fuel type (Bond and Keeley, 2005), population/human influence (Bowman et al., 2011), topography (Rollins et al., 2002), distance to settlement and roads (Laurance et al., 2009), aspect (Allexander et al., 2006), rainfall, and temperature (McKenzie et al., 2004; Aldersley et al., 2011) have been identified to be related to forest fires. Urbieta et al. (2019) emphasized that, despite the increasing fire risk factors, fire can be controlled and reduced provided fire suppression resources (mainly aerial) are increased. This emphasizes the need to develop integrated early warning and response systems.

Efforts have been made world-over to develop fire danger rating systems/alert systems/information systems to disseminate timely information about potential fire danger to different stakeholders. The national fire danger rating system (NFDRS) of the United States of America uses fuel type, topography, and weather as inputs to estimate fire danger for next 24 h (USDA, 2022). The Canadian Forest Fire Weather Index (FWI) system provides maps of fire danger and fire weather indices on a daily basis (Natural Resource Canada, 2022). The European Forest Fire Information System (EFFIS) provides services related to wildland fire in Europe. EFFIS provides “fire danger” forecast upto 9 days in the European region (EFFIS, 2022). Forest Survey of India uses satellite data [SNPP (Suomi-National Polar-orbiting Partnership)-VIIRS] to identify and track large forest fires on a near real-time basis in India (FSI, 2022). In addition, the sophisticated numerical two-way coupled models such as WRF-FIRE have been developed and tested for operational fire spread simulation in different countries, i.e., Israel (Mandel et al., 2014), Greece (Giannaros et al., 2020), and USA (UCAR, 2022) to quickly (hourly basis) predict the wildfire spread once the ignition location is detected.

Even when the different fire early warning systems are available to provide the quickest possible warnings about the wildfire location, severity, and danger, often, it is difficult to deploy the resources quickly to contain the wildfire at short notice. Response time is further delayed when the terrain is complex. It is therefore the mid- to long-term (monthly and yearly) fire forecasts that are the key for pre-allocating resources at potentially vulnerable locations during the fire season, hence shortening the response time of the forest fire fighting authorities.

Time series analysis (TSA) of the forest fire is useful in forecasting fire events at specified time intervals. Forecasts may be based on multiple independent (causal) variables or simply on the fire time series itself (univariate). Fire forecasts considering independent variables are based on dynamic interrelationships of wildfires and associated independent variables with respect to time. Univariate and multivariate times series models have been tested and compared previously in numerous cross-disciplinary studies, for example, for estimating tourist demand (Preez and Witt, 2003) and forecasting the medical emergency department demand (Jones et al., 2009; Sarfo et al., 2015). Some studies have found the multivariate time series models superior to univariate models (Jones et al., 2009; Sarfo et al., 2015), whereas, others have recommended univariate models over the multivariate models, for example, Iwok and Okpe (2016), recommended univariate model while studying different time series variables from Nigeria’s gross domestic products. Sethi and Mittal (2020) found the univariate ARIMA model superior to the multivariate vector autoregression model (VAR) for the prediction of air quality index. TSA is a data-intensive process and requires several years of data to derive meaningful interpretations. As compared to multivariate model, the univariate model is less data demanding.

Modern remote sensing platforms (along with ground sensor networks) have made it possible to obtain such datasets at higher temporal and spatial resolutions. For example, NASA’S Fire Information for Resource Management System distributes fire point data on daily basis within 3 h of observations by MODIS and VIIRS satellite sensors at 1 km and 375 m spatial resolution, respectively (EARTHDATA, 2022). Similarly, the network of different ground meteorological stations world over has made it possible to obtain the high-frequency data, for example, of rainfall and temperature.

There are different methods in vogue that use the time series data for forecasting. For example, linear and multiple regressions, exponential smoothing, autoregressive Integrated moving average (ARIMA), logistic regression, and Probit models have been frequently used in carrying out fire predictions (Contreras et al., 2003; Preez and Witt, 2003; Preisler and Westerling, 2007; Prestemon et al., 2012; Briët et al., 2013; Huesca et al., 2014; Zhang et al., 2016; Lei, 2017; Ye et al., 2017; Feng and Shi, 2018; Santana et al., 2018; Slavia et al., 2019). The TSA involves analysis of data collected over a specified time interval, i.e., day, month, week, or year and is useful in near future forecasting. Among different time series forecasting models, the ARIMA model has been used extensively in cross-disciplinary domains (Mujumdar and Kumar, 1990; Wangdi et al., 2010; Sarfo et al., 2015; Zhou et al., 2018; Jesus et al., 2022).

The forecasts obtained through TSA models depends upon how the variables are defined. In the present research the “active fire counts” has been used as a dependent variable, thus, the forecast is “active fire count.” The active fire count data is useful in estimating fire susceptibility using the frequency ratio method (Lee and Pradhan, 2007; Pradhan, 2010; Biswas et al., 2015; Kale et al., 2017).

The present research is aimed at (i) forecasting active forest fire counts using the ARIMA model (univariate ARIMA and ARIMA model with regressors) in prominent forests of India and (ii) deriving probabilistic fire susceptibility map in GIS using fire frequency ratio (FFR) method.

Study area

The study has been carried out in Tropical Evergreen forests (TEGF), Tropical Moist Deciduous forests (TMDF), and Tropical Dry Deciduous forests (TDDF) located in different biogeographic zones of India (Figure 1). The Indian region can be divided into ten bio-geographic zones, i.e., Trans-Himalaya, Himalaya, Desert, Semi-Arid, Western Ghats, Deccan Peninsula, Gangetic Plain, Coasts, North-East, and the Islands (Rodgers and Panwar, 1988). High elevation regions of the Himalayas receive significant snowfall. Parts of western and central regions receive meager rainfall, whereas, eastern India is one of the highest monsoon rainfall regions of the world. The type of the forest is attributed mainly to climate, soil, and past treatment (Champion and Seth, 2005). The TEGF, TDDF, and TMDF are distributed in different proportions in the Western Ghats, Deccan Peninsula, Gangetic Plain, Himalaya, and North-East regions. The temperature and precipitation in these forests vary between 8 and 38^°C and 400 and 7,000 mm, respectively (Roy et al., 2015).

The fires in TMDF are generally sporadic and patchy (Kodandpani et al., 2008). These forests become prone to fire when disturbed. Once disturbed and fragmented, the grasses become prevalent which makes these forests susceptible to fire (Woods, 1989; Suresh et al., 1996; Freifelder et al., 1998; Kodandpani et al., 2008). Fire leads to a further spread of grasses which becomes heavy and continuous, particularly where there is an open canopy (Champion and Seth, 2005). The TMDF supports significant biodiversity as they are comparatively denser and moister than TDDF. It is therefore extremely important to understand the possible future trends of forest fires in these critical forests so as to take mitigative measures. TDDF are significantly burning forests of India. These forests have associated grasses (particularly in central India) and exhibit seasonality. Leaves are shed in specific seasons, which usually start with the onset of summers. During extreme climate, due to increased temperature and dryness, significant burning happens in these forests. TEGF are comparatively lesser burnt forests due to the availability of moisture even during the fire season. These forests are confined to narrow strips in Western Ghats, North-East India, and Andaman and Nicobar Islands (Roy et al., 2015).

Forest fire generally follows a well-defined timeline in India. The fire season in India is mostly from February to June (Satendra and Kaushik, 2014; Kale et al., 2017). Every year wildfire peaks in the month of March/April and then decline at the start of the rainy season (June). Such seasonality continues year after year. However, geographical differences in the duration and peak of the season do exist, requiring in-depth analyses such as the one presented here.

Materials and methods

Different datasets were used in this research for forecasting the active forest fire counts [using univariate ARIMA model and ARIMA model with regressors (NAU, 2020)] and estimating FFR.

The dependent variable, i.e., active fire locations (2003–2017) was sourced from MODIS level 6 active fire product (MODIS, 2022). This depicted fires burning in 1 km pixel at the time of satellite overpass. Independent variables, i.e., monthly dry days (days without rainfall) and maximum average monthly temperature were estimated based on the daily rainfall grid (25 km × 25 km) and temperature grid (1° × 1°) from 2003 to 2017 obtained from India Meteorological Department (IMD; Srivastava et al., 2009; Pai et al., 2014).

Other datasets, i.e., the vegetation type map of India (scale 1:50,000) providing details of natural and semi-natural land-use and land-cover systems and the natural vegetation (Roy et al., 2015) were used to extract the TDDF, TMDF, and TEGF of India; elevation (90 m resolution) was sourced from SRTM¹ and used to derive the slope map. These datasets were used as inputs for FFR estimation in addition to temperature and dry days. Biogeographic zones of India (Rodgers and Panwar, 1988; Rodgers et al., 2000) were used for zone-wise segregation of active fire and FFR.

R software package was used to run the ARIMA model (primary modules used were Forecast (Hyndman and Khandakar, 2008; Hyndman et al., 2022) and TSPred (Salles and Ogasawara, 2022). Different online resources in relation to TSA were also referred (Hyndman, 2022b; STACKOVERFLOW, 2022).

All the datasets were resampled to 25 km × 25 km grid in GIS. The forest types, i.e., TMDF, TDDF, and TEGF were extracted from the vegetation type map based on the majority class occurring in each grid. The active fire counts were extracted for each forest type on monthly basis for all the grids and all the years. The average number of dry days and maximum average temperature were extracted for TDDF, TMDF, and TEGF on a monthly basis for all the years at grid level. The geographic coordinate system and WGS 84 datum were used in the present research.

The daily rainfall and maximum average temperature data had gaps, i.e., non-availability of data for particular days or a particular month for a specified grid. These gaps were filled either by taking the average of previous and next day observations or by taking the average of two neighboring grids belonging to the same forest type and terrain conditions. Such gaps, however, were not significant and hence filling them did not affect the overall trend.

The forest type-wise monthly active fire counts, the total number of dry days, and the maximum average temperature were used as inputs to forecast the active fire counts. We avoided 0 fire counts in the model and replaced all 0 values with 1. For TDDF and TMDF and TEGF 2%, and 1% of all the observations had “no fire” incidence (0 fire counts), respectively. Further, this occurred only during the rainy season. For TEGF 3%, no fire incidences occurred during the fire season, and this was because these forests were comparatively wetter during the fire season than other forest types.

Forecasting using ARIMA model is based on the premise that the time series data is stationary, i.e., it has no trend and seasonality. The ARIMA model has three components, i.e., “AR,” “I,” and “MA.”. The “AR” or autoregression component provides information about regression of time series data with itself in different time lags. It gives clues about how significantly the data is auto-correlated, and thus the impact of past data on future forecasting could be understood (Jones et al., 2009). The “AR” analysis is also helpful in understanding whether the data is seasonal. For seasonal time series, the “AR” exhibits sinusoidal patterns. For non-stationary data, significant autocorrelation is observed even for the larger time lags, whereas, for stationary time series, autocorrelation dies down-within the initial few time lags (Mujumdar, 2022).

The “I” or integrated term provide information about how many “difference” terms are required to make the data stationary in case the time series data is non-stationary. The “difference” terms mean the difference between the current time and the past time for all the time lags.

The “MA” provides information about the number of “moving average” terms to be used. The “MA” is the average of the fixed set of previous time-period values when it is moved throughout the data sets to obtain a new time series. The “MA” can be carried out multiple times, i.e., the output of the first “MA” can again be subjected to a second “MA” and so on, till stationarity is achieved (Mujumdar, 2022).

The model belonging to ARMA family may be written as follows (Mujumdar and Kumar, 1990).

where Y(t) (t = 1, 2,….) is the time series used in the modeling; m₁ is the number of AR terms; ϕ_j is the jth AR parameter; m₂ is the number of MA terms; θ_j is the jth MA parameter; C is a constant; and w (t) (t = 1,2,…..) is the residual series.

Kashyap and Rao (1976) argued that time series differencing, which is an essential component of ARIMA, causes the variance to increase continuously, and hence such models cannot be used for data simulation. They may, however, be used for one-step ahead forecasting.

Identification of optimum “AR,” “I,” and “MA” terms is important to achieve stationarity. This process is time-consuming and at times may not converge for a higher denomination of “AR,” “I,” and “MA” terms, thus we used the Hyndman and Khandakar (2008) algorithm for automatic identification of these terms and avoided manual iterations. Further Hyndman and Khandakar (2008) believed that it was better to make as few differences as possible to avoid widening of the prediction intervals. Thus, they preferred the unit root test for defining the difference terms over minimizing the AICc (Corrected Akaike’s Information Criterion), which tends to lead to over differencing. AIC has, however, been used to select the orders of the AR and MA components in auto.arima model (Hyndman, 2022a). In this context, the concerns of Kashyap and Rao (1976) seem to have been addressed. The forecast module of “R” open-source package developed by Hyndman et al. (2022) was used to forecast the active fire counts.

Model calibration

The monthly forest type-wise time series data (2003–2017) were divided into two parts, i.e., “train data” (2003–2014) and “test data” (2015–2017). The “train data” was used for ARIMA model building, whereas, the “test data” was used to test the model performance (Walters, 2022).

The autocorrelation function (ACF) and partial autocorrelation function (PACF) of “train data” were plotted to investigate the “seasonality,” whereas, the cumulative periodogram was plotted to investigate the “periodicity” present in the time series. The ACF and PACF are the plots of time lags against autocorrelation values, whereas, periodogram analysis was carried out in the frequency domain and plotted between frequency and periodicity.

The Hyndman and Khandakar (2008) algorithm helped in automatically determining the number of “AR,” “I,” and “MA” terms for non-seasonal, i.e., p (AR), d (I), and q (MA) as well as seasonal components, i.e., P (AR), D (I), and Q (MA) of the data. The number of difference terms (d, where d is 0 ≤ d ≤ 2) were determined by repeated Kwiatkowski-Phillips-Schmidt-Shin (KPSS) test to achieve the “stationarity.” The algorithm does not search for every possible combination of different ARIMA terms, rather it uses a step-wise search to traverse through the model space. All the results were validated to confirm that the model residuals were white noise (significant autocorrelation was absent).

In the univariate model, only the monthly fire data from 2003 to 2017 were considered, whereas, in model with regressors, temperature and dry days (2003–2017) were also considered for forecasting forest active fire counts [in the Indian context dry days and temperature are important forest fire drivers (Kale et al., 2017)].

We used “xreg” argument under auto-arima function in R to consider regression errors for forecasting (Hyndman and Athanasopoulos, 2018). The overall functioning of ARIMA model with regressors is more or less same as the standard regression process.

Model validation

The best ARIMA model obtained was subjected to validation. This was achieved by observing the significant level of residuals of ARIMA model through plotting ACF, PACF, and cumulative periodogram, and by conducting the Portmanteau test (L-Jung Box test) (Box and Pierce, 1970; Ljung and Box, 1978; Harvey, 1993; Coghlan, 2018; Hyndman and Athanasopoulos, 2018). The insignificant ACF, PACF, and cumulative periodograms were desired for valid forecasts. The Portmanteau test was conducted based on the number of lags as suggested by Hyndman and Athanasopoulos (2018). Since data had seasonal components, the number of lags was equal to 2 m, where m is the period of seasonality [12 (monthly) in the present research]. The hypothesis was tested by analyzing p-values. The higher p-values (>0.05) depicted that the residuals were not distinguishable from white noise series, i.e., there was no serial correlation present (Hyndman and Athanasopoulos, 2018).

Post model validation, forest fires were forecasted for the year 2015 onward. The original (test data) vs. forecasted monthly fire were compared for the years 2015–2017 (36 values) by plotting them within the confidence bands of 80 and 95% level significance.

To estimate the forecasted fire at the grid level (to be used for FFR estimation), the percent monthly contribution of each grid of a particular forest type toward fire occurrence was averaged for 15 years (2003–2017) and the forecasted fire (for the year 2017) were segregated in the grids in proportion to their average percentage contribution of fire.

The FFR investigation for the year 2017 was carried out by taking different gridded themes into consideration. These included (maximum average monthly temperature, average number of dry days, monthly average of 15 years, i.e., 2003–2017), elevation, slope, and vegetation type (Figure 2).

The temperature, dry days, elevation, and slope grids were classified into 5 classes, i.e., 1. Low 2. Moderate 3. High 4. Higher, and 5. Highest based on the “Jenk’s natural breaks” in GIS. The theme-wise frequency ratio was estimated by taking the ratio of “fire ratio” and “grid ratio” (Lee and Pradhan, 2007; Pradhan, 2010; Biswas et al., 2015). The “fire ratio” is the ratio of fire that occurred in a particular class of a theme (for example, fire occurred in low elevation class in “elevation” theme) and the total number of fires that occurred in that particular theme (for example total fire occurred in the elevation theme). The “grid ratio” is the ratio of the total number of grids present in a particular class of a theme to the total number of grids present in that theme. The grid-wise FFR for each theme was integrated to estimate the final grid-wise FFR (Supplementary material 1), which depicted the fire susceptibility. The higher FFR depicts higher susceptibility toward forest fire.

The forecasted FFR (2017) was compared with the original FFR (2017) to understand the variations.

Results

Forest fire incidences

The TDDF, TMDF, and TEGF depicted almost similar patterns of fire occurrences during different months of the year. Fire incidences reached the peak during summers and the trough during rainy seasons. This pattern was repeated year after year (Figure 3). The forests achieved peak fire during March/April for all the studied years. No increasing fire trend was, however, observed for any of the studied forest types. The highest incidences of fire were observed in TMDF followed by TDDF. TMDF are prominent in the North-East regions of India, where shifting cultivation is prevalent in pockets, thus the anthropogenic fire incidences are rampant. The TDDF are prominent in the Deccan Peninsula and the semi-arid regions of India which are warmer and thus the fire incidences are high, particularly during extreme climate events. The TEGF had comparatively lesser fire events as forest floors were not as dry as TDDF and TMDF.

The maximum average temperature had significant seasonality with not much variation in patterns for TDDF and TMDF, whereas, fluctuations were observed in TEGF, particularly in the year 2009, when the temperature reached minimum among all the studied years (Figure 3). The number of dry days was maximum in TDDF, whereas, the range of dry days was smaller for TEGF as compared to TDDF and TMDF (Figure 3).

The ACF of forest fires exhibited a gradual decay, which indicated higher autocorrelation among the subsequent lags rather than the lags at the distal end of the time series. The ACF had a sinusoidal pattern which depicted seasonality. Significant autocorrelation (both positive and negative) was present even upto the time lag of 100, which is suggestive of strong seasonality (Figure 4). Clear PACF peaks were observed in TDDF, TMDF, and TEGF (Figure 4). This is suggestive of the number of AR terms that may be required in the ARIMA model. All the forests had significant periodicities, which were reflected in the cumulative periodogram plots (Figure 4). On many instances, the periodicities crossed the limit of the significance range, which was at 95% level. Periodicities in the data were subjected to spectral analysis in the frequency domain rather than the time domain. The “cumulative periodogram” depicted the presence of significant periodicities in all the forest types.

Autoregressive integrated moving average models and their intercomparison

A univariate ARIMA model for TDDF, TMDF, and TEGF was developed using auto.arima function. The best model obtained using Hyndman-Khandakar algorithm had p,d,q, and P,D,Q-values of 0,0,1, and 2,1,0 for TDDF; 0,0,1, and 2,1,1 for TMDF and 1,0,0 and 2,1,2 for TEGF, respectively (Table 1 and Supplementary material 2). After applying the model, the residuals had insignificant autocorrelation and partial autocorrelation for most of the time lags for all the forest types (Figure 5). The cumulative periodogram also depicted insignificant periodicities (Figure 5).

The Portmanteau test (L-Jung box test) depicted that the residuals were not distinguishable from white noise, i.e., there was no serial correlation present (P > 0.05) (Table 1). Forecasting was carried out using “train” data (Figure 6). The plot of the forecasted and original fire depicted good agreement (Figure 7). There was, however, considerable variation between the original and forecasted values for TDDF during the year 2015, and the peak original fires were not in either 80 or 95% confidence bands. This was revived in 2016 and then in 2017. The TMDF exhibited good agreement between original and forecasted fires for the years 2015 and 2017; however, for 2017 a shift was observed during the peak fire months. The TEGF had a better agreement between original and forecasted fires in all the studied years. During the year 2017, however, there was over prediction of fire. The non-fire season active fire counts were in good agreement with forecasted fire for all the years and all the forest types.

ARIMA model with regressors considered temperature and dry days as independent variables. The p,d,q and P,D,Q-values for these forests types were,0,0,1 and 2,1,0 for TDDF; 2,0,0 and 2,1,1 for TMDF and 1,0,0 and 2,1, 0 for TEGF, respectively (Table 1 and Supplementary material 2). The residuals had insignificant autocorrelation and partial autocorrelation for all the forest types and the cumulative periodogram depicted insignificant periodicities (Figure 8). Portmanteau test also depicted that residual series was not distinguishable from white noise (Table 1).

Not much variation in R² between original and forecasted forest fires was observed for the univariate ARIMA model and ARIMA model with regressors for all the forest types (Table 2). ARIMA model with regressors, however, had fluctuating minimum active fire counts during the non-fire season, which were smooth in the univariate model. This resulted in negative forecasted values (Figure 9). The overall trend of forecasted fire was more or less similar in both models (Figure 10).

On comparing the performance of the univariate ARIMA model and ARIMA model with regressors, no major difference was observed when overall model performance was concerned. A better agreement between forecasted and original active fire counts was observed for the year 1 for TDDF and TEGF when ARIMA model with regressors was used. For TMDF equal performance for year 1 was observed with both the models. Except for TDDF, the performance of both the models deteriorated from years 1 to 3. The agreement (R²) between the original and forecasted fires for the third year (2017) for TDDF and TMDF was superior in the univariate ARIMA model, whereas, for TEGF it was superior in the ARIMA model with regressors (Table 2). Among all the studied forests the third year (2017), forecast for TMDF was superior with the univariate model (Table 2).

Forest FFR provided information about fire susceptibility. For most of the forests, FFR ranged between 0 and 10. Higher FFR was depicted in the month of January, October, and December. North-East and Deccan Peninsula regions frequently depicted higher FFR in different months of the year. March and April were the prominent fire months; however, FFR did not increase significantly as the fire was well spread rather than clumped in specified regions. In May, FFR reduced in the central Deccan Peninsula and parts of North-East, which again increased in June in the central Deccan Peninsula. In October, high FFR was reported in the same region.

The original and forecasted FFRs were in good agreement for different months of the year 2017 (Figure 10). In March and May, there was underestimation of FFR in the central region of the Deccan peninsula, whereas, overestimation was observed in North-East in May. Some parts of the Deccan Peninsula had higher FFRs in October, which were also depicted in the forecasted FFR. In November, overestimation and in December underestimation were observed in the North-East region. For high and very high categories, original and forecasted FFRs were in good agreement (Figure 11). The FFRs were not forecasted for rainy season, i.e., July to September.

Discussion

The time series forecasting provides the flexibility to investigate the univariate as well as multivariate time series where the dynamic relationship of different variables are investigated. We investigated the univariate ARIMA model and ARIMA model with regressors. The univariate ARIMA is a single equation model which relates the future and past observations of a given time series, whereas, the ARIMA model with regressors is also a single equation model but also takes into consideration explanatory variables. It is based on the assumption that explanatory variables affect the response variables but not vice versa.

A multivariate ARIMA model is an n equation n variable linear model that relates to past values of different variables; thus, the variables relate dynamically in a closed loop system (Jones et al., 2009). As dry days and temperature used in the present study are causal variables of forest fire and vice versa is not true, i.e., forest fires do not cause variations in seasonal temperature and number of dry days (at least at regional scales), we found a single equation model with exogenous variable more appropriate for our research.

Even when multivariate time series models provide greater flexibility to investigate the dynamic relationships in different variables time series, on many occasions, it is practically difficult to find seamless data for multiple variables for a significant amount of time. This is where the univariate models are found useful as the data requirement is limited just to the single variable in question. As we had seamless time series available for dry days, temperature, and active fire, we experimented both with the univariate model and the model with regressors. We obtained near similar results with both the models and finally preferred the univariate model for forecasting forest fires considering the computational gains and parsimony. It would, however, be important to investigate the twofold modeling approach across different geographical set-ups to validate the findings of the present research.

It is worthwhile to argue for not using the ‘‘count’’ time series model even when the fire data was a ‘‘count’’ data. ‘‘Count’’ time series models are generally used for rare events. Typically, such models assume that the events are independent of time and thus they are memory-less.² Wildfire has a well-defined seasonality in India. Further, the fire events are not “independent,” for example, during an extreme climate event, fire incidences are significant, but in the immediate next year fires are diminished mostly due to less availability of fuel (as it was consumed by fire in the previous year and could not be replenished through regeneration). It is therefore we preferred the regular ARIMA model to “count” time series models.

Some studies, for example by Albertson et al. (2009), used the probit model to consider the chance of wildfire breaking out on a given day. Such models are useful in predicting the presence or absence of fire on a given day. The authors could correctly predict five out of seven fires in a typical year. Our study was based on the monthly time series data, and hence such a model was of limited use; however, we believe that the model may be useful in estimating the probability of occurrence of fire on a particular day, i.e., on a holiday (as studied by the authors), particularly, in the set-up where forests have significant anthropogenic pressure. Kadir et al. (2020) used the ARIMA model to predict the hotspots in the Riau province of Indonesia. They could accurately forecast the hotspots for upto 5 months using the ARMA model. They have not used any explanatory variables for the occurrence of fire. We also believe that the univariate model may be extremely useful when monthly forest fire forecasting is the objective as the computation cost of regressors can be avoided. In the present research, the performance of univariate and model with regressors was near similar; however, for 1 year ahead, the forecast model with regressors was found slightly superior. The model with regressors may further be useful in studying the time series at weekly or daily time steps as the effect of climate variables on early or late occurrence of fire season can be studied objectively.

Viedma et al. (2018) used the longitudinal negative bionomial and the ZINB mixed models to model the number of fires. The main purpose of using these models was to handle the overdispersion (variance greater than the mean) or excess of 0 values present in the data. They experimented both with univariate as well as models involving explanatory variables. In the present research, no overdispersion or excess of 0 values was observed, and hence the ZINB model was not of much use.

Different studies have used auto.arima algorithm for the identification of the best ARIMA model for forecasting (Choudhary et al., 2022; Sharmin et al., 2022). Our study as well as the auto.arima algorithm worked well, and the p,d,q, and P,D,Q parameters obtained could provide a comparable active fire forecast with that of original active fire counts.

The fire forecast in the present research was based on 15 years of available active fire data from MODIS. During this time, forests owing to climate extremities encountered many fluctuations in fires, particularly in the years 2004, 2009, and 2012, thus the fire incidences considered in the present research were a mix of normal and extreme years, and no separate scenario was formed for normal and extreme years. This is justifiable because 1. Extreme events have become so frequent that they need to be included in TSA for meaningful forecasting and 2. It is important to consider all available data to have continuity in time series.

Biswas et al. (2015) found the fire susceptibility in Myanmar in the protected and non-protected areas and suggested that human activity explained most of the variance in burn probability. They suggested that fuel composition is one of the most important factors in determining fire susceptibility. Our study is in partial agreement with theirs. The forest fires in India are mainly anthropogenic. In some pockets of North-East India, shifting cultivation is a major driver of forest fires in the fire season. This has resulted in higher fire susceptibilities in these regions. We believe that climatic factors, i.e., temperature and dry days form conducive conditions for the fire to occur. It was observed that for lower FFRs, percent average contribution of vegetation type was maximum (25%), whereas, other causal drivers, i.e., elevation, slope, temperature, and dry days contributed 20, 22, 17, and 16%, respectively. As against this for higher FFRs, the temperature and dry days were the main causal drivers. In the month of January, the dry days and temperature together contributed around average 70% of the total FFR in the highest (15–20) FFR region. A similar trend was observed for all the months from February to April, which are the prominent fire months in India. This clearly indicates that the dry days and temperature are important fire drivers in India, and the forests remain less susceptible (despite the availability of fuel) to fire when these drivers are not significant contributors.

Among all the drivers studied, the average FFR from January to April (prominent fire months) was maximum for dry days followed by the slope. Except for elevation, all the drivers depicted decreasing trend from January to April and clumped around the value of 1 in the month of April. This indicates that during the start of the fire season, different drivers contribute maximum in defining the fire susceptibility, whereas, during the late fire season they have a minimum contribution.

In this research, we have explored the possibility of using FFR as an indicator by taking forecasted forest fires as the input. The forecasted active fire counts obtained were at the level of forest types and not at the level of 25 km × 25 km grid. This is because there may be instances when there are no incidences of fire at the grid level, thus the continuous data are not available, and the forecasting may not be meaningful.

We found both univariate ARIMA model and ARIMA model with regressors as potential tools to forecast forest fires. We preferred the univariate model due to parsimony. There is, however, a need to further investigate the model performance using high temporal and spatial resolution datasets.

The TSA is a data-intensive science, and better results are expected when fire data for a significant number of years is available. This also helps in validating the model with a larger test dataset. The correct MODIS active fire data is available since 2003 and there is a need to compile the forest fire data prior to 2003 from different available sources to further improve the forecasts.

Data availability statement

Publicly available datasets were analyzed in this study. This data can be found here: https://firms.modaps.eosdis.nasa.gov/download/.

Author contributions

All authors listed have made a substantial, direct, and intellectual contribution to the work, and approved it for publication.

Acknowledgments

We acknowledge the support received by P. Bhawani related to data format conversions.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/ffgc.2022.882685/full#supplementary-material

Footnotes

1. ^ https://lpdaac.usgs.gov/
2. ^ http://utdallas.edu/∼pxb054000/code/count-examples/ECTS-I-2010.pdf

References