Quantifying Forest Cover Loss as a Response to Drought and Dieback of Norway Spruce and Evaluating Sensitivity of Various Vegetation Indices Using Remote Sensing

: The Norway spruce is one of the most important tree species in Europe. This tree species has been put under considerable pressure due to the ongoing impacts of climate change. Meanwhile, frequent droughts and pest outbreaks are reported as the main reason for its dieback, resulting in severe forest cover loss. Such was the case with Norway spruce forests within the Kopaonik National Park (NP) in Serbia. This study aims to quantify, spatially and temporally, forest cover loss and to evaluate the sensitivity of various vegetation indices (VIs) in detecting drought-induced response and predicting the dieback of Norway spruce due to long-lasting drought effects in the Kopaonik NP. For this purpose, we downloaded and processed a large number of Landsat 7 (ETM+), Landsat 8 (OLI), and Sentinel 2 (MSI) satellite imagery acquired from 2009 to 2022. Our results revealed that forest cover loss was mainly driven by severe drought in 2011 and 2012, which was later significantly influenced by bark beetle outbreaks. Furthermore, various VIs proved to be very useful in monitoring and predicting forest health status. In summary, the drought-induced response detected using various VIs provides valuable insights into the dynamics of forest cover change, with implications for monitoring and conservation efforts of Norway spruce forests in the Kopaonik NP.


Introduction
In Europe (excluding Russia), from 1990 to 2020, according to GFRA 2020 (Global Forest Resources Assessment), the area under forest increased by about 17 million ha, while in certain areas of Europe, the process of deforestation has been more evident in recent years [1].In this context, due to various disturbances, such as adverse climate events, insect pests, epiphytocia, etc., one of the most endangered species in Europe is Norway spruce (Picea abies (L.) Karst.)[2][3][4][5][6][7][8].This is very important as Norway spruce, from an ecological and economic point of view, is one of Europe's most important conifer species.It can be found in various areas, occurring from the boreal forests in Northern Europe to subalpine areas of the Alps and the Carpathian Mountains [9].In the Balkans, unlike the rest of Europe, Norway spruce rarely occurs at lower altitudes, and it is usually found in high mountains and mountain ranges where the climate is cold and humid [10].In due to long-lasting drought effects in the Kopaonik NP.Based on our research goals, we formulated two hypotheses:

Hypothesis 1 (H1).
There is a significant relationship between drought-induced dieback of Norway spruce and forest cover loss.
Hypothesis 2 (H2).Different VIs showed varying levels of sensitivity in the prediction of droughtinduced dieback of Norway spruce, with some VIs providing more accurate prediction of forest cover loss than others.

Study Area
The study area (Figure 1) is situated within the Kopaonik National Park (NP) in southern Serbia, which gained its current status in the year 1981 due to its biodiversity, rich flora and fauna, and great cultural and historical importance [69,70].The park stretches across 11,969.04ha of land in low-populated areas of municipalities Brus and Raška, mainly on mountain Kopaonik (2017 m.a.s.l) [71,72], of which 7427.24 ha is covered by forests [73].Higher parts of the mountain are mainly covered with pure or mixed conifer stands of Norway spruce (Picea abies (L.) Karst.) and Silver fir (Abies alba (L.) Mill.), with or without European beech (Fagus sylvatica L.), which in addition to Austrian pine (Pinus nigra Arn.) and oak species (Quercus spp.), dominates the lower parts of the mountain [73].Such species distribution is mainly driven by the wide altitudinal range, namely by site-specific ecological conditions of different altitude levels.Generally, the climate in the Kopaonik NP is characterized as subalpine [70], with an average annual temperature of 4.1 • C and an average annual precipitation of 1040.1 mm (climatic sequence 1991-2020) [74].By comparing the last two climatic sequences (1961-1990 and 1991-2020), it can be found that the average annual temperature in the Kopaonik NP increased by 1.4 • C, and the average annual precipitation increased by 119.3 mm [75].As significant devitalization and dieback of trees are reported more frequently in pure stands and less in mixed stands of Norway spruce, we narrowed the research area down to 2385.72 ha of such forests, using forestry stand maps provided by the Kopaonik NP.A major component of the research area is located in the area under the protection regime of the second degree, where, according to Ðor dević et al. [73], limited and strictly controlled use of natural resources and activities is established to the extent that it does not endanger natural habitats.

Data Collection
To evaluate the impact of drought on the forest cover loss at Mt. Kopaonik (Apendix A), we downloaded Landsat 7 (ETM+), Landsat 8 (OLI) Level 1, and Sentinel-2A/2B (MSI) Level 1C satellite imagery (from 2009 to 2022) using the U.S. Geological Survey Earth Ex-

