Testing NDVI, tree cover density and land cover type as fuel indicators in the wildfire spread capacity index (WSCI): case of Montenegro

This paper introduces the fuel properties into the wildfire spreading capacity index (WSCI), which originally relies on multi-criteria of social, hydro-meteorological, and geo-physical character of the context. Normalized Difference Vegetation Index (NDVI), Tree Cover Density (TCD), and land cover type are launched as indicators of fuel properties of the vegetated surfaces being indexed. The materials and software utilized here belong to a variety of open sources. CORINE Land Cover (CLC), Open Street Map (OSM), TCD via Copernicus high resolution data, and multispectral satellite images via Landsat 8 (Semi-Automatic Classification PluginSCP) are utilized as raw materials in a workflow in QGIS software. The study area is Montenegro, a Mediterranean territory prone to wildfire risk. Following the inventory stage, the indexing method relies on a normalizing procedure in QGIS and the assignment of weighted impact factor to each criterion via analytical hierarchy process (AHP). The WSCI value is derived as the sum of the products between the normalized class and the respective weighted impact factor of each criterion. The validation phase relies on ROC curve of WSCI values utilizing burned areas as positive cases. The results show high wildfire spreading capacities of vegetated areas of significant socio-economic and environmental profile such as the coastal zone along Adriatic Sea and the trans-boundary regions as part of Balkan Green Belt. Besides the methodological improvements the results of this work deliver tangible outputs in support of forest fire risk reduction means within disaster risk management and fire safety agendas in Montenegro.


Introduction Introduction Introduction Introduction
The projected future scenarios of climate change report significant increase of mean temperature values between 3.5 to 7 °C by the end of 21 st century in the south-eastern Europe (Lelieveld et al., 2012). For example, the projection being produced for the period between 2070 and 2099 foresee a decline of 70% in the precipitation records during the hot seasons in some areas of southern Europe. The forest fires are expected to become more frequent and further stretched in time as they are dependent on high temperatures and drought (Adger et al., 2007). Furthermore, the Mediterranean region historically is vulnerable to desertification, due to climate diversity and anthropogenic land consumption (Beguería et al., 2007). The correlation between wildfire events and both decreased precipitation levels and drought imply for a rising frequency and severity of forest fires. The risk is expected to be increased simultaneously in three directions; total number of days per year with fire risk, the elongated season of fire risk, and increase in the number of extreme events during the wildfire seasons.
The western Balkans are specifically mentioned for a remarkable wildfire risk rise in the mountainous areas (Moriondo et al., 2006). Thus, this study brings the case of Montenegro as a western Balkan country with a geography consisting of considerable vegetated mountainous areas. Analysis on the tendencies of daily temperature extremes between 1980-2010 has shown significant increase on both the annual number of summer days recording daily maximum temperature above 25 °C, and the annual number of tropical nights of nigh maximum temperature above 20•C in Montenegro (Burić et al., 2014).
The complexity of wildfire phenomenon is not a very recent finding, as its intermingled profile is previously stated in literature (Chapin et al., 2008). A complicated network of cause-effect cycles that are present in the spatial context underline its multi-layered character (Levin et al., 2016). Thus, simultaneous consideration of hydro-climatic, anthropogenic and vegetation factors in the landscape is helpful in a comprehensive understanding of forest fire regimes (Costa et al., 2011).
In this context, this study follows a multi-criteria approach as developed by Hysa and Baskaya (2019), which originally considers at the same time the social, environmental, and geophysical properties of the territory. Here, we bring the fuel properties along with the others as they are determinant to both wildfire ignition and spreading risk. The updated model is compared with the previous one which did not consider the fuel properties due to lack of data. The test is based on ROC curve analysis via SPSS software, where the positive cases are relying on burned surfaces. Corine land cover (CLC) data serves as the raw material for extracting the burned geometries.
In total there are 14 criteria shortlisted as having considerable implication with wildfire spreading capacity, the spatial information of which are provided through a diverse set of open source geospatial data. The method is both time efficient and cost-free, being reproducible to other study areas. The results of the study aim to serve as reference material to the institutions responsible of disaster risk reduction and fire management at both national and regional levels.

