THE ACCURACY OF GILILAWA DARAT WILDFIRE SPREAD ESTIMATION USING BURN SEVERITY AND WRF-SFIRE MODEL

In this research, the origin of the wildfire area is assessed by using the potential of burn severity and the WRF-SFIRE model. This research focuses on the mountainous savanna region, by taking the case of the Gililawa Darat wildfire event. The most accurate index among Sentinel-2B optical burn severity indices, i.e. dNBR, MIRBI, dMIRBI, CSI, dCSI, NDVI, dNDVI, EVI, and dEVI and among Landsat-8 OLI/TIRS thermal and mixed burn severity indices, i.e. LST, dLST, LST/EVI, and d(LST/EVI) was used to map the areas with low burn severity, an indication generally found at origin area. A series of fire spread simulation from these areas was conducted using WRF-SFIRE to assess the accuracy of each simulation in reproducing the burned area. The burn severity accuracy assessment showed that dNBR and dCSI indices had the highest value of Overall Accuracy and Kappa Hat Coefficient, i.e. 91.67% and 0.889 (almost perfect agreement). However, dNBR was the most suitable index for mapping burn severity in the region due to its goodness-of-fit measure for linear regression model with the R-squared value of 0.7856. The assessment of thermal and mixed burn severity indices based on Landsat-8 OLI/TIRS resulted in LST, LST/EVI, and d(LST/EVI) gained the same overall accuracy of 58.33% and Kappa Hat Coefficient of 0.444 indicating moderate agreement, whereas dLST performed poorer than these indices. However, it is not recommended to use these burn severity indices in the region due to the nonlinearity of severity level with the index value. According to WRF-SFIRE simulations result, it was found that fire ignition started from low burn severity area coordinates which have a distance of 0 to 334 metres from the origin area resulting in fire area witan h overall accuracy value range from 77.04% to 81.90% and Kappa Hat Coefficient value range from 0.536 to 0.626. The simulation from the origin area resulted in an overall accuracy of 81.57%, a Kappa hat coefficient of 0.613, underestimated burned area ratio of 0.08, overestimated burned area ratio of 0.23, and a backing fire perimeter difference ratio of 0.4 to the reference.


INTRODUCTION
Wildfires have many environmental, social, and economic consequences (Jhariya and Raj, 2014).Wildfire investigations must be conducted to determine the origin and causes of fires so that fire events do not reoccur or become less common (Bureau of Indian Affairs, 2019).In fire investigation, the observations of damage after the fire are frequently independent of the path taken by the fire, making it difficult to determine where the fire started (Gorbett et al., 2015).When the flames became too intense, it became difficult to read any patterns.Observing prominent and easily visible fire patterns from a distance appears best suited for areas with high fire intensity (Simeoni et al., 2017).
Burn severity which displays changes in objects and fuel when affected by a fire can be used as a fire pattern indicator.Satellite observation has a synoptic view capability to identify burn severity from a distance.This capability can view the fire patterns easier than close observation in the field.Remote sensing burn severity estimation can overcome extensive field surveys' cost and logistical constraints (Franco, et al., 2020).
The various fire pattern indicators, combined, analysed within fire behaviour, and compared with witness statements, can lead to the ignition area (National Wildfire Coordinating Group, 2016).Weather Research and Forecasting Model (WRF) with a Spread Fire model (WRF-SFire) can calculate the average and direction of the fire using surface wind data from interpolated atmospheric models and land characteristics such as the type, amount/mass of land cover, moisture content, and elevation data (Lai et al., 2020).
Burn severity of a wildfire is essential for post-fire assessment and wildland management decisions (Barkley, 2019), and WRF-SFIRE can assist authorities in issuing fire bans and allocating firefighting and fire prevention resources (Mandel et al., 2014).Despite the current common uses, the potential use of the methods to determine wildfire origin areas was assessed in this research.The proposed method was assessed using the Gililawa Darat Island wildfire data, which occurred on 1 st August 2018.

Study Area
The research focused on large wildfires that burned Gililawa Darat island (Figure 1) on 1 st August 2018 in a dry season.The fire was reported at around 19.00 Central Indonesia Time (CIT) and extinguished at 03.10 CIT on 2 nd August 2018 (Balai Taman Nasional Komodo, 2018).

