Operational Forest-Fire Spread Forecasting Using the WRF-SFIRE Model

Abstract

In the present research, the open-source WRF-SFIRE model has been used to carry out surface forest fire spread forecasting in the North Sikkim region of the Indian Himalayas. Global forecast system (GFS)-based hourly forecasted weather model data obtained through the National Centers for Environmental Prediction (NCEP) at 0.25 degree resolution were used to provide the initial conditions for running WRF-SFIRE. A landuse–landcover map at 1:10,000 scale was used to define fuel parameters for different vegetation types. The fuel parameters, i.e., fuel depth and fuel load, were collected from 23 sample plots (0.1 ha each) laid down in the study area. Samples of different categories of forest fuels were measured for their wet and dry weights to obtain the fuel load. The vegetation specific surface area-to-volume ratio was referenced from the literature. The atmospheric data were downscaled using nested domains in the WRF model to capture fire–atmosphere interactions at a finer resolution (40 m). VIIRS satellite sensor-based fire alert (375 m spatial resolution) was used as ignition initiation point for the fire spread forecasting, whereas the forecasted hourly weather data (time synchronized with the fire alert) were used for dynamic forest-fire spread forecasting. The forecasted burnt area (1.72 km²) was validated against the satellite-based burnt area (1.07 km²) obtained through Sentinel 2 satellite data. The shapes of the original and forecasted burnt areas matched well. Based on the various simulation studies conducted, an operational fire spread forecasting system, i.e., Sikkim Wildfire Forecasting and Monitoring System (SWFMS), has been developed to facilitate firefighting agencies to issue early warnings and carry out strategic firefighting.

1. Introduction

Forecasting forest-fire spread is important for providing early warnings to the stakeholders. Such forecasts are based mainly on the meteorological conditions, fuel properties, and topography. Wind controls the intensity and direction of the fire spread. The variability of fuels (dead or alive, fine or coarse, surface or canopy) and their moisture conditions determine the fire spread rate in horizontal and vertical directions. Topography governs the wind vector, vegetation type distribution, and moisture variability. The combined effect of these determines the fire spread in time and space. Except topography and vegetation type, other parameters, i.e., meteorology and fuel conditions, are dynamic, and thus capturing their variability frequently provides important insights about the spread of forest fires.

Atmospheric phenomena and fire spread operate at distinct scales; nevertheless, both are coupled. Atmospheric fluxes affect the fire spread, and the fluxes coming out of the burning fire control the ambient wind speed and direction (particularly for large fires). Such a phenomenon is not captured in many of the popular non-coupled models (such as FARSITE), and thus the accuracy of fire spread forecasting is hampered with the progression of time [1,2]. This is where coupled models such as WRF-SFIRE are important, where the impacts of the atmosphere on the fire and vice versa are dynamically interlinked and the fire spread’s speed and direction can be appropriately determined on a near real-time basis.

Another important distinction between coupled and non-coupled models is that non-coupled models are ‘kinematic’, which means they only capture changes in the rate of the fire spread and not the forces that are responsible for it. Kinematic models have not been successful in representing complex mountain airflows and wildfires evolved through dynamic fire–atmosphere interactions [3,4]. These forces are estimated in a coupled model and are useful to uncover the physical bases of fire phenomena [1]. New-generation wildfire and atmospheric models were reviewed by Bakhshaii and Johnson [2]. They opined that developments in satellite remote sensing and high-performance computing (HPC) will provide enormous opportunities for the operational deployment of such models. Some high-resolution models such as FIRETEC and WFDS use computational fluid dynamics (CFD) technology to carry out fire spread simulation in 3D, which is computationally intensive. The demand for computation is so high that such models run slower than real time even when executed on a parallel supercomputer [5].

WRF-SFIRE is a community managed open-source model used in wildfire spread simulation [6]. It was developed by the University of Colorado Denver and NCAR [7] and is foundation for operational systems such as the Colorado Fire Prediction System (CO-FPS) [8] and Israel national wildfire simulation and danger system (MATASH) [9]. Similarly, the coupled WRF-Fire model has been used for operational implementation and the evaluation of a rapid-response fire spread forecasting system named IRIS [10].

The Weather Research Forecasting (WRF) model provides the forecast of different atmospheric variables including temperature, wind speed and direction, and relative humidity. It operates at coarser scales for atmospheric variables; however, for fire spread simulation, it can operate even at a decameter level [2]. High-resolution fire spread simulation is achieved through the use of multiple nests, which enables step-wise downscaling of the atmospheric data to the innermost nest (the fire nest). The innermost atmospheric nest is then again divided into a high-resolution fire mesh. In numerous studies, a 1:10 ratio of innermost atmosphere nest to fire mesh has been considered; however, there are instances when a finer ratio has been used [11,12]. In the fire mesh, the fire spread is initiated based on the ignition location and an equation governing the rate and direction of fire spread. This particular phenomenon is governed by the SFIRE.

In WRF-SFIRE, the fire spread rate mechanism proposed by Rothermel is considered. The Rothermel model is quasi-empirical and based on the heat balance model developed by Frandsen [13]. According to Rothermel, the rate of fire spread is the ratio between the heat flux received from the source and the heat required for ignition by the potential fuel [14,15]. The input parameters to estimate the rate of fire spread include fuel particle properties (heat content and mineral content), fuel array arrangements (surface area-to-volume ratio (SAVR), oven-dry fuel load, fuel bed depth, and dead fuel moisture of extinction), and environmental values (fuel moisture, mid-flame wind speed, and slope steepness) [15].

Andrews [15] reviewed the Rothermel equation and provided the details of different sub-equations and corrections made for the equation mainly by Albini [16]. The only differences in SFIRE compared to the Rothermel [14] model are that SFIRE subtracts the moisture from the fuel loads in its computation rather than using upfront values; limits the slope and the wind speed; and reduces the wind to the mid-flame height [5].

Apart from adopting the fire spread rate mechanism from the Rothermel model, the WRF-SFIRE model uses the level set function for propagating the wildfire from cell to cell in the finer-resolution fire mesh. Unlike particle-based approaches where fire particles are traced for their evolution, in the level set method, the Eulerian approach is followed and the evolution of the level set function is governed by the partial differential equation [5].