Study area
The territory of Montenegro is selected as the study area of this research due to its relatively small surface area, the presence of complex vegetation layers, and facing a considerable forest fire risk. Montenegro is a country within the western Balkan region in the south-eastern Europe ( Figure 1). The total surface area of the country is 13812 km 2 (Rajović and Bulatović, 2013). Despite its relatively small area, it is not rare to face a high level of diversity of geomorphology and vegetation type within a short distance within the territory (Wraber, 1983;Frankl et al., 2016).
Montenegro is reported to have faced a significant increase in the total area of the vegetated surfaces during the second half of the 20 th century. This is related to the post-50s migration from rural lands to urbanized centres being triggered by the process of industrialization (Nyssen et al., 2014). The distribution pattern of the vegetation layer is defined by (i) the distance to the Mediterranean Sea and coastal areas and (ii) the elevation above the sea level (regarding the temperate climate). According to the final report of "The First National Forest Inventory of Montenegro", 59.5% (about 826 000 ha) of the territory consists of forested surfaces and another 9.9 % (137 500 ha) of forestland. In other words, both classes cover almost 70% (964 262 ha) of the total territory of Montenegro. Referring to the same report, while the high forests consist of 51.1% (with 253.1 m /ha), the remaining 48.9% (with 62.6 m /ha) of the forest surfaces belong to the coppice forest cover (Dees et al., 2013). The predominant forest surfaces are portrayed by terminate canopy and by a considerable number of species of ground flora, shrubs and lower trees as shown in Figure 2 (Curovic et al., 2011). Andelic et al. (2012) report that within the forest fund of Montenegro there are registered 68 tree species in total. Among them 57 belong to the broadleaves and the remaining 11 species are conifers. Yet the existing forests include considerable degraded areas as well. The anthropogenic presence and social behavioural patterns have a great effect on forest degradation. Proximity to settlements and the adjacency to the transportation network and forest degradation has a strong correlation as common human practices like firewood harvesting are more frequent.

Forest fires in Montenegro
The western Balkans countries are highlighted as geographies that will struggle with forest fires soon if adequate preventive measures are not developed in time (Nemeth, 2015). According to the report of Forest Fire Country Studies of Montenegro the records kept by the national Forest Administration for the period 2003-2007, the average surface area of territories that faced forest fire was 4800 ha. While the mean number of fires was 53. According to the records between 2000 and 2015 in Montenegro there have been occurred at least 1000 large scale wildfires. The total estimated burned forest surface is approximately 15300 ha, and around a total mass of 500000 m 3 of timber were either lost or deteriorated. The Environmental Protection Agency of Montenegro (EPAM) reported that during the 2009 fire season the region which has been affected the most by forest fires is Bijelo Polje, while the least affected is Danilovgrad (Figure 3).
Within the National Strategy for Emergency Situations, fires are categorized as primary hazards at the national territorial scale. Furthermore, forest fires remain an important sub-category of this priority hazard type. Records show that considerable deterioration and material losses have been caused by forest fires in some complex forest habitats with unique properties or related areas with the National Parks and protected areas. The Black Forest area and Barno Lake basin Mlinski streams, Zabojsko zones around the lake, Crna Poda, Tara canyon and Susica river basin are among remarkable habitats which have been significantly affected by wildfires (Kovacevic, 2011).

Methodology and the workflow
This study relies on the method developed by Hysa and Baskaya (2019) for indexing broadleaved forest surfaces by their wildfire ignition probability index (WIPI) and wildfire spreading capacity index (WSCI). The original method did not include different types of vegetation or other criteria regarding fuel type due to data unavailability. Instead, in this study certain fuel properties are introduced as significant input in the spreading behaviour of wildfire. The procedure of the research builds on a workflow of five stages: preliminary, inventory, refinement, indexing, and validation.
I. First, the vegetated surfaces within the national borders are extracted from CORINE Land Cover data (CLC) as the core areas of analysis. A set of regular point's grid (500 m) which overlaps with the vegetated surfaces is targeted as the main layer of reference points. This set will serve as pivot points to be loaded the inventory measurements for each criterion. II.
An important step of the study is the inventory phase. It follows the inventory method developed by Hysa et al. (2017) which includes a procedure of recording the properties of the burned surfaces based on a variety of criteria. At this stage, the point's layer is calculated new field of attributes per each criterion via QGIS field calculator. III.
The loaded values of criteria are further normalized into a common range of values between 0 and 1 due to the diversity among individual range of absolute values. Furthermore, the refinement stage includes the assignment of a weighted impact factor for each criterion. The method used here is the analytical hierarchy process (AHP) pairwise comparison. AHP is used in multi-criteria assessment of natural hazards like flood and landslide as well (Stefanidis and Stathis, 2013). IV.
The main stage of the workflow is the indexing one. The WSCI values are calculated as a summation of the products between the normalized values of each criterion and the weighted impact factor. Thus, the WSCI value remains within a range between 0 and 1. The higher the WSCI value the higher the wildfire spreading capacity of a certain location within the vegetated surface. V.
The final stage consists of a tentative for validation of the method. The results of the study are compared with surfaces which according to CLC data of 2018 are reported as burned during recent years. Different from our previous work (Hysa and Teqja, 2020), at this stage we utilized the receiver operating characteristic (ROC) curve test. The ROC graph and the area under curve (AUC) value deliver quantifiable information about the reliability of the model. This analysis is performed in SPSS software.

