Zagros Grass Index – a New Vegetation Index to Enhance Fire Fuel Mapping: Case Study in Zagros Mountains

1. Introduction

Forest fires pose significant concerns for the environment, economy, and human safety in a majority of forested regions globally [1,2,3,4]. Ecologically, fire serves as a pivotal factor influencing vegetation diversity and dynamics over time and space [1,2,3,4,5]. Authorities, including civil protection agencies, governments, local authorities, and forestry corps, are compelled to effectively manage forest fires and establish preparedness strategies to preserve biome services and ensure citizen safety [6,7,8,9]. While forest fires can arise naturally due to dry weather, volcanic eruptions, or lightning, human activities stand out as the predominant factor, particularly during periods of heightened water stress [10,11].

The forest areas and rangelands along the western and northern expanse of the Zagros Mountains (Iraq) chain have faced a considerable number of fires since 2005. The forests in Marivan and Paveh, situated in the Kurdistan and Kermanshah provinces of western Iran, respectively, are particularly affected [12,13]. Furthermore, the Province of Sulaymaniyah and Halabja in the Kurdistan Region (KR) of northern Iraq has also witnessed a notable surge in fire incidents in recent years, especially after 2008 [10]. According to both official reports and conducted studies, human activities are identified as the most frequent ignition sources in these forested areas [13,14]. It is also reported that more than 90% of fires in the European Union (EU) are human-caused [15,16]. Given the valuable opportunity presented by satellite-based indices for monitoring diverse Earth phenomena, Remote Sensing (RS) data and Geographic Information System (GIS) technology have become crucial tools for natural resource managers and researchers across government agencies, conservation organizations, and industry [17,18,19]. The integration of RS data/techniques and GIS facilitates the efficient and accurate analysis of wildfire dynamics, enabling informed decision-making processes for fire management and mitigation strategies [17,18]. The occurrence of fires, including their severity and duration, is intricately correlated to vegetation condition [10,17,18], including critical dynamic factors such as Fuel Moisture Content (FMC) and Fuel Temperature (FT) [20,21]. Several studies have been conducted worldwide to classify and map land cover due its wide and important role in natural resources management [22], agriculture management [23], biodiversity conservation [24], among others. Amongst them, a range of researchers have tried to differentiate and map vegetation species such as trees, shrubs, and grass species using RS data and techniques. Some studies have used Light Detection and Ranging (LiDAR) data with trees height information to differentiate these species since these vegetation types have different height [23]. However, LiDAR data is not available everywhere and it is expensive [25]. Deep learning (DL) algorithms, on the other hand, has also been used for land cover classification in numerous recent studies [26]. Accordingly, advanced DL techniques with high-resolution satellite data have better performance than traditional methods for classifying land cover and detecting objects. Although Dl algorithms achieved high accuracy, it needs more diverse training data to be efficient in different situations. However, most of them are applied on high resolution satellite images [26]. In the work of Saah et al. [22], DeepLabV3+, a semantic segmentation-based DL method, was employed to categorize three types of vegetation land covers (trees, shrubs, and grass) utilizing solely Sentinel-2 RGB images.