Data Collection
To evaluate the impact of drought on the forest cover loss at Mt. Kopaonik (Appendix A), we downloaded Landsat 7 (ETM+), Landsat 8 (OLI) Level 1, and Sentinel-2A/2B (MSI) Level 1C satellite imagery (from 2009 to 2022) using the U.S. Geological Survey Earth Explorer website (https://earthexplorer.usgs.gov,accessed 11 January 2024) and the Semi-Automatic Classification v.7.10.11-Matera(SCP) plugin [76] from the QGIS v.3.22.6 Białowie ża (OSGeo, Chicago, IL, US) software (Tables 1 and 2).The 2009 to 2022 time period was selected to ensure that the state of vegetation in pre-drought (2009), drought (2011 and 2012), and post-drought (2013-2022) periods when severe pest outbreaks occurred was analyzed in order to obtain a complete picture of how Norway spruce is responding to the adverse effects of climate change.We selected only the cloud-free imagery acquired during the growing season, which, in our case, included imagery acquired only in July and August (except for one image from June).The 2010 imagery was not downloaded because, in all available Landsat 7 (ETM+) data, the images covering most of our research area were covered with clouds.

Data Processing
The downloaded Landsat 7 (ETM+) and Landsat 8 (OLI) MS bands, R, G, B, NIR, SWIR1, and SWIR2, including Sentinel-2 (MSI) Level-1C MS bands, B, G, R, VRE, VRE2, VRE3, NIR, NIR2, SWIR2, and SWIR3, were automatically processed using the SCP plugin by converting them from DN [Landsat] and scaled top of atmosphere (TOA) reflectance [Sentinel] into the TOA reflectance to reduce the inter-scene variability through a normalization for solar irradiance.Atmospheric correction of all images was carried out using an image-based technique called Dark Object Subtraction (DOS1) [77], as cited in [76].Ordinary least squares regression (OLS) equations from Roy et al. [78] were used to normalize the reflectance of one Landsat sensor to the other (ETM+ to OLI).Before applying the pan-sharpening Brovey Transform technique [79] using the SCP plugin, as recommended by Rahaman et al. [80], we calculated individual relationships of Landsat 7 (ETM+) and Landsat 8 (OLI) R, G, B, and NIR bands with the PAN band using regression analysis with R Studio v.4.3.2 (Posit, PBC, Vienna, Austria) [81] and a raster [82] package.The results showed weak relationships between those variables for several Landsat 7 (ETM+) MS bands, R: r 2 -0.44,G: 0.62, and B: 0.41, except for NIR: r 2 -0.90.Because of the possible distortion of spectral data that might occur after pan-sharpening these MS bands, which may produce misleading conclusions in time series analysis of vegetation indices (VIs), we only used original MS Landsat 7 (ETM+) bands.On the contrary, Landsat 8 (OLI) bands showed a strong relationship with the PAN bands R: r 2 -0.99,G: 0.99, and B: 0.99, except for NIR: r 2 -0.54.As such, we used pan-sharpened Landsat 8 (OLI) MS bands (R, G, and B) in forest cover loss analysis for the years 2013 and 2014.

Forest Cover Loss Analysis
The land cover classification was carried out using the Supervised (semi-automatic) classification, which involves identifying materials in the image according to their spectral signatures by drawing the Regions of Interest (ROIs-Training Areas) over the homogeneous area of an image.For the sake of precise drawing, we used high-resolution imagery of the year 2022, provided in Google Earth Pro v.7.3.6.9345-r0(Google, Mountain View, CA, USA), overlaid with different MS band composites of downloaded imagery.Of all tested MS band composites, the so-called "agricultural composite" (SWIR1-NIR-B) and the "short-wave infrared composite" (SWIR2-SWIR-R) performed best in underling the difference between stands dominated by conifer or deciduous trees.In this way, we excluded stands dominated by deciduous trees from our analysis.Finally, we drew eleven reliable and constant ROIs for all years analyzed, six for forest cover (average area 6.17 ha) and five for non-forest cover (average area 7.24 ha), which were evenly distributed all over the area.Forest cover included all canopy undisturbed stands, while non-forest cover included forest glades, meadows, bare lands, and small artificial objects.After drawing all the ROIs, they were dissolved to form two land cover macro classes.Using the Land Cover Signature (LCS) classification in the SCP plugin [76], we defined spectral thresholds for each ROI signature (a minimum value and a maximum value of each MS band), defying the spectral region of each land cover macro class.Spectral thresholds were calculated for all years separately to avoid misclassification of land cover due to inter-year variability in the vegetation spectral characteristics.Pixels that were not classified in either of the two macro classes, that is, pixels found inside overlapping regions or outside any spectral region, were classified using the Minimum Distance algorithm [76,83].In this way, Euclidean distance was calculated between the spectral signatures of every pixel in the image and ROI spectral signatures, thus assigning each pixel to the class of the spectral signature that was closest.After the land cover classification, the final raster processing was conducted using the Postprocessing group of tools in the SCP plugin, which included, to a certain extent, the correction of incorrectly classified pixels and the merging of rasterized polylines and polygons of roads and other artificial objects into classification rasters, whose incorrect classification may contribute to the misinterpretation of the results.Using the Accuracy function in the SCP plugin, the accuracy assessment of the produced maps (classification rasters) was performed with the calculation of an error (confusion) matrix by comparing produced map information with reference data [84], which was, in our case, high-resolution imagery provided in Google Earth Pro v.7.3.6.9345-r0(Google, Mountain View, CA, USA) (CNES/Airbus, Maxar Technologies, etc.).Given that each produced map contained more than 100,000 pixels checking the classification accuracy of all of them would be impractical from several points of view.Therefore, a stratified random sampling method was used for this research.The total sample number was calculated for each analyzed year separately (from 2013 to 2022) by applying Equation (1) [85,86]: where n is the number of samples (ROIs), S(Ô) is the standard error of the estimated Overall Accuracy that we would like to achieve (here used as 0.01), W i is the mapped proportion of the area of map class, and S i is the standard deviation of stratum (values proposed by Olofsson et al. [86]).Sample size allocation to strata (map classes) of each analyzed year was calculated as an average number of proportional and equal sample size allocations previously calculated for each stratum.The random distribution of samples for each map class was conducted using the SCP tool Multiple ROI creation (to create stratified random points).The process of labeling (assigning) the sample units to each macro class was carried out using Google Earth Pro v.7.3.6.9345-r0(Google, Mountain View, CA, USA) and upon its completion, the data were exported into KMZ format, which was finally converted into the shapefile (.shp) format to match SCP Accuracy tool requirements for the calculation of accuracy quantitative measures, such as Error Matrix, Overall Accuracy (OA), Producer's Accuracy (PA), and User's Accuracy (UA) [86].Forest cover loss was calculated as the absolute and relative difference between the surface area of forest cover (ha) in the reference year (2013) and all other years consecutively.The cumulative forest cover loss dynamics were calculated on a fragment level, as an average area change of all of them, excluding non-forest areas existent in 2013.Land cover classification results visualization was conducted using R Studio v.4.3.2 (Posit, PBC, Vienna, Austria) [81] and a raster [82] package, and the sf [87], RColorBrewer [88], ggplot2 [89], ggpmisc [90], patchwork [91], and gt [92] packages.

Evaluation of VI Sensitivity in Detecting and Predicting Drought Effects in Norway Spruce Forests
To examine the state of forest health and vitality pre-drought and during the drought period (2009-2014) that preceded forest cover loss, we selected multiple VIs from different groups, such as Typical VIs, Water VIs, and wetness and greenness components of the Tasseled Cap (TC) transformation (Table 3).
The selection of VIs was based on their sensitivity in detecting various vegetation properties.For example, Typical VIs are well known for assessing photosynthetic activity, forest health status, and detecting forest stressors such as pest outbreaks [43,51,55,[93][94][95][96].On the other hand, Water VIs primarily provide a quantitative measure of water content in various tree species, early detection of water stress, and assessment of drought impacts on forested areas [48,51,52,57,[96][97][98][99].Tasseled Cap (TC) transformation components are selected as they compress multispectral data into a few bands associated with physical scene characteristics with minimal information loss [100], thus sharing or having greater sensitivity in detecting various vegetation properties of both Typical VIs and Water VIs [42,101,102].
A calculation of mean values was segregated on the spatial level to areas where forest cover loss occurred and to areas where it did not.The spatial distribution of those areas was taken from the land cover classification (LCC) rasters.For this purpose, we selected the years 2015 (when the bark beetle outbreak started) and 2017 (when the bark beetle outbreak reached its peak), excluding non-forest areas that were present in the LCC raster from the year 2014.By this means, data preparation was made for an analysis whose only purpose was to determine if there was an association between the spatial distribution of forest and forest cover loss and variation in VI values.Therefore, we used Cohen's d [108] to measure the effectiveness of the VIs in forest cover loss detection by determining whether or not there is a statistically significant difference between VI values in the areas of forest and non-forest cover (forest cover loss), and how large that difference is Equation (2).Calculation of Cohen's d was carried out using the R Studio v.4.3.2 (Posit, PBC, Vienna, Austria) [81] and lsr [109] packages.The effect size was classified using the Sawilowsky scale [110], where 0.1 represents a very small effect size, 0.2 a small effect size, 0.5 a medium effect size, 0.8 a large effect size, 1.2 a very large, and 2.0 a huge effect size.For this research, we only considered a medium, large, very large, and huge effect size sufficient to predict forest cover loss.As such, according to Cohen's U3 [108], 69.1%, 79.8%, 88.5%, and 97.7%, respectively, of the lower-meaned land cover class areas are exceeded by the average VI value in the higher-meaned land cover class area: where X 1 is the mean of the first group, X 2 is the mean of the second group, SD 2 1 is the standard deviation of the first group, and SD 2 2 is the standard deviation of the second group.

Forest Cover Loss
Concerning land cover change, the analysis showed (Figure 2) that forest cover loss in the Kopaonik NP rose continuously until the end of the observed period.Until 2017, forest cover loss nearly doubled each year (Figure 3).From 2017 to 2019, forest cover loss increased roughly by 0.5%, and from 2019 to 2020 by 0.14%.The 2021 forest cover loss was identical to that of 2020, while in 2022, forest cover loss increased by 0.28%.In total, 5.75% of forest area was lost on behalf of the non-forest cover (without Norway spruce trees, and not counting the saplings).
In 2014 and 2015, forest cover loss occurred in relatively small areas throughout the analyzed area (Table 4 and Figure 3).In those years, the maximum area of forest cover loss in terms of fragment level-later referred to as forest glades-did not exceed 1 ha, and this area was quite uniform, as indicated by a low standard deviation.Later, in 2016 and 2017, forest glades emerged in new areas, while existent areas spread out, mainly in the west and northwest, which was particularly evident in 2017 (Figure 3).The mean area of forest glades was almost identical to those from 2014 and 2015, although the maximum area  4).Moreover, the spread of existing forest glades continued in the western and northwestern parts of the research area (Table 5), which was confirmed by the increased standard deviation.The 2018 results are particularly interesting, as the maximum area of forest glades and standard deviation nearly doubled compared to the previous year (Table 4).The spread of existing forest glades continued in the mentioned areas (Figure 3), while its area increase, as indicated by mean, maximum, and standard deviation values, was the largest in the observed period (Table 4).Interestingly, the results of 2019, 2020, 2021, and 2022 do not significantly differ from 2018 (Table 4), as only a small number of forest glades emerged across the analyzed area, including a small-scale continuous spread of existing glades in the west, northwest, and northeast, which gradually increased until the end of the observed period (Figure 3, Table 4).
by the average VI value in the higher-meaned land cover class area: where  ̅ 1 is the mean of the first group,  ̅ 2 is the mean of the second group,  1 2 is the standard deviation of the first group, and  2 2 is the standard deviation of the second group.

Forest Cover Loss
Concerning land cover change, the analysis showed (Figure 2) that forest cover loss in the Kopaonik NP rose continuously until the end of the observed period.Until 2017, forest cover loss nearly doubled each year (Figure 3).From 2017 to 2019, forest cover loss increased roughly by 0.5%, and from 2019 to 2020 by 0.14%.The 2021 forest cover loss was identical to that of 2020, while in 2022, forest cover loss increased by 0.28%.In total, 5.75% of forest area was lost on behalf of the non-forest cover (without Norway spruce trees, and not counting the saplings).In 2014 and 2015, forest cover loss occurred in relatively small areas throughout the analyzed area (Table 4 and Figure 3).In those years, the maximum area of forest cover loss in terms of fragment level-later referred to as forest glades-did not exceed 1 ha, and this area was quite uniform, as indicated by a low standard deviation.Later, in 2016 Quantitative accuracy assessments showed that the OA was quite high for all years, indicating that the produced maps have the required accuracy (Table 5).However, this should be taken with caution, as the dominance of forest cover pixels (Table 6) influences the calculation of the OA.Qualitative accuracy assessment revealed some classification  Quantitative accuracy assessments showed that the OA was quite high for all years, indicating that the produced maps have the required accuracy (Table 5).However, this should be taken with caution, as the dominance of forest cover pixels (Table 6) influences the calculation of the OA.Qualitative accuracy assessment revealed some classification errors in those years, such as misclassification of non-forest cover as forest cover and, rarely, vice-versa.Such errors occurred in areas where individual trees were scattered all across non-forest patches and where the area of the non-forest cover was inconsiderably small (less than 100 m 2 ).It is worth noticing that the PA for non-forest cover was significantly lower than for forest cover (Table 5).Such results indicate that numerous non-forest areas were misclassified as forest cover, especially in the years 2017, 2020, and 2021 (Table 6), which was also pointed out by the qualitative assessment of LCC in those years.

Evaluation of VI Sensitivity in Detecting Responses to Drought and Predicting the Dieback of Norway Spruce
The results of the VI time series analysis (Figure 4) showed that a large number of VIs used to reflect the Norway spruce response to the severe drought that occurred in 2012, such as the DSWI, NDVI, TCG, TCW, EVI, TVI, and SAVI.
Regardless of drought, the values of other VIs, such as the MSI, NDMI, and NMDI, remained almost constant across all years analyzed.On the other hand, for all VIs used, responses to drought occurring during 2011 were non-existent.As for the post-drought period (2013-2014), the values of VIs regained the pre-drought state.
The analysis showed that the DSWI, NDVI, and TVI in the pre-drought period (2009) did not predict the spatial distribution of forest and non-forest cover (forest cover loss) in 2015 (d < 0.50) (Figure 5).On the other hand, the MSI, NDMI, and NMDI showed a medium effect in predicting the spatial distribution of forest and non-forest cover (forest cover loss) in 2015.As for the EVI, SAVI, TCG, and TCW methods showed a large effect (d≈0.80) in predicting the spatial distribution of forest and non-forest cover (forest cover loss) in 2015.
For the drought period (2011-2012), we also found a large effect for the MSI, NDMI, NMDI, and TCW in predicting forest and non-forest cover (forest cover loss) in 2015, where the strength of such effect slightly increased in 2012 (Figure 5).Moreover, a medium effect in predicting the spatial distribution of forest and non-forest cover (forest cover loss) in 2015 was shown using the DSWI, EVI, SAVI, and TCG in 2011, while, in 2012, only the SAVI and TCG maintained their predicting ability.
The analysis results of the post-drought period (2013-2014) are quite interesting, as contrary to the drought period, the MSI and NDMI now have a medium effect, while the EVI, SAVI, TCG, and TCW have a large effect in predicting the spatial distribution of forest and non-forest cover (forest cover loss) in 2015.The strength of such an effect for the above-mentioned VIs increased to a very large effect (d > 1.20) and reached its peak in 2014.Again, in the same year, the MSI and NDMI showed a large effect in predicting the spatial distribution of forest and non-forest cover (forest cover loss) in 2015.
VIs used to reflect the Norway spruce response to the severe drought that occurred in 2012, such as the DSWI, NDVI, TCG, TCW, EVI, TVI, and SAVI.
Regardless of drought, the values of other VIs, such as the MSI, NDMI, and NMDI, remained almost constant across all years analyzed.On the other hand, for all VIs used, responses to drought occurring during 2011 were non-existent.As for the post-drought period (2013-2014), the values of VIs regained the pre-drought state.The analysis showed that the DSWI, NDVI, and TVI in the pre-drought period (2009) did not predict the spatial distribution of forest and non-forest cover (forest cover loss) in 2015 (d < 0.50) (Figure 5).On the other hand, the MSI, NDMI, and NMDI showed a medium effect in predicting the spatial distribution of forest and non-forest cover (forest cover loss) in 2015.As for the EVI, SAVI, TCG, and TCW methods showed a large effect (d≈0.80) in predicting the spatial distribution of forest and non-forest cover (forest cover loss) in 2015.
For the drought period (2011-2012), we also found a large effect for the MSI, NDMI, NMDI, and TCW in predicting forest and non-forest cover (forest cover loss) in 2015, where the strength of such effect slightly increased in 2012 (Figure 5).Moreover, a medium effect in predicting the spatial distribution of forest and non-forest cover (forest cover loss) in 2015 was shown using the DSWI, EVI, SAVI, and TCG in 2011, while, in 2012, only the SAVI and TCG maintained their predicting ability.
The analysis results of the post-drought period (2013-2014) are quite interesting, as contrary to the drought period, the MSI and NDMI now have a medium effect, while the EVI, SAVI, TCG, and TCW have a large effect in predicting the spatial distribution of  Further analysis showed, in contrast to the previous one, that the NDVI, EVI, SAVI, TVI, and TCG in the pre-drought (2009) and drought (2011-2012) periods did not predict the spatial distribution of forest and non-forest cover (forest cover loss) in 2017 (d < 0.50) (Figure 6).Exceptions were the MSI, NDMI, DSWI, NMDI, and TCW in 2009 and 2011, which showed a medium effect in predicting the spatial distribution of forest and nonforest cover (forest cover loss) in 2017.Starting from 2012 (the year of severe drought), we found a large effect of the MSI, NDMI, NMDI, and TCW (except for DSWI) in predicting the spatial distribution of forest and non-forest cover (forest cover loss) in 2017 (Figure 6).The strength of such a relationship remained constant in 2013, a post-drought year, except for the NMDI, where it slightly increased.In the same year, the NDVI, EVI, SAVI, and TVI, which previously showed no effect in the pre-drought (2009) and drought (2011-2012) periods (d < 0.50), reached a medium effect in predicting the spatial distribution of forest and non-forest cover (forest cover loss) in 2017.
Interestingly, in 2014, another post-drought year, the EVI, SAVI, and TCG showed no effect, while the NDVI and TVI retained a medium effect in predicting the spatial dis- It is important to note that the NDVI and TVI, through all years analyzed (2009-2014), did not predict the spatial distribution of forest and non-forest cover (forest cover loss) in 2015.Still, the DSWI in 2014 showed a medium effect in predicting the spatial distribution of forest and non-forest cover (forest cover loss) in 2015.
Further analysis showed, in contrast to the previous one, that the NDVI, EVI, SAVI, TVI, and TCG in the pre-drought (2009) and drought (2011-2012) periods did not predict the spatial distribution of forest and non-forest cover (forest cover loss) in 2017 (d < 0.50) (Figure 6).Exceptions were the MSI, NDMI, DSWI, NMDI, and TCW in 2009 and 2011, which showed a medium effect in predicting the spatial distribution of forest and non-forest cover (forest cover loss) in 2017.Starting from 2012 (the year of severe drought), we found a large effect of the MSI, NDMI, NMDI, and TCW (except for DSWI) in predicting the spatial distribution of forest and non-forest cover (forest cover loss) in 2017 (Figure 6).The strength of such a relationship remained constant in 2013, a post-drought year, except for the NMDI, where it slightly increased.In the same year, the NDVI, EVI, SAVI, and TVI, which previously showed no effect in the pre-drought (2009) and drought (2011-2012) periods (d < 0.50), reached a medium effect in predicting the spatial distribution of forest and non-forest cover (forest cover loss) in 2017.
(Figure 6).Exceptions were the MSI, NDMI, DSWI, NMDI, and TCW in 2009 and 2011, which showed a medium effect in predicting the spatial distribution of forest and nonforest cover (forest cover loss) in 2017.Starting from 2012 (the year of severe drought), we found a large effect of the MSI, NDMI, NMDI, and TCW (except for DSWI) in predicting the spatial distribution of forest and non-forest cover (forest cover loss) in 2017 (Figure 6).The strength of such a relationship remained constant in 2013, a postdrought year, except for the NMDI, where it slightly i Figure 6.Evaluation of VI sensitivity in predicting forest cover loss in 2017 (gray-no effect, green to red gradient color scale-medium to very large effect).
Interestingly, in 2014, another post-drought year, the EVI, SAVI, and TCG showed no effect, while the NDVI and TVI retained a medium effect in predicting the spatial distribution of forest and non-forest cover (forest cover loss) in 2017.It is important to note that the DSWI had a large effect, and the MSI, NMDI, and TCW had a very large effect in predicting the spatial distribution of forest and non-forest cover (forest cover loss) in 2017.

Forest Cover Loss
In the example of the Kopaonik NP, it can be seen from the results of this study that Landsat 8 (OLI) and Sentinel 2A/2B (MSI) satellite imagery can be used, with satisfactory accuracy, in the mapping of small forest cover losses.Moreover, the high UA for non-forest cover (Table 5) also indicates satisfactory accuracy, as most pixels classified as non-forest cover represent the real state in the field.Nevertheless, both quantitative and qualitative accuracy assessments showed some minor drawbacks.For example, the lower PA for non-forest cover (Table 5) may indicate the impossibility of correctly classifying areas smaller than 10 × 10 or 15 × 15 m due to spatial resolution limitations of both sensors used (Sentinel 2 MSI up to 10 m and Landsat 8 OLI up to 15 m).Such was the case with KC et al.'s [111] land cover classification of Rupandehi District, Nepal, where barren land was classified as neighboring water bodies due to its small size.Sometimes, in an area of one pixel, we can find many different types of land cover, which significantly alter pixel spectral signature; thus, in the classification process, pixels can be assigned to the wrong land cover class.Inter-seasonal variation in vegetation photosynthetic activity and the current health status of forest cover may also alter its spectral signature, for example, to be similar to the neighboring non-forest cover (grassland or underbrush).This was observed by Forsythe et al. [56] to be the main reason for lower PA values in some classification results.The combination of both events surely contributed to the classification errors.Nevertheless, such errors can be ignored, as the undetected loss of several trees does not represent a significant error from the forestry management point of view.
Considering that 5.75% of the pure Norway spruce forest in the Kopaonik NP ceased to exist in the post-drought period (Figure 2), it is hard to attribute such a state exclusively to the drought effects.Kesić et al. [28] came to the same conclusion, claiming that soil acidification and monodominance of Norway spruce at Mt. Kopaonik were other possible reasons for its dieback.However, the nearly double increase in forest cover loss during 2015 and 2016 (Figure 2) can be easily attributed to the effects of pest outbreaks.As reported by Matović et al. [31] and Stojanović et al. [33], in those years, there was a huge outbreak of I. typographus and P. chalcographus, which, at that moment, acted like a primary pest.However, it should be taken into account that bark beetle outbreaks in Norway spruce forests are a consequence of adverse climatic effects, such as drought, as their defensive mechanisms are weakened when affected by summer drought [112][113][114][115]. Spatial-temporal expansion of forest glades in 2015, 2016, and 2017, which previously emerged over small areas in 2014 (Figure 3 and Table 4), clearly indicate bark beetle activity.Such a trend continued in later years with less intensity, following a decline in bark beetle outbreaks.However, a few questions arise.Is the forest cover loss a result of a single factor or the interaction of several factors?Are particular stands more susceptible to drought than others?The to these questions should be sought through the implementation of various long-term multidisciplinary research projects in these forests.

Evaluation of VI Sensitivity in Detecting Responses to Drought and Predicting the Dieback of Norway Spruce
Although the NDVI, EVI, TVI, SAVI, TCG, DSWI, and TCW revealed a large-scale drop in vegetation vigor and canopy water content all over the analyzed area, that is, the response of Norway spruce to severe drought occurred in 2012 (Figure 4), not all VIs predicted forest cover loss in 2015 (Figure 5).Besides TCW, Cohen's d showed that other VIs, which did not show any response of Norway spruce to severe drought in 2012 (MSI, NDMI, and NMDI), had large and very large effects in predicting forest cover loss in 2015.A similar result was found for 2011, which was a year with less severe drought occurrence.Although the MSI [46] and NDMI [49] are considered to be highly effective in assessments of moisture stress in plants, this was not the case in our study.Based on such results, we can assume that NIR-SWIR1 ratio-based Water Vis, such as the MSI and NDMI, indicated only different soil water retaining capacities in areas where forest cover loss occurred and where it did not.We found the base for this assumption in a conclusion in Welikhe et al.'s research [104], where it was reported that MSI is strongly correlated to soil moisture at 20 cm depth.On the other hand, in a review study, Le et al. [49] summarized findings from other studies [98,116,117], concluding that the NDWI method (in our research named NDMI) yielded unsatisfactory results when applied to forest objects for water stress monitoring.Worth noticing is the large effect of the pre-drought (2009) results of the EVI, SAVI, TCG, and TCW in predicting the forest cover loss in 2015, as such a state points to pre-drought differences, and possibly the susceptibility of different Norway spruce populations, or their respective habitats, to drought events in the Kopaonik NP.The cause of this may be found in the research of Rehschuh et al. [118], in which they reported that Norway spruce trees growing on shallow, well-drained soil expressed a relatively higher drought sensitivity compared to trees from a site with deep, silty soil.The practically non-existent ability to predict the forest cover loss in 2015, with the post-drought data (2013 and 2014) using the NDVI and its modified version TVI, should not be considered unusual.Although these VIs showed strong sensitivity in the detection of Norway spruce response to severe drought, they cannot be used in predicting forest cover loss, as they do not exhibit any statistically significant difference between VI values in the area of forest and non-forest cover (forest cover loss).As such, we agree with Le et al.'s [49] conclusion stating that the NDVI cannot be effectively used in the early detection of drought effects.On the contrary, other "drought-sensitive" VIs, such as the EVI, SAVI, TCG, and TCW, showed a large (2013) to very large effect (2014) in predicting forest cover loss in 2015, indicating that the post-drought period is crucial in predicting drought effects, as it can strongly suggest where forest cover loss might occur.In contrast, these VIs, except for the TCW, did not perform well in predicting forest cover loss in 2017 (Figure 6), indicating that the primary cause of Norway spruce dieback after 2015 was mainly driven by pest outbreaks.As seen in Figure 2, forest cover loss doubled from 2015 to 2017.Such a finding goes in line with an earlier report from Matović et al. [31], where it was stated that, in those years, bark beetle began to act as a primary pest.What challenges this conclusion is a post-drought medium (2013) to a large effect (2014) of the DSWI and a large (2013) to a very large effect (2014) of the TCW in predicting forest cover loss in 2017 (Figure 6), which may indicate a direct influence of drought on the loss of forest cover in 2017.Nevertheless, so-called Water VIs (MSI, NDMI, and NMDI) performed almost the same as for 2015 forest cover loss prediction-having a large (2012) to very large effect (2013 and 2014) in predicting forest cover loss.Considering these results together with previous conclusions, where we stated that such results only indicated different soil water retaining capacities in areas where forest cover loss occurred and where it did not, we can only confirm such assumptions.

Implication for Conservation of Norway Spruce Stands in the Kopaonik NP
As indicated by the results, severe drought greatly impacts forest cover loss in Norway spruce stands in the Kopaonik NP.Although severe drought has not occurred since 2012, according to Miletić et al. [37], such events may occur more often in the future.Accordingly, we can only expect that forest cover loss will continue to rise.However, we did not take into account several other reasons, which surely had or may have a great impact on forest cover loss.In their study in the Kopaonik NP, Matović et al. [31] found that devitalization and dieback of Norway spruce trees were more pronounced in structurally and age-homogeneous stands.As such, within areas of protective regimes, it should be legally enabled to implement adequate forest management measures that will support structural and age differentiation.Furthermore, the introduction of complementary species, such as Silver fir (Abies alba Mill.) and European beech (Fagus sylvatica L.), to improve stability and overall resistance of Norway spruce stands should not be neglected, as the monodominance of one species, such as Norway spruce, leads to instability and reduced tolerance to pests and adverse climatic events, as was proven in our and many other studies [119][120][121].Regarding the deforested areas, support should be provided through tree planting.As Tanovski et al. [116] proposed, this should involve using reproductive material of known origin with adaptive properties suitable for the environmental conditions of the regeneration site.

Methodological Limitations of the Used Methodology in Detecting Responses to Drought and Predicting the Dieback of Norway Spruce
The main reason for some previously proven VIs, such as the NDWI, MSI, and NDMI [38,116], exhibiting low performance in monitoring and predicting the health status of Norway spruce forests in the Kopaonik NP may be the lower spatial or spectral resolution of the imagery used.Although, from 2015 until 2022, higher spatial and spectral resolution Sentinel 2 (MSI) imagery was used for land cover mapping, strong sensitivity in predicting forest cover loss using lower spatial and spectral resolution Landsat 7 (ETM+) and Landsat 8 (OLI) imagery was simply impossible due to various factors.For example, one pixel in Landsat 7 (ETM+) and Landsat 8 (OLI) imagery may have mixed spectral values, as it, in a spatial manner, contains up to three pixels from Sentinel 2 (MSI) imagery, which may include distinct land cover types.A similar problem was reported in Abdollahnejad et al.'s [42] study, which points out that lower-resolution satellite imagery has limited use; that is, it could be used only in studies where sample sizes are not less than the spatial resolution of used imagery.Taking into account that Sentinel 2 (MSI) has been in orbit since June 23, 2015, such shortcomings, in the context of this study, could not be overcome.Another problem lies in the low temporal resolution and unavailability of cloud-free Landsat 7 (ETM+) and Landsat 8 (OLI) imagery during the entire growing season in Kopaonik NP.If more were available, coupled with ground-measured meteorological data, it would be easier to determine which drought levels or their cumulative effects, along the growing season, trigger Norway spruce dieback in the future.On the other hand, the usage of very-high spatial and spectral resolution imagery, such as Pléiades 1A/1B, QuickBird, SPOT 6/7, WorldView-2, etc., could provide precise and clear answers even at the single tree level.However, their high cost was a limiting factor in the framework of this study.Considering that other factors, such as stand and terrain characteristics, play a significant role in Norway's spruce dieback [28,31,118], future analyses should include these factors.This can be achieved by employing machine learning methods to provide more accurate and reliable results.

Conclusions
Based on the obtained results, remote sensing proved to be a powerful tool in analyzing the spatial-temporal trends of forest cover loss and the sensitivity of VIs in detecting responses to drought and predicting the dieback of Norway spruce in the Kopaonik NP.Our analysis revealed a consistent increase in forest cover loss in the Kopaonik NP until the end of the observed period (2013-2022), which was driven by severe drought in 2011, especially in 2012.Regarding rates of loss, forest cover loss nearly doubled each year, reaching its peak in 2018.Such a state was significantly influenced by emerging bark beetle outbreaks, which acted as a primary pest from 2015 to 2017.Overall, approximately 5.75% of forest cover was lost on behalf of non-forest cover.It was also found that forest cover loss first started in 2014 and 2015, occurring in small areas, but later expanded, particularly in the western and northwestern parts of our study area.
Further analysis revealed that several VIs, including the DSWI, NDVI, EVI, TVI, SAVI, TCG, and TCW, proved to be sensitive in detecting the response of Norway spruce to severe drought in 2012.However, the values of other VIs remained relatively constant across all years, with no observable indication of the response of Norway spruce to drought.The predictive capabilities of VIs for forest cover loss varied across different periods, with some VIs showing a significant effect during drought and post-drought periods (EVI, SAVI, TCG, and TCW).Notably, the NDVI, TVI, and DSWI were unable to predict forest cover loss or exhibited low performance in predicting forest cover loss in 2015 and 2017, respectively.The significant effect of the pre-drought data of the EVI, SAVI, TCG, and TCW in predicting forest cover loss in 2015 indicates that some parts of Norway's spruce forests are more susceptible to drought than others, possibly due to their differing capacities to retain soil water.Unfortunately, the unavailability of higher spatial and spectral resolution Sentinel 2 (MSI) imagery before 2015 significantly reduced the performance of the used VIs in monitoring and predicting the health status of Norway spruce forests.Nevertheless, the determined sensitivity of VIs to detect the drought-induced response of Norway spruce forests provides valuable insights into the dynamics of forest cover change, with implications for monitoring and conservation efforts in the Kopaonik NP.Accordingly, stakeholders in forestry should seriously consider implementing remote sensing methods within existing forestry management methods for daily use.Data Availability Statement: The data will be made available by the authors upon request.

Conflicts of Interest:
Authors Sr dan Simović and Mirko Dugalić are employed by the Public Enterprise National Park "Kopaonik".The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.The Public Enterprise National Park "Kopaonik" had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Figure 1 .
Figure 1.Study area (location of the Kopaonik NP).

Figure 2 .
Figure 2. Annual forest cover loss after drought in the Kopaonik NP.

Figure 2 . 22 Figure 3 .
Figure 2. Annual forest cover loss after drought in the Kopaonik NP.

Figure 4 .
Figure 4. Time series showing differences in mean values of VIs for forest and non-forest cover (forest cover loss) taken from 2015 (left) and 2017 (right) LCC.

Figure 4 .
Figure 4. Time series showing differences in mean values of VIs for forest and non-forest cover (forest cover loss) taken from 2015 (left) and 2017 (right) LCC.

Figure 5 .
Figure 5. Evaluation of VI sensitivity in predicting forest cover loss in 2015 (gray -no effect, green to red gradient color scale -medium to very large effect).

Figure 5 .
Figure 5. Evaluation of VI sensitivity in predicting forest cover loss in 2015 (gray-no effect, green to red gradient color scale-medium to very large effect).

Table 1 .
The imagery used for the evaluation of drought effects on the forest cover loss.

Table 2 .
The imagery used for mapping forest cover loss.

Table 3 .
The VIs used for the evaluation of drought effects on forest cover loss.

Table 4 .
Dynamics of cumulative forest cover loss at forest glade (fragment) level.

Table 5 .
Accuracy quantitative measures of LCC.

Table 6 .
The error matrix populated by estimated proportions of area.
W i -mapped area proportions (%), A m,i -mapped area.