Raw materials and data accessibility
The method used in this study requires a variety of geospatial data. They are selected as open access data and are derived via different data sources. CLC data which are available via European Environment Agency (EEA) Copernicus program, provides information about the vegetated surfaces, urbanized area, agricultural lands, and water surfaces. CLC data are provided at a spatial scale of 1:100000 where the minimum mapping unit is 25 ha for surfaces and 100m for linear elements in the landscape (Büttner, 2014). Hydro-meteorological information regarding solar radiation, precipitation, maximum temperature, and wind speed rely on raster data (resolution of 30 seconds) derived from WorldClim database (Fick and Hijmans, 2017). They are counted as monthly average values of July (between 1970-2000). Geophysical properties of the terrain such as slope, aspect, and altitude are provided by digital elevation model (DEM) via EEA.
Besides, the social, environmental, and physical aspects, in this study there is a need for fuel properties such as vegetation type, tree cover density (TCD), and normalized difference vegetation index (NDVI) which are introduced into the model as three new fuel criteria as shown in Figure 4. Vegetation type is informed by CLC data. Whereas, TCD is downloaded as a raster data from Copernicus portal (EEA).
Besides, NDVI data provide useful geospatial information about the land-surface properties including the vegetation layer. The relevancy of NDVI value to the wildfire regimes is already reported in literature (Leon et al., 2012). NDVI indexing which rely on the reflectance and absorbance capacity of vegetated surfaces is calculated as the ratio between the difference to the sum between the near-infrared (NIR) and red (RED) ranges of electromagnetic spectrum (EMS) (Chen et al., 2004).    The multi-spectral satellite images (Sentinel2) being available at 10 m resolution, are derived via the Semi-automatic classification plugin (SCP) in QGIS (Congedo, 2016;Leroux et al., 2018). Then, within the Sentinel2 layers the band 4 (RED) and band 8 (NIR) are used in the calculation of NDVI index. The layers are selected within the timeframe of the fire season and providing the minimum cloud coverage to maximize the reliability of the calculation. At this step, the raster images used here belongs to 25 th of July 2015 and have a total cloud cover less than 1%.