In contrast, there are other methods that rely on Vegetation Indices (VIs) that are generated from multispectral images [27]. VIs allows to extract valuable information from the spectral characteristics of plants, encompassing biochemical characteristics, environmental factors, and soil properties [27]. These indices play a crucial role in estimating vegetation biomass, canopy height, percentage vegetation cover changes, plant health, and Leaf Area Index (LAI). Additionally, they aid in distinguishing between soil and vegetation and mitigating atmospheric and topographical influences when feasible [28]. Meshesha et al. [29] found a strong correlation between forage biomass and spectral indices by employing Sentinel-2 Normalized Difference Vegetation Index (NDVI) and Enhanced Vegetation Index (EVI), integrated with ground sampling in Harshin district, Ethiopia, to develop a forage forecasting model. In Fakhri et al. [30], a novel vegetation index-based workflow was introduced, and, within it, the multi-objective particle swarm optimization (MOPSO) algorithm was applied to optimize a set of broadband VIs to reach both objectives of greenness estimation and vegetation/non-vegetation classification in a small area in Zagros sparse woodlands. A new index was also developed by Qian et al. [31] and applied in Beijing, China, merging spectral and texture features to differentiate tree from grass in urban areas, at a detailed level, using high resolution GeoEye-1 imagery. Another study conducted in northwest (NW) Russia utilized hyperspectral data and vegetation phenology to differentiate tree species [32]. It was concluded that classification using multispectral data effectively improves accuracy compared to using a single hyperspectral image. In another study area in Zagros [33], a study was conducted to generate accurate land cover map for the Shirvan County forests, a part of Zagros forests in Western Iran using Sentinel-2 derived NDVI, Google Earth imagery, and field data for protective management. The study proved that the Support Vector Machine (SVM) algorithm had the highest accuracy for the classification of Sentinel-2 data with an overall accuracy of 81.33%. In a greater area, for the entire Zagros Mountains, a new empirical model was introduced for mapping land cover for the whole Zagros mountains using Sentinel-2 derived NDVI [34]. Despite being so challenging, this study has effectively mapped the land cover (agriculture, build-up area, wooded area, plantation, bare soil, water, and rangeland) for its wide study area. In Iraq, on the other hand, some studies were also conducted for mapping land cover in KR using VIs and DL approaches [35,36]. Although, several studies have employed satellite images for mapping land cover, a minority of them are dealing with differentiating grass species from wood species, and if they do, they are focusing on urban areas using high resolution images [22,23,24,25,26,27,28,29,30,31,32,33,34,35,36]. However, the grass species have not been disregarded completely, and it is mostly classified as a member of rangeland class with other members such as shrubs. The land cover, vegetation dynamic, topography, distance from population center and road, besides many other static and dynamic factors, are considered as key factors in the fire susceptibility assessment [37,38,39]. Deriving reliable information on fuel types is regarded as a major factor since the fires need fuel to happen and to propagate [37,38]. Actually, in studying fuel types, grass species have not been looked only as an ecological factor, but as the most flammable fuel type [37,38]. Notably, the NDVI has been widely employed to estimate vegetation phenology, as well as its quality and growth condition [40,41,42]. NDVI, serving as an index of vegetation growth and coverage, finds extensive use in describing spatio-temporal characteristics of land use and land cover (LULC), including percent vegetation coverage [41,42,43,44]. However, the NDVI cannot differentiate between tree, shrub, and grass because of their similar spectral characteristics [22,37]. Using multi-temporal NDVI integrated with phenological information of vegetation covers is effectively helpful to differentiate vegetation cover [27,31,32]. Regarding this characteristic, the NDVI has been widely used in fire susceptibility studies to represent vegetation dynamic and condition [37,38,39].

In this research, we introduce a novel index, the Zagros Grass Index (ZGI) [45]. The ZGI was developed utilizing the Moderate Resolution Imaging Spectroradiometer (MODIS) NDVI to specifically identify dry grass masses, known to be highly susceptible to fires. This was achieved by integrating the phenological traits of both forest trees and grass species found in Marivan and Sarvabad, located in the western part of Iran. The primary aim of the ZGI is to serve as an easily accessible supplemental tool for detecting and mapping dry grass mass. Designed as an accessible tool for mapping dry grass, the ZGI's utility is further explored in this study by extending the analysis across the Zagros Mountain forests in western Iran and northeastern Iraq from 2013 to 2022. The goal is to assess the ZGI’s applicability, scalability, and generalizability to broader areas, enhancing fire management strategies in these regions.

2. Materials and Methods

2.1. The Study Area