The major challenge in running the WRF-SFIRE model is the availability of fuel data. Globally, primary fuel data are available only for few countries, for example Canada, the USA, and Australia. Such data are the outcome of considerable research carried out towards fuel characterization in these countries. This has resulted in the development of fuel models that can depict fuel classes and the associated fuel parameters [17,18,19,20,21,22]. Considering the variability in fuel parameters in different countries, it is not straightforward to directly use these models in areas other than those for which they were developed. It is therefore important to critically analyze region-specific fuel data to ensure the better representation of fuel in the WRF-SFIRE model.

Primary fuel data are traditionally collected through sampling the fuels (including line transect and fixed plot-based methods) and carrying out different field and laboratory-based investigations to characterize them. More recently, photograph-based fuel load estimation methods have also been developed [23,24,25,26,27,28,29,30,31,32,33,34,35,36]. Firemon (Fire Effects Monitoring and Inventory System) is an agency-independent initiative for characterizing ecosystem changes using plot-level sampling (https://www.frames.gov/firemon, accessed on 24 June 2024). Firemon provides comprehensive information about various aspects of a sample plot design and elements of the fuels that are to be investigated.

The planer transact method is based on the probability-proportional-to-size concept [37]. It is relatively fast and simple to use [28,29,30]. Fixed-area plot methods are based on the frequency concept to sample the fuels [38]. Fixed-area plot methods are considered accurate and unbiased for sampling, including for sampling fuels [39,40,41]. The photoload method, although easy to operate, may not capture all forest types adequately. The forest of a particular region may be different from the forests for which photoload references are available. This method warrants additional study before its widespread implementation across other vegetation types [37].

Different fuel parameters, i.e., fuel load, height/depth, moisture level, and SAVR, are considered important for fuel characterization. Fuel size class-based estimation has been carried out in many previous studies [42,43,44,45]. In these, the fuels are classified based on their diameter into four categories, i.e., 1 h (hour) (0–0.6 cm), 10 h (0.6–2.5 cm), 100 h (2.5–7.6 cm), and 1000 h (>7.6 cm) fuels. Fuel parameters vary significantly within fuel size classes, and thus the fuel size affects the fire spread pattern. The geometry of fuels, their size, shape, and SAVR are important for knowing the ignitability, heat transfer, and resident time of the spreading wildfire [46]. Other fuel properties, i.e., particle density, mineral content, and heat content, are generally used as constants [15].

Fuel models describe fuel categories and the associated fuel parameters. In numerous studies, landuse–landcover (LULC) maps with fuel parameters as attributes have been used as inputs for forest-fire spread simulation [47,48]. Keane [49] reviewed the approaches, methods, and systems used to describe wildland surface fuel loading. In high-biodiversity countries such as India, around 100 vegetation/landuse types have been mapped [50]. The fuel variability in terms of fine and coarse fuels is also high, and it is therefore important to develop a fuel model for India which can be readily used for forest-fire spread forecasting. To initiate this, there is a need to develop/adopt a fuel sampling approach to capture the significant fuel variability.

Satellite remote sensing is found to be extremely useful in detecting fuel variability based on the differences in reflectance patterns recorded in different bands. Numerous studies have been carried out to estimate fuel load and fuel moisture using remote sensing techniques [49,51,52,53,54,55]. Remote sensing has proved useful in live-canopy fuel characterization [56,57,58,59]; however, challenges lie in characterizing dead below-canopy surface fuels, which generally are available up to 2 m high from the forest floor. Efforts have previously been made to estimate surface fuel loads using satellite remote sensing-based inputs [60,61]. Reich et al. [61] used Landsat TM bands and forest class, elevation, slope, and aspect as independent variables to model forest fuel loading. Their models could account for 55% to 72% of the variability in forest fuels. Terrestrial Lidar Systems (TLSs) have been used in surface (below-canopy) fuel modeling and biomass estimation; however, only few studies are available [62,63]. Deep learning methods have also been found to be useful in classifying surface fuel types using stand photographs and satellite time-series data [35]; nevertheless, the direct estimation of surface fuel loads up to 2 m high using remote sensing is still a research challenge for operational fuel load estimation.

In numerous studies, the available fuel models have been used to simulate wildfire in different regions by relating the forest type of that region with the available fuel parameters as attributes [9,47,64]. Depending upon the correct adoption of the fuel model, such efforts may be useful for gaining first-hand information on forest-fire spread. It is, however, important to consider localized fuel parameters to improve the simulation accuracy and validate the usefulness of the available fuel models in different regions.

In addition to fuel characterization, the determination of the ignition location is critical for fire spread forecasting. The MODIS and VIIRS sensors onboard the Terra and Aqua satellites provide active fire locations at 1 km and 375 m resolution, respectively; however, there is a significant time gap between the two successive passes of these satellites. In the present scenario, six observations of satellite-based fire alerts are provided everyday in India by MODIS and VIIRS sensors. This may result in missing short-duration fires. As of now, there is time delay of around 1 h to 1.5 h (depending upon the sensor and processing time) in fire alerts reaching the user after the fire location is detected by the satellite [65]. Moreover, the ignition location may not be determined in overcast or smoky conditions; nevertheless, considering the presence of ground forest staff, the ignition location can often be quickly determined before the pass of a satellite.

Recently, the NASA-funded Fire Information and Resource Management System (FIRMS) has announced that it can provide active fire data within seconds after a fire is detected by a satellite. This is being achieved through the direct broadcast of these data in the continental United States (https://www.earthdata.nasa.gov/learn/articles/firms-urt-data, accessed on 25 June 2024). Such a system is needed on a global scale for the better mitigation of wildfire spread. Geostationary satellites provide information about fire locations at frequent intervals. For example, the INSAT 3D satellite (https://www.mosdac.gov.in, accessed on 25 June 2024) provides information about a burning fire’s location every half an hour; however, such satellites, due to their coarse resolution, are more effective in monitoring large fires.

The WRF-SFIRE simulations are computationally complex because of the significant data volume and use of partial differential equations. To make the partial differential equations computationally tractable, finite-difference methods have been used. In order to rapidly provide forecasting results to firefighting agencies, it will inevitably be necessary to use an HPC system.

The operational simulation of wildfire spread in HPC requires user-friendly graphical user infaces and internet connectivity. A real-time HPC workflow with the components of a graphical front end, an HPC driven simulation and management system, and web-based visualization has been developed previously (https://github.com/openwfm/wrfxctrl; https://github.com/openwfm/wrfxpy; https://github.com/openwfm/wrfxweb, accessed on 16 May 2024). Python and Java scripts were used in developing this workflow [7]. There is, however, a need to develop a user-friendly operational fire spread forecasting system exclusively for India based on the requirements of its firefighting agencies in order to cater to their needs in relation to day-to-day fire management and the dissemination of early warnings to stakeholders.

A GIS is an integral component of an operational wildfire forecast system as it provides great flexibility in visualizing the fire spread in the form of a map and making decisions for fire control on near real-time basis. A user interface, uninterrupted internet connectivity, and a server (to publish the map) are the typical components of web-based GIS visualization.

The present research is aimed at (1) characterization of surface dead fuel using field-based samples; (2) near real-time fire spread forecasting in the northern region of the Sikkim Himalayas using the WRF-SFIRE model; and (3) development of an operational fire spread forecasting system.

2. Study Area

The present research was carried out in the Indian Himalayan state of Sikkim (Figure 1). Sikkim is a predominantly forested state with 47.08% of its total area covered by forests [66]. The elevation varies from 310 m to 8586 m within the geographical area of 7096 km² (https://sikkim.gov.in, accessed on 16 May 2024). South Sikkim, dominated by the Sal (Shorea robusta) and Teak (Tectona grandis) species, is highly vulnerable to forest fires. These fires are mostly ground fires; however, there are incidences of canopy fires, particularly in Oak forests [67]. The fire season starts from November and lasts until April [68]. Many incidences of fires have been reported earlier, mainly in Grassland, Open Scrub, Sal (Shorea robusta), and Teak (Tectona grandis) forests. The highest number of MODIS fire alerts since 2006 during fire season was received in the year 2024 in Sikkim, which accounted for 22% of all MODIS alerts received from 2006 to 2024 (https://fsiforestfire.gov.in, accessed on 25 June 2024). Around 50% of these alerts were in the South Sikkim region and the remaining 50% were in West and East Sikkim together. This is suggestive of the high vulnerability of regions other than South Sikkim towards forest fire.

In Sal and Teak forests, significant litter is deposited on the forest floor, which makes them vulnerable to burn during fire season. In northern and western regions, the hills are covered with grasses of significant height (1–2 m), which when dried are prone to burn.

Owing to the complex terrain of Sikkim, the task of controlling fires is highly challenging and requires significant planning. This is where fire spread forecasting is highly relevant to enable the Forest Department to deploy resources at strategic locations to safeguard lives and properties from spreading forest fires.

3. Materials and Methods

The present research has been carried out to forecast the spread of the fire reported in the northern Sikkim region on 27 January 2023 at 0802 h UTC. The VIIRS sensor-based fire alert received through the NASA FIRMS fire alert service was used as the ignition point for carrying out 24 h fire spread forecasting on a PARAM SEVA HPC system, India using the WRF-SFIRE model.

Fuel parameters (based on field samples, laboratory investigations, and the available literature), forecasted meteorological data at hourly intervals (0.25° gridded GFS data), a digital elevation model (DEM) (12.5 m, ALOS PALSAR data), and an LULC map (as a proxy to a fuel map) at a 1:10,000 scale (sourced from the Department of Science and Technology, Sikkim, India) were used as inputs to forecast the fire’s spread.

3.1. Fuel Sampling

Fuel samples were collected by laying 23 sample plots (0.1 ha each) based on a stratified random sampling approach in different fuel categories (11 in the Forest, 6 in Open Scrub, 5 in Dense Scrub, and 1 in Grassland/Grazing Land). Past fire records, elevation, slope, and approachability were considered for laying down the sample plots. The past fire occurrences (2003 onwards) were determined based on MODIS and VIIRS satellite-based fire alerts, having spatial resolutions of 1 km and 375 m, respectively. Each plot was divided into five 1 m × 1 m sub-plots (four corners and the center of the plot). The fuel samples were collected in the sub-plots, whereas the shrub inventory was made in the 5 m × 5 m area located at the center (Figure 2). The enumerations made in the sample plots included tree girth (>20 cm), tree height (m), average canopy base height (m), canopy length (m) estimations, and fuel bed depth measurements. Canopy base height is the height of the tree from the base until the start of the canopy (Figure 2) (https://www.frames.gov/firemon, accessed on 25 June 2024). A low base height increases the forest’s vulnerability to canopy fire. The canopy length provides information about the spread of the canopy, which was measured based on the ground projection of the canopy length. Even though live fuel data were collected and organized as a part of broader fire spread research, in the present research, only the dead surface fuels were considered for the surface fire spread simulation.

The surface (downed) dead fuels were collected inside each sub-plot. The downed dead fuels were segregated into four classes, i.e., 1 h, 10 h, 100 h (fine fuel), and 1000 h (coarse fuel), based on the fuel diameter. The fuel diameter varied between 0.6 cm and 7.6 cm (1 h—0–0.6 cm; 10 h—0.6–2.5 cm; 100 h—2.5–7.6 cm) [69]. Coarse woody debris (1000 h) with diameters ≥7.6 cm were measured within the 1 m × 1 m sample plot boundary and also inside 0.1 ha plot based on the visual observations of such logs [69,70]. These logs were measured for their CBH at the small and large ends. The coarse woody debris was categorized in 4 classes: class 1 was ‘no decay’, whereas class 4 was ‘most decayed’, and all other classes were graded between 1 and 4 [29]. Litter samples were collected in each sub-plot. No samples of the coarse woody debris were collected. Fuel depth was the average length of the fuel (considering all the sub-plots) from the surface to the height of the fuel (Table S1).

The hourly categorization of fuel signifies their saturation time, i.e., 1 h, 10 h, 100 h, and 1000 h. It is presumed that finer fuels become both wet and dry in a short period of time and reach equilibrium quickly when compared to coarser fuels. This is why fine fuels play a significant role in the initiation and spread of forest fires [70].

Samples were collected separately for each class of fuel in each subplot. Each sample was weighed for its wet weight in the field. The litter samples were collected after all of the fine-fuel category samples were collected from each sub-plot. This was to avoid the intermixing of fuel samples. In total, up to 20 dead and dried fuel samples were collected from each plot (5 samples each of 1 h, 10 h, 100 h, and litter). The fuels were placed in zip-lock bags and sent to the laboratory for further analysis.

The collected samples were oven-dried in the laboratory for 24 h to estimate the dried fuel load. The plot-wise average fuel load and fuel depth were estimated by averaging all of the sub-plot level observations, whereas the vegetation-type-wise fuel load and fuel depth were estimated by averaging the plot-wise fuel load and fuel depths within that particular vegetation type. The fuel array was analyzed considering the SAVR (fine dead fuels), oven-dry fuel load, and fuel moisture content of the surface dead fuels.

3.2. Surface Area-to-Volume Ratio (SAVR) Estimation

The SAVR estimation was carried out for fine (dead) fuels (which are mainly responsible for the initiation of burning) following the method proposed by Buffachi et al. [71]. The fine fuels were either cylindrical (veins of leaves, grass streaks) or flat (leaf lamina). The SAVR was estimated for the lamina, veins, and mid-veins by extracting these components from the selected dead leaves (1 large, 1 medium, and 1 small-sized intact leaf) and taking multiple observations of the diameter (veins and mid-veins) and thickness (lamina) using a vernier caliper. The mid-vein and vein diameter was measured at 10 different locations. Five observations on each side of the mid-vein were made to measure the vein diameter. The lamina was measured at 10 different locations on the leaves (5 on each side of the midrib). For grasses, the diameter was measured at the base and at the top for 10 randomly selected grass streaks. These observations were finally averaged to estimate the diameter of the grass and the SAVR was estimated using ratio 4/d, where d is the diameter (for cylindrical fine fuels). Similarly, the SAVR of the lamina was estimated using ratio 2/t, where t is the average lamina thickness. The SAVR was weighed by the actual weight of these components. The plot-wise average SAVR was estimated by averaging all of the sample SAVR observations of that plot.

3.3. Localized Fuel Model Development

The Anderson fuel model was modified based on the LULC classes and fuel conditions observed in the study area. The fuel classes obtained in Sikkim were mapped to the nearest Anderson fuel class accordingly: the (1) Grassland/Grazing Land, (2) Dense Scrub, (3) Open Scrub, and (4) Forest classes were mapped, respectively, to the Short Grass, Timber (grass and understory), Tall Grass, and Closed Timber litter classes of the Anderson fuel model. Although the present research was carried out in forested landscapes, occasionally there is a likelihood of fire spreading from Cropland to Forest areas; thus, the Cropland and Agriculture Plantation areas were also clubbed with class 1, and secondly, Barren Rocky areas were also clubbed into class 1 owing to the possibility of occasional fires in these regions.

The averaged fuel load, fuel depth, and fuel moisture data acquired through the field investigations were linked with the respective fuel classes. The vector LULC map at the 1:10,000 scale was converted to a raster fuel map at 40 m resolution. The fuel map was categorized into ‘Fuel’ and ‘No-Fuel’ classes. The Fuel classes were those with the potential to burn, whereas No-Fuel classes had little or no potential of burning, such as the Built-up (Rural), Built-up (Urban), Core Urban, Gullied/Ravenous, Hamlets and Dispersed Households, Lakes/Ponds, Mixed Settlement, Peri-Urban, Reservoir/Tanks, River/Stream/Drain, Sand, Snow/Glacial, Transportation, and Village classes. The North Sikkim area where the fire spread forecasting was carried out had the Forest, Open Scrub, Dense Scrub, and No-Fuel classes in its vicinity.

3.4. Meteorological Conditions

The GFS forecasted meteorological data (0.25°) was used for providing the initial meteorological conditions to the WRF-SFIRE model. Five days of forecasted data were downloaded from https://nomads.ncep.noaa.gov/pub/data/nccf/com/gfs/prod/ (accessed on 25 June 2024) and used in the fire spread forecasting depending on the fire ignition time. In order to carry out fire spread forecasting at finer scales, these data were dynamically downscaled to the level of the fire mesh using triple-nested WRF domains. The purpose of this was to carry out high-resolution fire spread investigations. The outer to inner domain resolutions were scaled in a 1:3 ratio. The outer domain resolution was 7.2 km, whereas the inner domain resolutions were 2.4 km and 800 m, respectively (Figure 3). The inner domain was further divided into a finer mesh of 40 m (innermost grid-to-mesh ratio of 1:20) for carrying out fire spread forecasting (

Table 1. Domain parameters.

	Domain 1	Domain 2	Domain 3
Grid size of the domain	218 × 242	76 × 79	76 × 76
Resolution (1:3)	7200 m	2400 m	800 m
Fire grid ratio	------	------	1:20
Fire grid resolution	------	------	40 m

). Domains were created using the WRF domain wizard.

The WRF-SFIRE model required different input parameters for running the simulation. These included fuel parameters (fuel load, fuel moisture, fuel depth, SAVR, total and effective mineral contents, and moisture extinction coefficient), ignition parameters (ignition location and time of fire initiation), model parameters (domain details, parent grid ratio, projection system used, initial fire spread area, etc.), and physics parameters (microphysics, planetary boundary layer (PBL), and land surface model (LSM)). These parameters were provided to the model through three different namelist files, i.e., namelist.fire, namelist.wps, and namelist.input (Figure 4).

The model run time was adjusted based on the fire ignition time and was parameterized for a 24 h simulation. Different microphysics, planetary boundary layers (PBLs), and land surface models (LSMs) were experimented with until acceptable accuracy was obtained. In the present research, the micro-physics–Thompson scheme; Longwave/Shortwave radiation physics–RRTMG Shortwave and Longwave schemes; the sf_sfclay_physics–MYNN scheme; the sf_surface_physics–RUC land surface model; and bl_pbl_physics–Mellor–Yamada Nakanishi Niino (MYNN) Level 2.5 and Level 3 schemes were used.

The total and effective mineral contents and the moisture extinction coefficient were adopted from the available literature. The fuel data also included scalars (fire wind height and combustion heat of dry fuel) that were common to all fuel categories. The Mercator projection system was used to carry out the fire spread simulation.

The simulation was carried out in the PARAM SEVA supercomputing environment, having a peak computing power of 838 TF (https://nsmindia.in/node/183, accessed on 16 May 2024). A total of 8 nodes were used during the simulation, with each node having 48 cores.

In WRF-SFIRE, the fire spread rate is estimated based on a modified Rothermel equation (Equation (1)) [5,14]:

$$\mathrm{S}={\mathrm{R}}_0\left(1+{\mathsf{\phi }}_{\mathrm{w}}+{\mathsf{\phi }}_{\mathrm{S}}\right)$$

(1)

where ${\mathrm{R}}_0$ is the spread rate in the absence of wind, ${\mathsf{\phi }}_{\mathrm{w}}$ is the wind factor, and ${\mathsf{\phi }}_{\mathrm{S}}$ is the slope factor.

Andrews [15] comprehensively explained the Rothermel surface fire spread model along with the adjustments made [16].

The no-slope, no-wind fire spread rate (ft/min) is estimated using Equation (2):

$$\mathrm{R}=\frac{{\mathrm{I}}_{\mathrm{R}}\mathsf{\xi }}{{\mathsf{\rho }}_{\mathrm{b}}\in {\mathrm{Q}}_{\mathrm{i}\mathrm{g}}}$$

(2)

where the following applies:

I_R = reaction intensity (Btu/ft2/min);

ξ = propagating flux ratio;

ρ_b = bulk density (lb/ft3);

∈ = effective heating number;

Q_ig = heat of preignition (Btu/lb).

Equation (2) defines the quasi-empirical nature of the model. Each component of the equation is a set of numerous equations having different associated constants. In finally estimating the R, the three important components that affect the numerator and denominator of the equation are the SAVR, oven-dry fuel load, and fuel bed depth [15]. In the present research, the oven-dry fuel load and fuel bed depth parameters were estimated through ground-based fuel samples, which resulted in a better representation of the fuels in the estimation of R.

The level set method was used to estimate fire propagation after ignition. The evolution of the level set function is governed by the partial differential equation (Equation (3)) [5,72]:

$$\frac{\partial \mathsf{\psi }}{\partial \mathrm{t}}+\mathrm{S}‖\nabla \mathsf{\psi }‖=0$$

(3)

where $\mathsf{\psi }$ is the level set function, s is the fire spread rate, and t is time.

Once the model is initiated, the wind from the atmospheric model is interpolated into the fire spread rate equation. Upon the detection of active fire, the level set function is activated, which leads to the advancement of the fire across time and space. The fuel fractions are estimated after the ignition occurs, which decreases exponentially. The sensible and latent heat fluxes are then calculated, and resultant heat flux densities are averaged up to one atmospheric model cell. These densities are then fed back to the atmospheric model [5].

3.5. Fire Spread Forecast Validation

The fire spread forecast was validated through a satellite-based burnt area assessment using the Sentinel 2 satellite dataset, which has a spatial resolution of 10 m. The actual burnt area was delineated using visual interpretation method using the Sentinel 2 image of 29 January 2023 (post fire). The Shortwave Infrared band (SWIR) along with the Near Infrared (NIR) and Visible bands were considered for delineation of the burnt area. This was then compared with the forecasted burnt area. Data on the wind speed and directions (on an hourly basis) with similar timings to those of the fire spread were used to analyze the progression of the forest fire in comparison to the wind velocity and direction.

The fire wind vectors were derived based on the fire-mesh-level fire wind observations (UF and VF) (https://www.esri.com/arcgis-blog/products/product/analytics/displaying-speed-and-direction-symbology-from-u-and-v-vectors/, accessed on 25 June 2024).

Magnitude and direction of wind was estimated using Equations (4) and (5) respectively.

$$\mathrm{M}\mathrm{a}\mathrm{g}\mathrm{n}\mathrm{i}\mathrm{t}\mathrm{u}\mathrm{d}\mathrm{e}=\sqrt{{\mathrm{u}}^2+{\mathrm{v}}^2}$$

(4)

$$\mathrm{D}\mathrm{i}\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}=\mathrm{a}\mathrm{t}\mathrm{a}\mathrm{n}2(\mathrm{v},\mathrm{u})$$

(5)

where u and v are horizontal and vertical wind components, respectively. It was assumed that in no-wind conditions, the wind vectors were directed towards the north and afterwards changed direction as per the values obtained through Equation (5). In the present research, oceanographic wind direction conventions were followed to define the wind direction.

3.6. Operationalization of Forest-Fire Spread Forecasting

The fire spread simulation was operationalized by developing a web-based software system for fire spread simulation utilising remote HPC through internet. The open-source tools and technologies used were React (17.0.2), Leaflet (3.2.1), Java (11.0.23), Springboot (2.5.4), and Geoserver (2.23.2). React was used to build the user interface, Leaflet is a JavaScript library that was used to develop interactive maps, and Geoserver was used to share the geospatial data, whereas Postgres SQL was used as a database. In this system, the query is fired first to the OFRA (orchestration framework), through which it is directed to the connected HPC system (Figure 5). The user is supposed to provide the fire ignition location, the date and time of ignition, and the desired forecast duration (up to 48 h). The forest-fire ignition location is routed through the Forest Survey of India (FSI) server, which is based on the MODIS and VIIRS-based active fire alerts. Users also have the flexibility to feed in ignition details received from sources other than satellites.

Once the parameters are provided by the user, the namelist.input and namelist.wps files are programmatically modified and sent as a zipped file to the OFRA. The scripts were developed to carry out simulation via the HPC using static datasets (fuel model and DEM) and dynamic datasets (atmospheric (GFS) data). Once the simulation is fired, its status is displayed through the OFRA, and after the completion of the simulation, the fire spread area is converted to geotif format and sent back to the user through the OFRA, where it can be displayed in the GIS on the dashboard and the potential fire spread up to the next 48 h can be investigated in relation to other GIS layers, i.e., road network, important government office locations, settlement, etc. The user can then send SMS alerts to the registered stakeholders in the specified buffer region around the potential fire spread area.

4. Results

The fuel parameters obtained through field-based investigations and gathered from previous studies were linked with the fuel map to carry out the fire spread simulation. Hourly analysis depicted the pattern of fire spread in the region. The impact of wind speed and direction and the terrain was analyzed on an hourly basis. An operational forest-fire spread forecasting tool was finally developed to carry out forest-fire spread simulation via remote HPC using simple internet connectivity.

The highest fuel load estimated was for the Forest, followed by Open Scrub, Dense Scrub, and Grazing Land/Grassland. Open Scrub had the highest average fuel depth of 0.91 m, followed by Dense Scrub, Grazing Land/Grassland, and Forest. The Forest had a lower fuel depth; however, due to the presence of more coarse fuel and significant compact litter, the fuel load was highest for Forest. The SAVR was maximum for Open Scrub and minimum for Forest. The SAVR is an extremely important parameter for estimating fire spread. Fine fuels have a high SAVR when compared to that of coarse fuels; therefore, areas dominated by grasses burn rapidly when compared to areas that have coarser fuel types. The present study is part of ongoing research where field sampling investigations are still ongoing. The samples related to the SAVR have been collected and are being analyzed for different vegetation types. Due to this, reference SAVR observations from previous studies have been considered in the present research. The SAVR is different for different grasses. Overholt et al. [73] and Cheney et al. [74] reported SAVRs of 9770 m⁻¹ and 12,240 m⁻¹ for grass species Eriachne burkittii and Themeda australis, respectively. In the present research, an SAVR of 12,300 m⁻¹ has been considered for the Open Scrub class as the surface fuel of this class is covered with grasses that are morphologically similar to Themeda australis (for the purpose of fire spread simulation). The SAVRs for the Grazing Land/Grassland, Dense Scrub, and Forest were adopted from Scott and Burgan [19] (

Table 2. Fuel parameters used in this study.

Fuel Class	Fuel Load (kg/m2)	Fuel Depth (m)	Surface Area-to-Volume Ratio (ft−1)
1. Grassland/Grazing Land/Barren Rocky/Cropland/Agriculture Plantation	0.16	0.305	1631
2. Dense Scrub	0.234	0.61	1371
3. Open Scrub	0.243	0.91	3750
4. Forest	0.270	0.061	1224

). The fuel load and fuel depth observations were based on the data obtained through field sampling.

In addition, other parameters, i.e., fuel moisture extinction, fire wind height, combustion heat of dry fuel, fuel density, and total and effective mineral contents, were used as simulation inputs (

Table 3. Model input parameters used in this study.

Fuel Property	Values	Reference
Fire wind height (m)	6.096	[5] (https://wiki.openwfm.org/wiki/Namelist.fire, accessed on 16 May 2024), [47]
Combustion heat of dry fuel (J/kg)	17.433 × 106	(https://wiki.openwfm.org/wiki/Namelist.fire, accessed on 16 May 2024)
Fuel moisture of extinction (kg/kg)	Grassland/Grazing Land/Barren Rocky: 0.12 Dense Scrub: 0.15 Open Scrub: 0.25 Forest: 0.30	(https://github.com/openwfm/WRF-SFIRE/blob/master/test/em_fire/namelist.fire_fmc, accessed on 16 May 2024), [16,17]
Fuel density (lb/ft−3)	32 (as fuels were considered solid)	(https://wiki.openwfm.org/wiki/Namelist.fire, accessed on 16 May 2024) [15]
Total mineral content (fraction)	0.055	[15]
Effective mineral content (fraction)	0.010	[15]
Wind reduction factor (Windrf)	Grazing Land/Grassland: 0.36 Dense Scrub: 0.36 Open Scrub: 0.44 Forest: 0.36	[75]

). These have generally been used as constants in numerous previous studies. The wind reduction factor (windrf) data were adopted from Baughman and Albini [75]. The BEHAVE system measures wind at 6.096 m and log-interpolates it to mid-flame height. A similar effect is created when the wind speed is multiplied with windrf (https://wiki.openwfm.org/wiki/Vertical_wind_interpolation, Accessed 20 May 2024), [5]).

The fuel model (map) used in the present research was developed by linking the fuel parameters with the modified LULC map (Figure 6).

After the detection of the ignition point, forecasting was carried out on a PARAM SEVA HPC system with a theoretical peak floating-point performance of 838 TFLOPS. The computer nodes included 2X Intel Xeon Cascadelake 8268, with 24 cores and 2.9 Ghz processors per node. In total, 240 cores were used to carry out the 24 h forecasting, which was completed in 2 h 11 min. (

Table 4. Simulation specifications.

Number of cores	240
Fire simulation duration	24 h
Fire simulation date and time (IST)	27 January 2023, 1330 h
Ignition point location	88.67, 27.63
Computational time	2 h 11 min

).

The simulation resulted in hourly forecasted geotiff files of 24 h, which were converted to vector files (.shp files) and investigated for each hour.

The ignition location was the epicenter of the fire’s spread. The burning started in the Open Scrub. The fire moved towards the northern, southeastern, and eastern directions in the first hour. It further spread in the second and third hour in the same directions. In the fourth hour, the fire moved northwards as well as southwards. From the fourth hour to the sixth hour, the movement of the fire was southwards and southeastwards. From the seventh hour to the ninth hour, the fire movement was southwards as well as westwards. From the 10th hour to the 12th hour, the fire moved mainly towards the south and southwest direction. From the 13th hour onwards, the fire spread was mainly towards the southwest direction. From 15 h onwards, there was not much fire spread was obsrved (Figure 7). The forecasted burnt area (1.72 km²) was validated against the satellite-based burnt area (1.07 km²) obtained through Sentinel 2 data. There was a good match between the shapes of the real and forecasted fire spreads, barring the spread in the southwest direction.

Out of the total burnt area after 24 h, 76%, 16%, and 8% of the area belonged to the Open Scrub, Dense Scrub, and Forest respectively. The surrounding region had a significant amount of Dense Scrub; however, whenever Open Scrub was available, it was the preferred burning class. Around 98% of the total burnt area was burnt in the first 15 h itself (Figure 8). Thus, it was forecasted that the fire was likely to remain for mainly up to 15 h, and from 15 to 24 h, there would only be a marginal burn.

Around 85% of the Open Scrub and 91% of the Dense Scrub were burnt in the first seven hours (Table S2). This was mainly due to the fact that the Dense Scrub was localized mostly in the eastern part and had an elevation variation of 400 m, compared to the Open Scrub, which had an elevation variation of 1700 m and was geographically well spread; thus, despite it being a preferred burning class, it took more time to burn than the Dense Scrub (Figure 9). Along the western side, the Dense Scrub was located in complex terrain and hence was saved from burning. The burning of the Dense Scrub and Open Scrub peaked in the second hour and gradually declined thereafter. For the Forest, the burning peaked in the ninth hour, declined in the tenth hour, again increased in the eleventh hour, and then declined. The fire reached the Forest areas in the second hour in the northern direction; however, it could not burn aggressively due to the lesser SAVR, higher moisture of extinction, and low wind intermixing inside the Forests.

The wind maps depict the speed and direction of wind. The wind speed varied between 0 and 8 m/s. The 0 wind speed indicates calm wind. The arrows depict the direction of the wind. The straight arrow depicts the north direction and other arrows depict the directions in which the wind was flowing.

The ignition was initiated in the Open Scrub. This area was composed of tall grasses (up to 2 m) and its SAVR was higher (3750 ft⁻¹) than the SAVR of the surrounding fuel classes. This resulted in a comparatively faster fire spread. Secondly, the grasses were compact and thus the heat transfer was efficient, and the area being open resulted in a proper mixing of air with the fuel. The fire wind had a velocity of 2 m/s in the northeast direction, which resulted in spread of the fire in the northeast direction covering an area of around 0.11 km² in the first hour. During this time, the fire encountered Dense Scrub having an SAVR of 1371 ft⁻¹, which was less than that of the Open Scrub; however, since the Open Scrub shared a large perimeter with the Dense Scrub and the direction of the wind was also towards the Dense Scrub, it could burn some of the Dense Scrub as well. The aspect of the region was southern and the slope was moderate (13° to 57°). The fire started at an elevation of 2700 m and could move 100 m southwards and 200 m northwards (Figure 10).

In the second hour, the fire moved towards the northern and eastern directions, and the increase in wind velocity and gentle sloped conditions resulted in a rapid increase in burnt area (0.41 km²). There was 400 m of northward movement and 400 m of southward movement of the fire. The high wind velocity eastwards resulted in a rapid burning of Dense Scrub; however, the fire moved southwards even when the wind velocity was eastward. This was due to the ample availability of fuel with high SAVR (i.e., Open Scrub). The higher fuel load of Open Scrub when compared to Dense scrub resulted in a larger fire spread rate in the Open Scrub, as high fuel load results in a higher release of heat, which accelerates the preheating of un-burnt fuel ahead in the fire spread path [1,76].

In the next hour, the fire could not reach more than 3200 m elevation northwards. This was just a 100 m northwards elevation gain from the previous hour. The fire, however, could move rapidly eastwards at a comparatively lower elevation, i.e., 2300 m to 2400 m. This movement was supported by the eastward flowing wind. The fuel was Open Scrub, and this kept burning until the barrier of the river was reached. There was not much burning of Dense Scrub, mainly due to high slopes and a reduction in wind velocity compared to the previous hour.

Subsequently, two-way movement, i.e., eastward and southward movement, of the fire was observed. The wind flew eastwards and supported the burning of Open Scrub. The entire perimeter of fire then kept moving southward. Some Dense Scrub in the northeast direction was also burnt. The fire moved southwards despite the highly sloped terrain. The fire moved southwards from 2300 m to 1900 m mainly due to fact that there was a continuous tract of Open Scrub in this region. The movement of wind was eastward; however, its velocity was less, i.e., 1–2 m/s, which could drive the fire eastwards barely until the next fuel type with a lower SAVR and a higher fuel moisture of extinction (i.e., Forest) was reached. The fire could not further spread eastwards; however, it kept spreading southward.

For the next few hours, the fire moved in the southwest direction despite the wind direction being north and northeastwards. This was because the wind was almost calm with a very low velocity and hence there was hardly any impact of wind on the fire’s spread; secondly, most of the fuel eastwards was already burnt and the Open Scrub in the east was surrounded by Forest areas with a low SAVR and a high fuel moisture of extinction, and the wind did not have the strength to carry the fire/heat eastwards. On the other hand, Open Scrub was available until the banks of the river southwards and hence fuel kept burning in the southeast direction. After 15 h of burning, there was no significant subsequent burning as most of the fuel was already burnt and the river barrier was reached southwards and Forest was reached westwards (Figure 10).

For 24 h, the mean wind velocity was 1.44 m/s in the region of the forecasted burn. Wind with a velocity of up to 2 m/s flew mainly in the east, followed by the northeast and north directions. The comparatively higher wind speeds (2–4 m/s) in the first three hours of ignition provided ample thrust for the spread of the fire in the northeast direction (Figure 11). The eastward high-speed winds were observed from the 21st to 24th hours; however, by that time, most of the Open Scrub and Dense Scrub fuel had been burnt and the fire reached the Forest.High fuel moisture of extinction and complex topography, prevented the further burning eastwards.

Besides wind and fuel properties, the complexity of terrain also plays an important role. The Dense Scrub present on the western side could not be burnt, mostly because it was present in the gorge and more moisture was retained, which prevented the burning. A similar fuel type present on the northeast side was burnt as the terrain was less complex.

Operational Fire Spread Forecasting System

The operational fire spread forecasting system ‘Sikkim Wildfire Forecasting and Monitoring System’ (SWFMS) provides flexibility to the user (Forest Department personnel) to carry out fire spread simulation on remote HPC machines using simple internet connectivity. The authorized user can access the system through their login and provide the required details, i.e., the ignition location, the fire start date and time, and the forecast duration. At present, forecasts of up to 48 h can be carried out (Figure 12). The information provided by the user is programmatically updated in the different namelist files and sent instantly as a zipped file to the OFRA. Simulation is activated after such information is received in the HPC system through the OFRA. The forecasted meteorological data are continuously downloaded and deleted so that a 48 h fire spread forecast can be carried out at any point of time. The fire ignition location is directly ingested into the system through the FSI application. The flexibility has been provided to manually enter the ignition location based on the fire detected by the ground forest staff. This results in initiation of the fire spread simulation in the shortest possible time, without waiting for a satellite to pass to obtain the ignition location.

The simulation outputs are geotif files depicting the fire spread area. The processed output for the defined forecast duration is diverted back to the user. The outputs are programmatically converted to the Web Mercator projection system for web visualization of the forest-fire spread. The user, after selecting the date of simulation, can visualize the simulation on the GIS dashboard along with different GIS layers (Figure 13). Once the fire spread is displayed in the GIS, SMS alerts can be sent to stakeholders in a user-specified buffer region (Figure 14). This is to reduce undesirable panic among the larger population. Since the SMS alerts are highly localized, they are extremely helpful for taking quick action to control the spreading fire.

5. Discussion

The fuel sampling in the present research provided the variability in fuel load, fuel moisture, and fuel depth. Although the sample plots were allocated considering the best possible representation of terrain, having a higher number of sample plots may help in stabilizing the variability of these fuel parameters. Secondly, the fuel model has been developed considering the LULC classes. An alternative approach may be to lay down a significant number of sample plots unless the variability in fuel properties is stabilized and then derive a classification considering the empirical models using satellite remote sensing. This approach, however, is costly and time-consuming.

It is important to have synchronization between the fire mesh, fuel map, and DEM resolutions. In the present research, the nest sizing was carried out to achieve a fire mesh (grid) resolution of around 40 m.

Grid resolution has an impact on the time required for the simulation [77]. If the resolution is coarser, the number of grid cells to be processed is less. When the resolution increases, the number of cells to be processed also increases for all the domains considered. This magnifies the computation task and increases the simulation time; however, it will likely improve the accuracy of the simulation. It is, however, not straightforward that reducing the grid size results in faster simulations. Considering the number of cores/nodes used, it is required to have a minimum number of cells to be processed per node; thus, if very small domains are used, the solution may not converge. Moreover, the innermost domain should be of sufficient size to cater for the geographical area, large enough to allow the spread of forest fire. In the current study, domains and grid resolutions were defined considering the potential fire spread, optimum simulation time, and accuracy of simulation. The innermost fire mesh (grid) resolution of 40 m was found to be suitable considering the heterogeneity of the LULC/vegetation types in the study area. The innermost domain size was considered based on an investigation of the historical fire spread area in the region and optimum number of cells required for the computation.

Secondly, the calculation of physical processes takes a lot of time for running the model and hence optimum parameters need to be defined considering the region where the simulation is being carried out. Rim et al. [78], while carrying out their WRF-SFIRE-based fire spread simulation, disabled the options for physical processes in order to improve the simulation time. In the present research after considering different physics parameters, the parameters that provided better simulation accuracy were considered.

Besides defining the fuel classes, it was extremely important to depict No-Fuel classes, which acted as barriers to the fire’s spread. The high-resolution and high-frequency ‘Fuel’ and ‘No-Fuel’ class determination is an important area of research. As the forest fire initiation is a high-resolution phenomenon, the 40 m resolution optimally represented the spread of forest fire. Low resolution fire mesh resulted in inappropriate representation of ‘No fuel’ classes in particular which acted as barriers to the fire spread, thus negatively impacting the accuracy of fire spread simulation. On the other hand, high resolution fire mesh was found computationally more demanding resulting in higher computation time.

Even though the temporal and spatial resolution of satellite datasets is increasing, converting the spectral information to thematic classes to regularly derive a high-resolution multi-class LULC product has not been carried out to date. Recently, Brown et al. [79] developed the ‘Dynamic World’ near real-time global 10 m LULC product, the first of its kind, which provides extremely up-to-date continuous LULC data. These data, however, provide information about eight spectrally distinct LULC classes and hence are of limited use for high-resolution fire spread simulations in regions that have a variety of LULC classes. It is, however, extremely important to develop automatic high-frequency, high-resolution LULC products, preferably using Artificial Intelligence (AI) and machine learning (ML)/deep learning (DL) approaches, which at present are not available. In the hilly terrain, landslides are frequent; however, as they are not dynamically captured in the LULC map (Fuel–No-Fuel map), may affect the accuracy of forest-fire spread forecasting. This is because the fire spread forecast model will consider a landslide as a fuel and will allow the fire to spread, whereas on ground, landslides may act as fire barriers. Such incidence was also observed in the present research while carrying out different experimental simulations.

6. Conclusions

Forest-fire spread forecasting is important for taking effective mitigation measures to safeguard lives and properties. It is particularly relevant in areas with complex terrain where resource deployment during emergencies is challenging.

In the present research, an operational Sikkim Wildfire Forecasting and Monitoring System (SWFMS) was developed to enable the user to carry out fire spread simulations via remote HPC computing facility using simple internet connectivity, visualize the fire’s spread in GIS, and send SMS alerts in the specified buffer.

The WRF-SFIRE model was found to be effective in forecasting wildfire spread in the hilly terrain of the Sikkim Himalayas. The fuel characterization carried out in the present research is an important step towards generating the regional primary fuel model.

High-resolution vegetation type/LULC data at frequent intervals needs to be developed to represent fuel barriers in effective ways. This is particularly important for hilly terrain where landslides are frequent. Landslides act as fire barriers, however, as they are not instantly updated in the LULC map, it hampers the simulation accuracy.

Fire spread forecasting provides a clear idea of the locations where intervention from the Forest Department is critical. This helps in strategizing resource allocation. This is particularly important for complex terrain like the Himalayas, where resource deployment and ground operations are extremely challenging.

Unlike other countries, i.e., the USA, Australia, Canada, Chile, and Portugal, where devastating wildfires have been encountered earlier, in India, fires of this magnitude have not been reported to date; however, as per the study carried out by Barik and Roy [80], in an Indian context, future fire seasons are likely to be longer by 3–61 days and pre-monsoon fires will become more intense in 55% of forests. These are alarming signs and are likely to result in more fire events in the future. In this scenario, it is important to develop a near real-time National Fire Spread Forecasting and Monitoring System that will help forest-firefighting authorities in strategic resource deployment to save lives, properties, and precious natural resources.