Figure 1. Study Area
Gililawa Darat isa small island with an area of 40 hectares (from 8° 27' 39.5''S to 8° 28' 39.3''S and from 119° 33' 5.1''E to 119° 34' 1.1''E.The island is a part of Komodo National Park in East Nusa Tenggara, and This island is included in the administrative area of the West Manggarai Regency.It has a hot and dry climate, characterized by savannah vegetation (Ramono et al., 2000).

Materials
In this study, field observation data were obtained from the investigation that was conducted by a Forensic Investigator Team from the Denpasar Branch of the Indonesian National Police Forensic Laboratory Centre and West Manggarai Regency Police on the 5th and 6th August 2018.The Composite Burn Index (CBI) was used to estimate burn severity in the field.The severity is divided into four levels: unburned, low, moderate, and high (Key and Benson, 2006).A total of 30 plots (20 m x 20 m) were sampled within the perimeter of the burned area.
Based on the closest sensing date to the wildfire event and cloud free image, the two best datasets suitable for the study were selected from sentinel-2B images with sensing date 29 th July 2018 and 8 th August 2018.In the same way, the two datasets of Landsat-8 Operational Land Imager (OLI)/ Thermal Infrared Sensor (TIRS)  (Anderson, 1982).The National Centers for Environmental Prediction (NCEP) Global Forecast System (GFS) 0.25degree global forecast grids historical archive data (ds084.1)was used as meteorological gridded data.
(5) 5. differenced Char Soil Index (dCSI). (6) The continuous burn severity indices dataset was stratified into severity levels based on threshold levels to simplify the burned area description and comparison.Using the mean of the median index values for each consecutive pair of fire severity classes: unburned-low, low-moderate, and moderate-high, threshold values were calculated (Tran et al., 2018).The thresholds were adjusted until the highest percentage of samples was correctly classified (Franco et al., 2020).To obtain optimal thresholds, several trials were conducted for each spectral index (Marino et al., 2016).The most suitable burn severity index was determined using error matrix and was used to map the areas with low burn severity, an indication generally found at origin area (The National Wildfire Coordinating Group, 2016;Simeoni et al., 2017).
A series of fire spread simulation from the low burn severity areas was conducted using WRF-SFIRE to assess the accuracy of each simulation in reproducing the burned area.Multiple domains were run simultaneously at different grid resolutions in a two-way nested run.The coarse-to-fine grid ratio was five, thereby, atmospheric initialization (about 25 km) could be scaled down to the fire grid resolution (about 200 m).The fire model (coupled with the innermost domain) ran with mesh step 20 m.The time step for the inner domain and the fire model was 1,2 s, while for the outermost domain d01 it was 150 s.Initial conditions for all domains and boundary conditions for the outermost domain d01 were extracted from the NCEP GDAS/FNL 0.25 Degree Global Tropospheric Analyses and Forecast Grids (ds083.3).Model initialization is placed at 09:00:00 on the 31st of July 2018 UTC, more than 24 hours prior to the fire ignition, while the simulation time window is 9 hours.The model obtained files containing FIRE_AREA variable and it was used to asses the accuracy of the model by comparing the fire perimeter of the model with the reference (the most suitable burn severity index map).

Result
In dNBR, MIRBI, dMIRBI, CSI, and dCSI map, the locations of unburned areas corresponded accurately to the actual situation which covered almost only southwest region of the island.dMBR and dCSI were very similar in distributing all burn severity levels at entire region.In MIRBI and dMIRBI maps, more low severity areas were identified than those of dNBR, CSI, and dCSI maps.The wrong location of unburned areas among the burned areas was identified in CSI map.NDVI, dNDVI, EVI, and dEVI showed the worst scattering of unburned and burned areas at wrong locations.However, NDVI and dEVI could recognize unburned areas better that dNDVI and EVI. Figure 2 shows difference in the distribution of burn severity levels when comparing the burn severity indices maps.A large difference is also observed between thermal burn severity indices, i.e.LST and dLST and mixed burn severity indices, i.e.LST/EVI and d(LST/EVI) which were based on Landsart-8 OLI/TIRS datasets as shown in Figure 3.In the LST and dLST map, large areas of moderate burn severity (orange) were clustered in the middle part of burned areas.
The unburned area was overestimated in LST and dLST map covering the south east part of the island which was burned.In dLST map, there were unburned and low burn severity  Burn severity estimated by dNBR and dCSI resulted in an overall accuracy of 91.67% and the user's accuracy for the four burn severity levels was from 75% to 100% as shown in Table 1, indicating that these indices were good for estimating burn severity.The Kappa Hat Coefficient of these two indices were 0.889 indicating the almost perfect agreement.dNBR and dCSI performed especially well for discriminating low severity areas with unburned, moderate, and high severity areas.dNBR and dCSI had slight difficulty in separating moderate (a producer's accuracy of 66.67%) and high (a user's accuracy of 75%) burn severity levels.The second most accurate map was the CSI map with an overall accuracy of 75%, it did not identify all severity levels areas perfectly.The accuracy assessment on the optical burn severity indices that were based on Sentinel-2B datasets, dNBR and dCSI had the highest accuracy level and could discriminate areas with low burn severity well.However, dNBR was more convincing to be used obtain fire ignition start coordinate inputs for WRF-SFIRE model by considering the regression analysis of training plots burn severity value that showed the highest R-squared value.Whereas the thermal and mixed burn severity indices that were based on Landsat-8 OLI/TIRS datasets was not satisfying enough to be used due to the incapability of all indices in discriminating low severity areas well.
Fire ignition start coordinate inputs of the WRF-SFIRE model were selected from 12 coordinates of the low burn severity area of dNBR index map representatively.Coordinates 5, 6, 8, and 9 was closer to coordinates 7 that represented the origin area that the other coordinates.This selection was carried out with due observance of cardinal and ordinal points to assess the trend of the final fire perimeter, as depicted in Figure 4.
Using dNBR map as a reference, the error matrix was arranged to evaluate the fire area agreement between reference and WRF-SFIRE estimation as shown in Figure 5.Total 3868 points coordinates that were represented all pixels in dNBR map were used to arrange the error matrix.The assessment was carried out by inspecting whether the coordinates covered by WRF-SFIRE fire area output or not.As shown in Table 3, even though the fire area result of the simulation which used the valid fire ignition start coordinate (coordinate 7) had a good agreement with reference (the overall accuracy of 81.57% and the Kappa hat coefficient of 0.626), this performance was still below that of the simulation which used a false fire ignition start at coordinate 9 (the overall accuracy of 81.90% and the Kappa hat coefficient of 0.613).However, the tendency of better agreement between reference and estimation for the simulations that used fire ignition start coordinates near the valid fire ignition start coordinate could be used to be initial assessment of the most potential valid fire ignition start coordinate.It was apparent that the simulation that used the fire ignition start coordinate closed to origin area had the overall accuracy around 80% and the Kappa Hat Coefficient above 0,5 (moderate agreement).Although simulation from coordinate 7 generated larger overestimated burned area than those of coordinate 9, however the underestimated area result was smaller.The overestimated burned areas of simulations using coordinates 11 and 12 were at the fire's west flank, which was propagated mainly by lateral fire in the southwest and was also contributed by backing fire in the south direction.These fire spreads were resulted from placing the fire ignition start coordinates 11 and 12 on the west flank of dNBR's burned area perimeter, therefore, the southwest wind blow caused lateral fire to spread outward from the fire perimeter of the dNBR.The incorrect placement of fire ignition starts at coordinate 5 also resulted in the overestimation of fire area.However, coordinate 5 should not be taken into consideration due to the location had not been burned yet at least until 16.00 UTC according to the information from National Park personnel.On the other hand, all simulations except those that used coordinate 1 and 2, resulted in underestimation of fire area at east region of the island.The model seemed to be unable to capture the lateral spread at the flank side of the fire.Moreover, all simulations could not reach the backing region that located at the south edge of the island.However, the fire simulations that started from the valid fire ignition start coordinate (coordinate 7) and the nearby coordinates (coordinate 4, 5, 6, 8, and 9) showed the shape and extent of backing region that was closer to the reference.Coordinate 4 could be excluded from the analysis due to it had not been burned yet at the beginning of the actual event, therefore coordinate 7 generated the closest backing fire region to the reference.Assessment of overestimated burned area, underestimated burned area and backing fire difference in Table 5.10 shows that the simulation from coordinate 7 gained the smaller backing fire difference and underestimated area than those of coordinate 9.There was a tendency of the farther the coordinate from fire origin area, the smaller the differences.
Based on the analysis of error matrix of burned area and fire vector spread particularly backing fire spread, the simulation that used the valid fire ignition start coordinate (coordinate 7) reproduced the fire area that was closest to the reference than those of using the false ignition start coordinate

Discussion
Moderate and high severity area spectral signature that overlapped at SWIR band as shown in Figure 5 was needed to be paid attention because it could bring out the inaccuracy problem, especially for indices that only employ SWIR bands, i.e., MIRBI and dMIRBI.These areas were indicated by little change between moderate and high severity, with only a little light-coloured soil exposure reflectance decrease could be identified.The spectral signature of moderate and high severity areas that overlapped at the Visible and NIR band also could affect the accuracy of indices based on these bands.NDVI, dNDVI, EVI, and dEVI which employ these bands performed worse than NIR-SWIR spectral indices, i.e., dNBR, CSI, and dCSI, indicating the SWIR band outperformed the Visible band in burn severity detection, especially in discriminating moderate area and high severity area.Although the decrease of Visible reflectance and SWIR reflectance depended on the increase of black char, the different widths of SWIR band reflectance were more significant than Visible band reflectance.Kumar and Roy (2018) simulated the spectral response of different fires and found that fires with a larger size and higher temperature had higher reflectance values in the longer wavelength bands (e.g.SWIR) and negligible contribution in the Red band.Frazer et al. ( 2010) also mentioned that the NIR-SWIR bands enabled burned areas to be more easily discriminated against than the NIR-visible bands.
Visible-NIR indices, i.e.NDVI, dNDVI, EVI and dEVI showed poor performance.Since soils may contain organic materials and chemical constituents that could alter the reflectivity of the red band, using red band values does not model areas of exposed soil well.Also, burn severity and vegetation cover results may overlap, causing the findings to be overestimated (Wheeler, 2022).Hence the NDVI and dNDVI values appear not to be correlated to either burn severity and vegetation cover.Although, EVI as a modified version of NDVI that has been adapted for soil color and aerosol scattering and has improved sensitivity to areas which are rich in biomass (Ellsworth, 2012) could not outperformed NDVI due to the dry condition of the region.In addition, gray, which belongs to the blue spectral range, can be reflected from ash and char (Lewis et al., 2021) also could not be identified due to very little grey ash in the burned area.
Figure 7 depicts the increase of DN (Digital Number) values in band 10 Thermal Infrared 1 at the burned area.This increase was similar to the result of research conducted by Kafy et al. (2021), which showed an increase in Land Surface Temperature due to reducing vegetation cover.However, there were some problems using only this TIRS1 data for burn severity assessment in Gililawa Darat.First, the order of increasing DN values did not correspond directly with the increasing order of burn severity level.The DN values of the high severity area were not the highest; instead, those of the low severity area were the highest, and the DN value of the high severity area was only one level above the unburned area.Second, all burn severity classes overlapped, indicating poor accuracy of all severity classes.Consequently, if only Land Surface Temperature (LST) values were used to assess the burn severity at the burned area, there would be incorrectness in assigning high severity and low severity areas.Integrating LST with other burn severity such as EVI (Enhanced Vegetation Index) for burn severity assessment would probably produce better accuracy.
Green: Unburned, Yellow: Low, Orange: Moderate, Red: High All thermal indices had low Overall capacity and Kappa hat coefficient, indicating little agreement between burn severity and LST.The LST increase was detected in the burned area, but no increase linearity between LST and burn severity.Indeed, the decrease of vegetation in the burned area and the appearance of lower emissivity coverage (ash, char, and mineral soils) lead to a significant increase in the LST.However, aspect and illumination geometry also increase LST contrast in the different burn severity areas, so the better-illuminated slopes have higher contrast (Vlassova et al., 2014).Therefore, the nonlinearity of increase between LST and burn severity may be due to topographic influence, i.e., aspect and altitudinal differences.Thermal metrics are best to assess soil burn severity in locations with homogeneity in topography and moisture (Fernández-García et al, 2018), but they could not work well in Gililawa Darat that has topography heterogeneity.The previous research by Quintano et al. (2015) confirmed that topographic variables (elevation and aspect) related to LST contrast.Figure 6.5 depicts the relationship between burn severity, topographic variables, and LST.The more severe area in the high position, which faced north, tended to have higher LST.Indeed, integrating thermal with optical bands (LST/EVI and dLST/EVI) performed better than thermal metrics, but they did not exceed optical indices accuracy.This result was similar with the research result of Fernández-García et al. (2018).
This research result showed that the using of error matrix method alone to assess the agreement between burned area between reference and estimation (Salis et al., 2016;Giannaros et al., 2020) was not reliable for determining the correct fire ignition start coordinate.However, it could be used as a preliminary assessment to refine the low burned areas into the more potential fire origin area.
In this research, some simulations showed overestimation and underestimation of fire area.According to Rim et al. (2018) and Giannaros et al. (2020), the overestimated burned area of WRF-SFIRE output could be caused by overestimating wind speed due to the coarse atmospheric-fire model mesh.However, the overestimation in some simulations in this research was not caused by coarse mesh due to the overestimation did not occur evenly at all fire area.It seemed that the simulated fire spread shifted from the actual fire perimeter, therefore, apart from having an overestimated area at one flank, an underestimated area is also found at the opposite flank.Furthermore, Mandel et al. (2011) did not observe the overlarge fire spread at 400-metre resolution, whereas this research used a finer resolution (200metre) that could not possibly cause the overestimation of fire area.
On the other hand, the large underestimated burned area output could be attributed to an incorrect fire ignition start coordinate as mentioned by Giannaros et al. (2020) or to the incapacity of WRF-SFIRE to capture the extensive mid to long-range spotting.Simpsons et al. (2015) pointed out that on steep slopes, a wildfire spread model which was based on fuel properties, wind, and slope, but without a mid-to longrange spotting model, would understate the downwind extent of atypical lateral fire spread in a direction transverse to the background wind.During severe fires in rough terrain, vorticity-driven lateral spread (VLS) might occur.Wind and topography could interact to create lee slope eddies (rotors) that conveyed firebrands across the direction of the ambient wind, sparking many spot fires near the fire flanks.These spot fires subsequently merged with other spot fires and the fire edge, significantly increasing the rate of spread (ROS) of a fire flank (Hilton et al., 2019;Storey, 2021a).Spotting could accelerate dynamic spread caused by interactions between wind, steep slopes, and active fire, resulting in rapid lateral fire spread over ridges perpendicular to the primary wind direction.Most existing operational fire spread models had a severe weakness in capturing the influence of spots.(Storey, 2021b).Based on previous research by Simpson et al. (2016), lee fire whirls, which spread the fire laterally across the ridge, required steep lee-facing slopes, sufficient wind speed, and proper wind direction relative to the terrain aspect.Slope steepness at 20° required wind speed between 5 and 7.5 m s−1 and between 10 and 15 m s−1, 30° steepness needed wind speed between 7.5 and 10 m s−1, and steepness of 40° would allow pyrogenic vorticity, which developed and drove lateral spread if the wind speed was around 7.5 m s−1.Background wind is required to be angled between 20° and 30° relative to the terrain aspect.Because the wind direction was between 350° and 355° when the fire spread almost reached the east region, fire whirls would occur at lee ridge with terrain aspect between 10° -25°.Figure 6.11 shows the area where pyrogenic vorticity and fire whirls might occur when the fire spread reached (based on fire spread of valid fire ignition start coordinate).However, without the spotting model, the lateral spread was too slow to capture east region resulting in underestimation of the east flank of the fire area.
In addition, the simulation was also unable to capture the backing region of the fire area as close as the actual fire spread.This underestimation was most likely due to a simplified approach to backing fire ROS in the Rothermel model.The simplified analytical rate of spread (ROS) was used by empirical models to estimate the propagation of a fire as a function of time (Bakhshaii and Johnson, 2018).When the wind component normal to the fire line was present, the constant no-wind rate of spread was still used as a backing fire ROS (Coen, 2013).At the same time, there is no general backing fire ROS, and its evolution needs to be modelled with different spread mechanisms (Kochanski et al., 2013b).Kochanski et al. (2013a and2013b) and Peace (2014) identified that the Rothermel default no-wind ROS (which affects the flanks and backing fire spread) validates poorly for their study, and the same is likely to apply in the cases observed in this research.It could be understood that the model was mainly designed for wildfire risk management and firefighting (Mandel et al., 2014;Bakhshaii and Johnson, 2018); therefore, it focuses on the most active progression of the fire perimeter meanwhile backing fire which is the slowest spread among the fire's spread vectors (Coen et al., 2013) got less attention.
WRF-SFIRE could be customized and applied for simulating wildfire spread in Gililawa Darat island by providing the model with several datasets, i.e.National Digital Elevation Model (DEMNAS) from from DEMNAS BIG website at http: //tides.big.go.id/DEMNAS/, fire behavior fuel model in accordance with specific condition in the region using NDVI index to map the existing vegetations, and static geographical datasets and NCEP GDAS/FNL 0.25 Degree Global Tropospheric Analyses and Forecast Grids (ds083.3)meteorological datasets as initial and boundary conditions input from National Center for Atmospheric Research (NCAR).A series of WRF-SFIRE simulation using selected low burn severity coordinates of dNBR burn severity index map as fire ignition start coordinates showed that simulations from fire ignition start coordinates closed to origin area had the overall accuracy around 80% and the Kappa Hat Coefficient above 0,5 (moderate agreement).

Suggestion
A spotting model for fire flanks at steep slope must be developed in order to improve the functionality of WRF-Fire for covering the underestimation of the fire flank.An advanced research on backing fire spread can be conducted on improving model formulations in WRF-SFIRE in order to be able to approach the actual backing fire spread.The possibility of using the surface wind from geostationary satellites data with high temporal and spatial resolution, such as Himawari-8 data was need to be researched to substitute the wind derived from forecast or analysis meteorological data in order to approach the real condition
mistakenly assigned at the wrong location, i.e. the areas in the west and south east part of the island.The mixed burn severity indices (LST/EVI and d(LST/EVI)) showed the better distribution of burn severity levels but there were still large unburned areas that were identified at wrong locations.
CSI had a substantial level of agreement with Kappa Hat Coefficient of 0.667.MIRBI and dMIRBI had the same overall accuracy of 58.33% and Kappa Hat Coefficient of 0.667 0.444 indicating moderate agreement.NDVI, dNDVI, EVI and dEVI showed little accuracy with overall accuracy values ranged from 16.6% to 33.33% and Kappa Hat Coefficient values ranged from -0.111 to 0.111 indicating no agreement and slight agreement.Among the Landsat-8 OLI/TIRS thermal and mixed indices in Table 2, i.e.LST, dLST, LST/EVI, and d(LST/EVI), dLST showed the lowest overall accuracy with the value of 50% and Kappa Hat Coefficient of 0.333 indicating fair agreement.LST, LST/EVI, and d(LST/EVI) showed the same overall accuracy of 58.33% and Kappa Hat Coefficient of 0.444 indicating moderate agreement.

Figure 4 .
Figure 4.The dNBR Low Severity Area Coordinates Selected as WRF-SFIRE Fire Ignition Start Coordinates

Figure 8 .
Figure 8.The Region with Possible Pyrogenic Vorticity and Fire Whirls

Table 1 .
Accuracy Assessment of the Optical Burn Severity Indices Based on Sentinel-2B Datasets U: unburned, L: Low severity, M: Moderate severity, and H: High severity.

Table 2 .
Accuracy Assessment of the Thermal and Mixed Burn Severity Indices Based on Landsat-8 OLI/TIRS Datasets U: unburned, L: Low severity, M: Moderate severity, and H: High severity.

Table 3 .
Assessment of Burned and Unburned Points Coordinates Agreement of WRF-SFIRE Output with Those of dNBR