The study area is a part of Kordu-Zagrosian (KZ) forests of Iran and Iraq (longitude:$44°27\mathrm{\text{'}}37"\mathrm{E}-46°52\mathrm{\text{'}}20"\mathrm{E}$; and latitude:$34°30\mathrm{\text{'}}55"\mathrm{N}-36°33\mathrm{\text{'}}5.827"\mathrm{N}$), in the northern Zagros Mountain chain with a mean elevation of 1287 m Above Mean Sea Level (AMSL) (Figure 1). The KZ region has been used in some studies and corresponds to a geographical and ecological zone encompassing parts of the Zagros Mountain range across western Iran, eastern Iraq, and southeastern Turkey [46,47].

In Iran, the study area encompasses the forest of Marivan and Sarvabad in Kurdistan Provinve, as well as Paveh, Javanrod, Ravansar, and Salas in Kermanshah Province, which are in the west of Iran (Figure 1). The study area also covers the vast area of Sulaymaniyah and Halabja provinces, in KR, northern Iraq. Most of the Iranian part of the study area (almost 90%) is located in 1000 m AMSL, while nearly 35% of the Iraq study area is higher than 1000 m AMSL, and the rest range from almost 200 m to 1000 m AMSL. The KZ forests are dominantly covered by Brant's oaks (Quercus brantii) species, in coppice and thin trunks [48], which resulted in the proximity of trees canopy and grass species (Figure 2). These forests are also mostly open canopy expand mostly between 750 m to 1700 m AMSL [46]. On the other hand, the Grass Species (GS), which are the most vulnerable and flammable fuel type [45,49], extend all over the forests, but in different density, type, and growth pattern. The majority of these species dry in early summer (Figure 2) [10], except for a little species (sub-alpine vegetation) which are mostly in very high elevations (over 2000 m AMSL) [46] where trees don’t grow and forest fires have not been occurred there.

According to official reports, annual rainfall of Sulaymaniyah and Halabja ranges between 375- 724 mm with a semi-arid continental weather regime. It means that it is cold and wet in winter and very hot and dry in summer [50]. Significantly, the summer months from June to September are very hot and dry, except for the mountainous areas (e.g., Hawraman, Qandil, Penjwen, etc.). In July and August, the hottest months, mean temperature are 39º- 43º Celsius and often reach nearly 50º Celsius [51]. The eastern part of the study area, located in Iran, has more rainfall, with a mean annual precipitation of 700-991 mm with a standard deviation of 200 mm [52].

2.2. Data Sources

In this study, we employed two datasets, satellite imagery and field data, to examine vegetation dynamics (Table 1).

The primary satellite dataset used was the 250 m MODIS vegetation product (MOD13Q1) 16-day time series, covering the period between 2013 and 2022. This data provided access to a wider period, regarding our study data set, while other satellite images (e.g., 10 m Sentinel-2 images) cannot provide it, despite their higher spatial resolution. Furthermore, the Shuttle Radar Topography Mission’s Digital Elevation Model (SRTM DEM) was also used [53]. The SRTM DEM is free to download and provides a spatial resolution of 30 m [53]. The second dataset, provided by the Department of Natural Resources and Watershed (DNRWK) of Kurdistan and Kermanshah Provinces, offers in situ information on fire incidents [54]. Another data such as the border of the study area, towns, water bodies, and cultivated areas were provided by the municipalities of the targeted areas.

2.3. Zagros Grass Index

In this study, we aim to employ our previously established methodology [45] to detect and recognize grass-covered surfaces, a primary fuel source for fires, in the expanded study areas which are referred as Kurdo-Zagrosian Forests (KZF). Expanding our methodology, we also refine the proposed index regarding the elevation of the different parts of the study area, considering the effect of elevation on vegetation phenology. Distinct spectral reflectance emitted by vegetation throughout different seasons is effectively characterized by phenological metrics such as the Start of Season (SOS), End of Season (EOS), and Maximum of the Season (MAX) [55]. Figure 3 presents the flowchart that is followed to create ZGI maps. In this context, RS-based VIs, particularly the NDVI, have proven to be valuable tools for vegetation detection and extraction of phenological metrics [56]. Historically, the MODIS and the Advanced Very High-Resolution Radiometer (AVHRR) products have been extensively used for estimating these phenological metrics [56,57,58,59]. The SRTM DEM was additionally employed to delineate the study area based on elevation, as it significantly influences the phenology of both tree and grass species.

Consequently, the elevation affects SOS, MAX, and EOS of both tree and grass species [60,61]. Accordingly, the study area was classified into two major classes: i) the area over 1000 m AMSL, and, ii) the area under 1000 m AMSL. Figure 4 illustrates the variation in NDVI changes between areas above 1000 m AMSL and those below 1000 m AMSL.

Accordingly, MAX and EOS are different for these two classes, and they happen sooner for the lower areas (almost 15 days). Furthermore, the NDVI range for these two classes is not the same and the overall NDVI is higher for the area over 1000 m AMSL, due to the higher volume of forest area in this area, while the area under 1000 m AMSL are mostly rangelands.

Apart from the effect of elevation on phenological metrics, the tree and grass species in both elevation regions have different phenological behavior. Figure 5 depicts the overall NDVI changes over time in KZF for either forest areas, encompassing both species, and rangeland, predominantly covered by grass species.

Regarding the dependency of phenological metrics on both elevation and vegetation species, the proposed phenological scenario is illustrated in Figure 6.

Accordingly, Grass-SOS (GSOS) and Tree-SOS (TSOS) are almost the same for both elevation range, while the Grass-MAX (GMAX) and Tree-MAX (TMAX) for higher elevation areas is later than lower elevation areas. The EOS, on the other hand, differs either from trees to grass, or from higher area to lower area. Grass-EOS (GEOS) start in the Mid-June for higher areas, while it starts from late of May for lower areas. Tree-EOS (TEOS) is almost the same for the both elevation range. The SOS and EOS don’t affect our subject.

According to the phenological scenario (Figure 6), there are life and green Tree Species (TS) and dry and dead Grass Species (GS) after Mid-June. Therefore, the greenness after June is only from TS. Consequently, the positive NDVI values which display green areas are only from TS and the NDVI values from grass species are almost zero since they have dried.

From a mathematical perspective, the subtraction of GESO’s NDVI from GMAX’s NDVI can effectively emphasize the presence of grass. This subtraction delineates the decline in NDVI from the TMAX and GMAX to the GEOS, thereby elucidating the extent of grass coverage. It is important to note that this decline solely pertains to non-TS, as the trees retain their greenness. Accordingly, the proposed index will be defined as equations 1 and 2:

For open forest area in each year:

$${\mathbf{Z}\mathbf{G}\mathbf{I}}^{\mathit{o}\mathit{p}\mathit{e}\mathit{n}}={\mathit{N}\mathit{D}\mathit{V}\mathit{I}}^{\mathbf{G}\mathbf{M}\mathbf{A}\mathbf{X}}-{\mathit{N}\mathit{D}\mathit{V}\mathit{I}}^{\mathbf{G}\mathbf{E}\mathbf{O}\mathbf{S}}$$

(1)

For dense forest area in each year:

$${\mathit{Z}\mathit{G}\mathit{I}}^{\mathit{d}\mathit{e}\mathit{n}\mathit{s}\mathit{e}}=\mathrm{Average}\mathrm{value}\mathrm{of}\mathrm{ZGIopen}=\frac{1}{\mathit{N}}{\sum }_{\mathit{n}=1}^{\mathit{N}}{\mathit{Z}\mathit{G}\mathit{I}}_{\mathit{n}}^{\mathit{o}\mathit{p}\mathit{e}\mathit{n}}$$

(2)

Where, ${NDVI}^{\mathrm{G}\mathrm{M}\mathrm{A}\mathrm{X}}$ corresponds to the NDVI value of the GMAX which is the date within the NDVI gives the maximum value due to the attendance of both TS and GS. GMAX is 129^th day for the areas over 1000 m, and 113^th day for the area under 1000 m AMSL. ${NDVI}^{\mathrm{G}\mathrm{E}\mathrm{O}\mathrm{S}}$ corresponds to the NDVI of the GEOS which corresponds to the date within the NDVI gives its minimum value while the TS are still green, but GS have dried [47]. ${NDVI}^{\mathrm{G}\mathrm{E}\mathrm{O}\mathrm{S}}$ is the NDVI of 257th day of a year in this study. It can also be any other date after mid-June and before early October. ${ZGI}^{open}$ is the ZGI of open forests, and ${ZGI}^{dense}$ is the ZGI of dense forests pixels, which is calculated by averaging the ZGI of neighbor open forest pixels. N is the number of all pixels which are recognized as open forests in a year, and n is the pixel counter.

The ZGI was applied to verified natural surfaces (forests and rangelands) using a designated mask, which excluded the artificial areas (urban, water bodies, rural, and cultivated regions) from the study area. The mask has been created from pre-provided maps prepared by DNRWK, then updated manually using QGIS software base maps. Artificial areas are named non-forest areas and represented as black areas in the resulted maps. The areas affected by the fires have been maximized and repositioned within the margins of the maps to offer a more detailed representation of the ZGI status in the burned areas and their surrounding regions.

3. Results

This The procedural steps resulted in distinct ZGI maps for each year (Figure 7), with the fire occurrence locations provided by the administration of forest and watershed of Sulaymaniyah, Iraq, and from Kurdistan and Kermanshah provinces in Iran. Most of the data don’t encompass the exact locations of fire incidences, but rough information of affected areas has been provided, except for Marivan’s data which include a coarse Universal Transverse Mercator (UTM) coordinate (E, N) of fire occurrences. The black areas on the maps denote artificial regions. Besides, a number of reported fires have happened within agriculture fields which are observable within Non-Forested areas (Figure 7b3).

4. Discussion

The resulted maps revealed distinct patterns in the ZGI for each year, consistent with previous studies, that underscore the significant influence of climate conditions, such as rainfall and temperature, on the vegetation cover in semi-arid and semi-Mediterranean regions [30,34,62,63]. The ZGI thresholds have been defined using Natural Breaks (NB) method, which is based on actual values in the dataset rather than using predefined intervals [64]. Since the fires are mostly human-cause [13,14], deliberately or accidentally, the affected areas are majorly adjacent to artificial area where human are living or cultivating [39,65]. Despite the remarkable results of previous studies on land cover mapping, within the study area and regions with comparable climate conditions, particularly emphasizing grass species as ecological contributors [22,23,24,25,26,27,28,29,30,31,32,33,34,35,36], this study is focusing on grass as a potential source of dangerous dead biomass and the most flammable fuel type in our study area and other similar regions [37,38,39]. Spread and continuity of dead or living vegetation are primary factors in sustaining fires, where dry grasses - especially dense ones - has been proven to be the highest potential for fire propagation [37,68]. Consequently, the identification and mapping of fuel types plays a crucial role in defining risk conditions [38]. Fuel models involve the parameterization of different fuel types to estimate their fire behavior [38,69]. Various methodologies have been developed to generate and map fuel types according to the input data, intended use, and the scale of the study [20,37,68,69,70]. Despite numerous global studies characterizing different fuel types through detailed parameters, such as crown height, crown base height, vegetation coverage percentage, forest canopy density, crown density, canopy bulk density, number of trees per area, vertical and horizontal continuity, moisture content, live and dead fuel load, biomass [70], and using RS data/methods, this study focusing on dry grass mass [38,69]. It leverages phenological characteristics of grass and woody species within the study area. Furthermore, the KZ’s coppiced TS are more susceptible to ignition by fires originating from grasses, compared to conifers (e.g. pines) and other TS which have been addressed and categorized in other studies [38,48].

The resulting maps prove that the fires strongly overlap and follow the areas with higher ZGI values, specifically 0.24<ZGI<0.35 and 0.35<ZGI<0.46. The higher ZGI range (ZGI>0.46) belongs to mountainous areas (over 2000 m AMSL), mostly cover by sub-alpine vegetation species, and rarely experience fire [47].

Although, in some years with lower overall ZGI values, lower fires have been reported (e.g. Figure 7A, 7H and 7I) it cannot be inferred that poorer dry grass mass necessarily results into less fire, as seen in Figure 7C, since the overall ZGI doesn’t mean that dry grass mass is necessarily low all over the study area. Another reason can be the coppice structures of majority of forest trees that make them vulnerable even against fires ignited from poor grass condition (Figure 2) [71]. Apparently, the severity of fires should also be regarded to have a more comprehensive understanding on fire behavior. Reliable information on dry grass mass in fire season (late May- Late October) can also give a reliable view about the propagation rate and severity of the fires [37,72]. Although the fire severity and propagation rate have not been addressed directly in this study, ZGI can also help us with those because it can be inferable that the severity and propagation of fires could be potentially very high among dense forests with highly flammable grass mass (high ZGI values) [37,72]. Moreover, it is imperative to consider the various grass species, considering their distinct phenology, density, propagation, ignitibility, power of ignition, and firing duration.

The maps further indicate that within areas situated above 1000 m AMSL, the ZGI values are higher compared to those in lower elevations [10,30,33,34]. This observation aligns with the overall fire’s distribution across the study area, as depicted in Figure 7, affirming that the fire tends to be higher in regions above 1000 meters AMSL. It is also observed from the maps that the areas close to human associated areas (black areas) with lower ZGI (brown and green areas) have been rarely subjected to fires. Therefore, integration of human accessibility and ZGI, rather than higher ZGI values alone, emerges as a significant factor in fire incidences. Comparing the frequent application of NDVI and in fire susceptibility, as a vegetation detector in Zagros forests [10,30,33,34,39], ZGI seems to be a more suitable index for the study area since it is targeting the dry mass through pixel-wise temporal change detection of NDVI regarding the phenological characteristics of the grass and non-grass covers.

The Zagros Forest region experiences a semi-Mediterranean, and continental semi-arid climate regime, characterized by hot, dry summers and mild, wet winters [73]. Regarding the vast distribution of this type of climate regime on the world, this index can be generalized and can be used in other study areas, e.g., California (USA), Mediterranean Basin, and Southwestern Australia [74].

5. Conclusions

This Forest structures and the composition of grass species, in conjunction with human activities, emerge as significant contributors to the incidence of wildfires. Effectively mapping of grass species enhances our comprehension of fire susceptibility, their duration, and spatial distribution. The integration of RS data/techniques and GIS expedites a precise analysis of wildfire dynamics, providing reliable insights into the dynamic and static conditions of vegetation species. Moreover, phenological characteristics can be harnessed to formulate novel RS indices, exemplified by the ZGI applied in this study. The resulting maps exhibit a pronounced overlap between ZGI and areas accessible to human activities, suggesting a notable correlation between human-caused fires and regions with high ZGI values. The ZGI can be used jointly with other VIs (e.g. NDVI), or as an alternative index for providing helpful vegetation dynamic, flammable area, and fuel load among Zagros Mountains in Iran and Iraq.

Regarding the coarse resolution of MODIS images, the ZGI is only based on the overall phenological and spectral characteristics of trees and grass covers, ignoring texture differences, unlike urban areas utilizing high-resolution images. Therefore, for enhanced accuracy, future research may benefit from utilizing satellite images with higher spatial resolution (e.g. Sentinel-2, Landsat-8/9, GeoEye). The pursuit of further studies using satellite imagery with higher spatial resolution can yield more reliable results, offering improved insights into fire dynamics and vegetation characteristics at finer scales. The higher ZGI value, consequently, shows the drier mass. Nevertheless, this criterion alone is insufficient for predicting fires; multiple other factors (e.g. temperature, aspect, slope, distance from population center and road, etc.) must be considered as contributors to fire susceptibility.

However, this study lacks precise details such as coordinate system, areas affected, fire types, and causes. Some years did not have corresponding data on fires as well. For example, we didn’t access any recorded data related to the years before 2015, as well as, for 2021 for Kermanshah provinces cities (Paveh, Jwanro, Ravansar, and Salas). We also couldn’t find any data about Marivan and Sarvabad fires for 2018. Lack of awareness or attention from relevant administrations in obtaining accurate data on fire incidences, including field observations on fire severity, a dependable map of burned areas, and precise ignition points, has complicated the execution of our study. Regarding the semi-Mediterranean and continental semi-arid climate of the study area, characterized by hot, dry summers and mild, wet winters, this index can be generalized to other study areas, e.g., California (USA), Mediterranean Basin, and Southwestern Australia.