Multi-criteria inventory results
The first set of results in this study are derived from the inventory stage. As previously mentioned, at this stage the reference points are calculated the absolute measures regarding each criterion. Figure 5 shows the relative risk maps for each criterion represented by the respective measured value per each point location within the vegetated surface (among 40586 points in total). The absolute values are reclassified into ten classes via Jenks natural break optimization method, providing the minimum diversity within class and the maximum variance among ten classes (McMaster and McMaster, 2002).
Besides the series of maps presented in Figure 5, the inventory stage delivers numerical information as well. Table 1 presents the inventory values of each criterion reclassified into ten classes based on Jenks natural break as provided via layer properties symbology in QGIS. For example, the monthly precipitation records (average between 1970-2000) for July, show that the point location having the highest precipitation levels scores at 90 mm. While the median value is 71 mm, the least value remains at 48 mm.
The remarkable difference between the minimum and maximum values is evident while comparing the range of values among the listed criteria. The diversity in the range of values is inferable from the maximum values presented in the last column at Table 1. For example, while TCD (F2), precipitation (E2), and slope (P1) criteria are measured their respective top value at 99, 90, and 79.7, other criteria like NDVI (F3), altitude (P3), and wind speed (E4), scores at 0.88, 2.4, and 3.6 respectively. The diversity within the range of absolute values urges for a normalizing procedure before the indexing phase. While Jenks natural break is utilized just for visual representation of the results, the min-max normalizing technique is used here to reclassify each set of values per each criterion into a normalized range of values between 0 and 1 (Eq. 1). When the upper and lower bounds of the range are identifiable, this normalizing method is accepted to be the most appropriate normalizing alternative (Jain et al., 2005).
(1) Referring to the normalized values, the minimum and maximum values score 0 and 1 respectively. While the normalized median values considerably vary within the range between 0 and 1. The last rows of Table 1 show the WSCI index value as calculated via equation 2. The row of WSCI values range between 0.068 and 0.705, marking a median of 0.469. These results are represented both via their spatial coverage on the map of Montenegro (see Figure 6) and the histogram showing the frequency of the distribution (see Figure 7) of each WSCI value.
Besides the normalized inventory values, the other half of WSCI calculation relies on the weighted coefficient or the relative impact factor for each criterion. It is worth to restate here that this need relies on the assumption that not all criteria have equal impact on the wildfire spreading capacity of vegetated surfaces. As briefly explained in the methods section, the weights are assigned via AHP pairwise comparison method. AHP is applied in three steps utilizing the semi-automated tool developed by Goepel (2018). The pairwise comparison is based on literature review and brainstorm aiming to keep at minimum the consistency ratio (CR).
First, the main four categories (social, climatic-environment, physical, and fuel) are weighted in comparison with each other as presented in the 'IF (category)' row of Table 2. According to the results, the weighted value of fuel is 0.57. In other words, in means that 57% of the WSCI index value will be dependent on the criteria regarding fuel properties of the vegetated surfaces. While, social, climatic, and physical criteria score IF category values of 0.07, 0.24, and 0.12, respectively.  Figure 5. Figure 5. The territory of Montenegro including natural vegetation (a), and the relative risk map for each criterion; dist. to urban centers (b), dist. to main road (c), dist. to any road (d), solar radiation (e), precipitation (f), maximum temperature (g) wind speed (h), slope (i), aspect (j), elevation (k), distance to water (l), vegetation type-CLC (m), tree cover density (n), NDVI (o). Then, the criteria within each main category are weighted based on a pairwise comparison among each other, as shown in 'IF (within category)' row. Distance to main road (S4), wind speed (E4), distance to water (P4), and TCD (F2) score the highest IF within category values (0.67, 0.66, 0.57, and 0.54, respectively). The final weighted factor is calculated as the products of the former two, as shown in the bottom row 'IF (βj)' of the same table. Table 2 includes a further detail in the first row. It consists of relevancy factor for each criterion in reference with the wildfire spreading behaviour. Some criteria are directly related with the wildfire spreading capacity, while others are inversely related. For example, the larger the distance from urban centres, distance to road network, the slope, or the distance to water sources, the larger the wildfire spreading capacity of a certain location within the forest. In other words, the closer all these land uses the greater chances for a quicker firefighting mission and fire suppression.
On the other hand, some criteria are indirectly related with wildfire spreading behaviour. For instance, the larger the precipitation or the altitude values the lower the spreading chances of a wildfire at that specific location. The indexing stage consists of the field calculation procedure in QGIS using the WSCI equation (Eq. 2). The calculation of WSCI consists of the sum of the products between the normalized inventory measurement (Cj) and the weighted impact factor (βj) per each criterion: (2) where, βj is the weighted value of criterion j, and Cj is the normalized inventory value of criterion j.

WSCI indexing of vegetated surfaces in Montenegro
The indexing result of the workflow is presented in Figure 6. The map shows the final set of reference points being loaded the unique WSCI value according to the above-mentioned procedure. The index values are normalized and mapped after being reclassified into ten classes. The reclassification is following again the Jenks natural break, which is a useful optimization method in comparative studies since it visually highlights the differences among new classes. In the map, this difference is visualized via a colour palette of blue-yellowred gradient, where the red colour indicates high level of wildfire spreading capacities.
Referring to Figure 6, most of the reddish point are located either in the Adriatic coastal lands or along the south-eastern boundaries of the country. The former one is indicated mostly by the flammable vegetation type of a thick tree cover along the Mediterranean coast and the strong coastal winds. Both conditions, significantly contributes to the wildfire spreading capacities. While the second critical area is characterized by again dense tree cover layer being remotely located from the urban centres. While, these areas may have lower chances for forest fire ignition, they have significant wildfire spreading capacities if a forest fire occurs.
The reddish zones close to the south eastern and northern borders of Montenegro are crucial regarding the trans-boundary wildfire risk reduction. This is important as both neighbour countries, Albania and Bosnia-Herzegovina have considerable forested lands in the other side of their shared border with Montenegro. Moreover, the cross-boundary condition of wildfire risk is crucial in local scale as well, where the fragmentation in local administration systems and inefficiencies in risk planning must be aimed by neighbouring municipalities (Ager et al., 2018). According to Figure 6, there are critical zones shared between different local administrative regions in Montenegro. For instance, Kolasin, Andrrjevica, and Podgorica are sharing an area of relatively high levels of wildfire spreading capacities. While, Kolasin faces an opposite condition along the local administrative border shared with Savnik and Niksic. Figure 6. Figure 6. Figure 6. The frequency distribution of WSCI values among 40586 points locations is presented in Figure 7a. First, the references points layer (shp) is converted into a raster layer. Each pixel within the raster layer represents a single point and is loaded a single band value (WSCI index). The histogram chart is generated via the raster layer properties in QGIS. According to the histogram, the WSCI values being predominantly frequent are approximately at 0.56 (in normalized values it is 0.8), while a secondary peak is marked at 0.30 (normalized as 0.35). Figure 7a shows that considerable amount of point locations belongs to areas of high wildfire spread capacities.

Discussion Discussion Discussion
The validation stage relies on the comparative evaluation of the model results between the full set of points and the ones belonging to the burned surfaces. The information about the past fire events in this study is derived from CLC data, being documented under a specific land cover class entitled as 'burned area' (coded as 334). Since, one of the main goals of the study was to test the usefulness of introducing fuel properties (NDVI, TCD, and fuel type) into the WSCI, the validation procedure is split into two. WSCI values are simultaneously calculated as WSCI (including fuel properties) and WSCI_ESP (including just the environmental, social, and physical criteria). Then, for both cases, the WSCI results of the reference points which overlap with the burned surfaces (279 points out of 40586 in total) are compared to the values of the remaining ones, as previously shown in Figure 7b.
Moreover, the importance of cross-boundary fire risk among local administrative units necessitates analysis at municipality level as well. Figure 8 presents the box plot of WSCI and WSCI_ESP index values (not normalized) per each municipality in Montenegro. According to both box plots, the diversity among municipalities is higher in the case of WSCI_ESP values (Figure 8a). While the introduction of fuel properties into the model shifts the variance among municipalities towards a more stable state (see Figure 8b) especially regarding the lower and upper bounds of the range of values of each municipality.
Rozaje have the highest WSCI values (including the outliers in dots), based not only on the minimum and the maximum values, but also on the median and the set ranges of 3 rd and 4 th quartile. This fact is obvious in Figure 6, in which the territory of Rozaje is dominated by red colour. On the other hand, the territory of the capital municipality of Podgorica has the widest range of WSCI values in both plots, with no outliers. This is on the same line with the graphical information delivered by WSCI map (Figure 6), where the diversity among values is stretched into a gradient of red-yellow-blue palette within the borders of Podgorica. We tested further the introduction of fuel properties (WSCI) compared with the previous model (WSCI_ESP) by comparing their performance via ROC curve analysis as shown in Figure 9. The analysis is performed in SPSS. The records of burned surfaces are counted as positive cases of the test to calculate the sensitivity of the model. Figure 9 presents the result of the analysis in three means. First the ROC curve ( Figure  9a) shows that the introduction of fuel properties has decreased the sensitivity of the model. The overall model quality chart (Figure 9b) supports further this result, by indicating a quality coefficient of 0.55 in case of WSCI_ESP and 0.33 for WSCI. Furthermore, the area under the ROC curve (AUC) of WSCI_ESP is 63% larger than WSCI (Figure 9c), implying for a better performance of the model without including the fuel properties.
Even though the WSCI model sensitivity resulted to be low at national scale, this is not always the case at local scale. The sensitivity level of the updated model which includes the fuel properties, is diverse among different municipalities in Montenegro. Figure 10 presents a comparative ROC curve analysis of the WSCI model between the municipalities of Pluzine and Danilovgrad. According to ROC curve chart ( Figure 10a) the new model is significantly more sensitive for the territory of Pluzine in comparison to the area of Danilovgrad, scoring AUC values of 0.622 and 0.283 respectively (Figure 10b). While the overall model quality value is calculated as 0.57 for Pluzine, it stands at 0.01 for Danilovgrad (Figure 10c). Figure 9. Figure 9. Figure 9. Figure 9. Comparative values of ROC curve (a), overall model quality (b), and AUC between WSCI and WSCI_ESP (counting only the environmental, social and physical criteria). Figure 10. Figure 10. Figure 10. Figure 10. Comparative values of ROC curve (a), overall model quality (b), and AUC (c) between the municipality of Plizine and Danilovgrad At first look, the results of ROC curve analysis both at national and local levels can be considered contradictory with the objectives of the study or negative results. In other words, it can lead to a conclusion that the model is more successful without the inclusion of criteria about fuel properties. Yet, there are certain issues that must be discussed beyond the numerical results.
First, the proposed model is not aimed to be a mere predictive tool for the locations where future fire events may happen. But it is a tool for quantifying the fire spreading capacities after occurrence. So, among the vegetated surfaces which have been the target of this study, there may be areas of great wildfire spreading capacities. However, these areas may have not faced a fire occurrence yet. While other locations with a relatively lower WSCI value may have faced a fire ignition that have expanded into a wildfire burning at least a vegetated surface of 25ha (based on CLC evidences). This possible scenario implies that the location points within the burned surfaces do not have to score highest WSCI value since wildfire spreading is dependent to fire ignition. Thus, an integrated model that simultaneously considers the ignition and spreading criteria is crucial to be studied further.
Yet the model can be improved further. For example, the AHP is used to assign the relative impact factor for each criterion based on literature review and brainstorm. The later implies for a certain subjectivity which is reflected in the results as well. This can be considered a drawback especially when we count for the specificity of each study area in relation to driving forces of wildfire events. Thus, more rational method for assigning the relative impact factors must be developed. Future studies may apply the ROC curve analysis workflow in the opposite direction. Consequently, ROC curve analysis may assist in assigning the relative impact factors for each criterion in WSCI model. The AUC values of each criterion can be normalized as relative impact factor in WSCI calculation.

Conclusions Conclusions Conclusions
This paper presented an improved version of the GIS based method developed by Hysa and Baskaya (2019) for indexing the forest surfaces by their wildfire spreading capacity (WSCI). Relying on the original set of criteria, this study introduced three new factors regarding the fuel properties of the vegetated surfaces. Vegetation type in reference with flammability, the tree cover density, and the NDVI values were introduced as evidence about the composition of the fuel. The new criteria are informed by geospatial raw data being accessible as open source.
The process was stretched into five stages including the preliminary work, inventory measurements, refinement (normalizing and weighting via AHP), indexing via WSCI equation, and the validation based on CLC evidences and ROC curve analysis method. The workflow was modelled in Graphical Modeler in QGIS software enabling considerable reduction in processing time and easier reproducibility of the method to other study areas.
A crucial stage of the process is the validation procedure which at this point relies on a single evidence such as the burned areas being reported within CLC data. The reliability of the method could be extensively enhanced by comparing the model results with the historical data about the wildfire occurrences and their intensity within the study area. ROC curve analysis resulted suitable as a validation mean. Yet, further improvements of the method were suggested in the discussion part.
According to the results presented in this study, it can be concluded that the vegetated surfaces within Montenegro that are located along the coastal zone and along the south-eastern boundaries hold the highest wildfire spreading capacities. The territories along the Adriatic coast are critical as the majority of the touristic and year around social activities are located. An increased wildfire spreading risk is crucial to be properly managed during the pre-occurrence phase in these areas as the vulnerability is the highest. Similarly, the second hotspot being highlighted by the results of this study has a clear overlap with the protected natural reserves of both national and international importance. The trans-boundary forest surface shared with Albania in the south are crucial part of the extended Balkan Green Belt (Geidezis and Kreutz, 2012).
Finally, the results of the study are useful in support of disaster risk management and preventive measurements especially during the pre-occurrence phase. The WSCI map is helpful in highlighting hotspots within the territory of Montenegro that must be carefully managed by the respective institutions at both local and national levels.

Authors' Contributions Authors' Contributions Authors' Contributions Authors' Contributions
Conceptualization: AH; Methodology: AH; Software: AH; Calculations: AH; Literature: AH and VS; Data curation, AH; Writing-original draft preparation, AH; Writing-review and editing: AH and VS; Visualization: AH; funding acquisition, AH and VS.
The authors read and approved the final text of the manuscript.