Changes of Norway Spruce Health in the Białowie ˙za Forest (CE Europe) in 2013–2019 during a Bark Beetle Infestation, Studied with Landsat Imagery

: Among the largest disturbances affecting the health of spruce forests is the large-scale appearance of bark beetles. Knowledge on the spatial distribution of infected-spruce areas is vital for effective and sustainable forest management. Medium-spatial-resolution (20–30 m) satellite images are well-suited for spruce forest disturbance monitoring at a landscape and regional scale following bark beetle outbreaks. The aim of this study was to evaluate the health of a Norway spruce stand after a bark beetle outbreak based on Landsat 8 images and thematic and vector data, supplemented with selected climate variables. This research was conducted for a spruce stand in the Białowie˙za Forest District in 2013, 2015, 2017, and 2019. We hypothesised that the changes in spruce health would signiﬁcantly inﬂuence the NDVI distributions during the studied years. Our research revealed that the weather conditions in the period of May–September were beneﬁcial for beetle development and detrimental for the spruce stand, particularly in 2015, 2018, and 2019. SWIR-NIR-G and NDVI images showed a gradual deterioration in spruce health. The quantitative NDVI distributions varied; the minimum, mean, and median decreased; and the distribution shape of the index values changed over the studied years. An analysis of the spatial NDVI distributions revealed that the threshold NDVI value separating spruce stand areas in good and poor health was ca. 0.6. This study conﬁrmed the applicability of NDVI for monitoring alterations in spruce stands, and indicated that spatial NDVI distributions can provide valuable support in forest monitoring at a landscape scale, since medium-resolution, ready-to-use NDVI images are easily available from the Landsat archives, facilitating the routine assessment of stand health.


Introduction
Forests cover about one third of the Earth's land area and fulfil environmental, productive, and social functions [1,2]. These roles are interdependent and related to each other in time and space. They are also influenced by a variety of external factors. Among the largest disturbances affecting the health of forests in diverse geographical regions are natural hazards caused by the large-scale appearance of insects [3,4]. In Europe and North America, considerable environmental and economic losses are caused by bark beetle (Ips typographus L.) infestations.
Bark beetle is a permanent, integral element of the Norway spruce (Picea abies [L.] Karst.) ecosystem. In normally functioning spruce stands, bark beetle is an element of natural selection, where the beetles feed on trees weakened by abiotic factors. With a low bark beetle population in a region, there are not enough beetles to trigger an infestation of healthy trees [5]. This situation becomes complicated in the case of increasingly frequent drought conditions, induced by low precipitation and an increase in mean air temperatures [6,7]. These factors weaken the spruce stands by disturbing soil hydrological regimes and creating favourable conditions for increasing the number of bark beetle generations during a season. This results in the fast growth of bark beetle populations and a further increase in pressure on the trees. Under such conditions, the predispositions of spruce vegetation do not play a role in bark beetle infestation during the subsequent year, and the beetles start to infest all trees indiscriminately [5].
Bark beetle outbreaks can destroy large areas of spruce stands, thereby changing the ecosystem dynamics of forests and disturbing their structures and functions [8]. Therefore, comprehensive monitoring of the forest health and the identification of areas with extensive damage are necessary [9,10]. Remote sensing techniques that can support time-consuming field surveys are suitable for this purpose [11][12][13][14][15][16][17]. Knowledge on the spatial distribution of infected spruce areas is vital for forest managers to effectively plan appropriate countermeasures and predict the future dynamics of bark beetle infestations [5]. Medium-spatial-resolution (20-30 m) satellite images, such as those from Landsat, are wellsuited for spruce forest disturbance monitoring at a landscape and regional scale following I. typographus outbreaks, and are used in Europe and North America [12,[18][19][20][21][22][23][24][25].
To assess the extent of the spruce stand disturbance, the Normalized Difference Vegetation Index (NDVI) is widely applied. For example, Lausch et al. [26] observed that the best spectral indicators of spruce health in the Bavarian Forest (Germany) are concentrated in the Vis and NIR wavelength regions. Heurich et al. [27] proved that the NDVI is a suitable measure for separating deadwood areas and vital vegetation via objectorientated image analysis of the Bavarian Forest following a spruce bark beetle outbreak. Mišurec et al. [23], for spruce in the Ore Mountains (Czechia), noted that the NDVI detected not only severe forest damage, but also slight differences in the physiological status of trees, which did not exhibit high defoliation levels. Hart and Veblen [28] distinguished grey stands from green stands via NDVI in a beetle-infested spruce-fir zone of the eastern San Juan Mountains (USA).
Recently in Europe, Norway spruce bark beetle infestations have been noted in France, Germany, Austria, and Czechia. In 2012, a large-scale I. typographus outbreak was also detected in the Polish part of the Białowieża Forest [29]. Considering the high value of the Białowieża Forest ecosystem, comprehensive monitoring of forest dynamics supported by remote sensing techniques was conducted, including observations of the spruce stand [30]. Recent studies on this topic involve mostly the single-tree-level detection of spruce health and cover intra-annual bark beetle disturbances [7,16]. However, there is still a need for systematic, long-term observations of this area on a regional scale, which could support the identification of zones requiring more detailed research.
The aim of the study was, therefore, to evaluate the health of the spruce stand after a bark beetle outbreak in the Białowieża Forest District. The analysis was based on Landsat 8 imagery and thematic and vector data, supplemented by selected climate variables. The research was conducted for the Norway spruce stand in 2013, 2015, 2017, and 2019. To assess the ecosystem changes, we employed the visible (G-green, R-red), near-infrared (NIR), and short-wave infrared (SWIR) wavelengths, which can provide information about, respectively, leaf pigments, cell structure, and leaf water content. Furthermore, we verified the applicability of the NDVI for monitoring alterations in the spruce stands. We hypothesised that the changes in the health of Norway spruce would significantly influence the NDVI distributions over the studied years.

Study Area
The Białowieża Forest District, belonging to the Polish State Forests-National Forest Holding, is located in the Polish part of the Białowieża Forest and covers an area of 12,592.71 ha. This Forest District is divided into six forest units: Stoczek, Grudki, Batorówka, Teremiski, Nowe, and Podcerkiew. The Białowieża Forest District is situated in a temperate climatic zone. For the years 1983-2012, during the period May-September, the average daily air temperature was 15.3 • C, the number of days with an average daily temperature of at least 15 • C was 84, and the total precipitation was 353.8 mm (Institute of Meteorology and Water Management-National Research Institute, IMGW-PIB [31,32]).
A part of the Białowieża Forest District, with an area of 3800.3 ha, was investigated. This area was covered by tree stands (3663.2 ha) with 80% mesic forest habitat types. Within the area, tree stands where spruce was the dominant species with at least a 50% share were selected for this study. They covered 620.58 ha (Figure 1).

Data Sources and Data Processing
In this study, we used open access spatial and nonspatial data from diverse sources. The spatial data were processed and the maps were prepared using QGIS 3.6.0, a free and open source GIS program (https://qgis.org/en/site/). Maps were presented in the projection system EPSG:2180.
According to [26], spruces suffering from a bark beetle infestation can only be recognised starting from mid-summer. The autumnal satellite images thus allowed for the assessment of spruce stand disturbances caused by all generations of bark beetles in the current year. From the Landsat imagery archives [33], we selected four Landsat 8 scenes, without cloud cover over the study area, sensed on the following days: 19 October 2013, 2 October 2015, 8 November 2017, and 20 October 2019 (Table 1). Level-2 data from Landsat 8 satellite images were generated through the USGS EROS Center Science Processing Architecture (ESPA) On-demand interface using the LaSRC v. 1.4.1 application [34,35]. The obtained data included SR images for 1-7 bands of the Landsat 8 Operational Land Imager sensor (OLI) and metadata along with NDVI images derived from the Landsat 8 surface reflectance product using the Spectral_Indices v. 2.8.0 application [36,37]. The NDVI images were computed via the ESPA interface from SR bands as (NIR-R)/(NIR+R) or (Band 5-Band 4)/(Band 5+Band 4). The SR products were supplied in the Georeferenced Tagged Image File Format (GeoTIFF). For the study, the SR images from bands 3 (G, 525-600 nm), 4 (R, 630-680 nm), 5 (NIR, 845-885 nm), and 6 (SWIR, 1560-1660 nm) of the Landsat 8 OLI [38,39], and the above-mentioned NDVI images with a 30-m spatial resolution were utilised. Data from the OLI sensor were calibrated to greater than 5% uncertainty in terms of absolute spectral radiance and had an absolute geodetic accuracy better than 65 m with a circular error featuring horizontal accuracy with 90% confidence [40]. Table 1. Description of the data applied in the analysis of the studied spruce stand. To visualise the status of the spruce stand, SWIR-NIR-G (Band 6-Band 5-Band 3) false colour composites were prepared from the RS images ( Figure 2). Representative SWIR, NIR, and G band reflectance charts ( Figure 3) over a 1332 m long transect (the A-B line), as shown in Figure 2, were made for each study year using the QGIS plug-in Profile Tool v. 4.8.1. The spruce health was further assessed using images presenting the spatial distribution of the NDVI values ( Figure 4). For the NDVI images in the range of the studied spruce stand, the count and area of each unique value were measured using the QGIS raster processing algorithm "Raster layer unique values report". Then, descriptive statistics characterising the acquired NDVI distributions were calculated (Figure 5a, Table S1). These distributions were also checked for equality using a λ Kolmogorov-Smirnov test at p < 0.01. The distributions were compared collectively to find all statistically significant differences across the studied data set. Absolute and relative differences in the areas of the NDVI classes of the studied spruce stand were then calculated ( Figure 6, Tables S2 and S3).

Source Data Description Data Type and Application
Weather data for the Białowieża station for the period of 1 May through 30 September, 2013-2019, and the long-term averages for 1983-2012, were calculated based on the data retrieved from the resources of the Institute of Meteorology and Water Management-National Research Institute [31,32]; these data were processed (Table 1, Figure 5b). Research regarding bark beetle ecology (for review see, e.g., [41]) and phenology [42] indicates that the bark beetle population size and subsequent tree mortality are positively and significantly correlated with air temperature increases up to a threshold value, which ultimately impedes bark beetle development. The following climate variables were thus analysed: the average daily air temperature, the number of days with an average daily temperature of at least 15 • C, and precipitation. The thematic and spatial data for the tree stands were obtained from the Forest Data Bank [43]. The data package contained geometric vector layers, nonspatial attribute tables, and metadata (Table 1). The vector layers of the forest unit boundaries (administrative divisions of the Polish State Forests) were obtained directly from the Białowieża Forest District office. The absolute horizontal geodetic accuracy of the G_SUBAREA and G_COMPARTMENT layers was greater than 0.5 m, and for forest unit boundaries, greater than 50 m. Based on the G_SUBAREA layer, we calculated the share of tree species in the study area. Dominant species included (Figure 1b) Norway spruce Picea abies (L.) H.Karst 41.02%, Scots pine Pinus sylvestris L. 39.50%, and black alder Alnus glutinosa Gaertn. 9.35%, with pedunculate oak Quercus robur L. 4.93%, silver birch Betula pendula Roth 2.92%, common hornbeam Carpinus betulus L. 2.15%, and common ash Fraxinus excelsior L. in 0.13% of the tree stand area. Then, the shares of age categories relative to the area of the spruce stand were calculated: 20-40, 9.65%; 40-60, 22.20%; 60-80, 4.58%; 80-100, 20.27%; 100-120, 13.54%; 120-140, 9.12%; >140, 20.63% (Figure 1c). The average age of this spruce stand was 82.1 years. Such a high share of the oldest trees resulted from the policy of the Polish State Forests. Namely, the tree stands of at least 10% of trees of one species aged ≥100 years were excluded from economic (commercial) treatments [44].

Results
Landsat 8 imagery was utilised to assess the health of a spruce stand in the Białowieża Forest District in the years of 2013-2019. First, SWIR-NIR-G composites were analysed ( Figure 2). In the presented composites, SWIR, NIR, and G bands were assigned, respectively, to red, green, and blue colours for display purposes. On the maps, changes in the spruce stand selected for the study could thus be readily observed, revealing the cumulative effects from the beginning of infestation. In the colour composite applied, green indicates the limits of the healthy spruce stand. This colour resulted from the high reflectance in the NIR and, at the same time, 2-3 times lower reflectance in the SWIR. On the other hand, the spruce stand areas in a poor health showed different values of reflectance. We noted a decrease in the NIR accompanied by a considerable increase in the SWIR equalling, or even exceeding, the values observed in the NIR (Figure 3). Consequently, as shown in Figure 2, brown and red-brown colours can be observed in areas of the damaged spruce stand. The Landsat 8 image sensed on 19 October, 2013 presented the studied spruce stand in a relatively good health. This was confirmed by the dominance of green and dark green colours within the boundaries of the spruce stand. In 2015, the area of damaged spruce expanded, which was visible by an increase in brown-coloured spots. The results presented on the 2017 and 2019 maps revealed further deterioration in the health of the spruces, at which point brown colour predominated within the examined stand.
For a more detailed description of the spruce stand health, the Normalized Difference Vegetation Index was applied. The spatial distributions of the NDVI also clearly showed changes in the spruce health ( Figure 4). In 2013, the NDVI ranged from 0.521 to 0.923, with dominant values of 0.7-0.9 and a mean of almost 0.8 (Figure 5a, Tables S1 and S2). The distribution of values was skewed towards high values above ca. 0.75. In 2015, the NDVI decreased, which was indicative of the deterioration in the health of the spruces. The NDVI ranged from 0.433 to 0.854, with a mean close to 0.7 and dominant values of 0.6-0.8. The NDVI distribution, as of 2013, was also shifted towards high values: 75% of the stand area was characterised by an NDVI value above 0.64. On the 2017 map, the NDVI varied from 0.417 to 0.914, with a mean of ca. 0.64. The values from 0.5 to 0.7 predominated, and 75% of the results were below 0.69. In 2019, the NDVI ranged from 0.273 to 0.904, with a mean of 0.62, dominant values of 0.5-0.7, and 75% of the results below 0.66 (Figure 5a). The range of the NDVI values extended throughout the studied dates, mainly due to the decrease in the minimum. Moreover, the mean and median decreased, and the shape of the NDVI distribution changed, which was illustrated by kurtosis and skewness (Table S1). Within 6 years, the area of spruce stand with an NDVI ≥ 0.7 decreased by 444.87 ha, i.e., by over 80% (Figure 6, Tables S2 and S3). A λ Kolmogorov-Smirnov test at p < 0.01 confirmed the presence of statistically significant differences between the NDVI distributions, although in 2017 and 2019, these distributions were very similar (Figure 5a). Figure 7 shows areas with NDVI values in the range of 0.2-0.6 in relation to the age of the studied spruce stand. These maps revealed that low NDVI values were detected initially in the oldest spruce stands (over 140 years), located mainly at the northern part of the study area. In the following years, the zones of low NDVI covered other spruce stands, especially those over 80 years old.

Discussion
Bark beetle infestations are currently a problem in Norway spruce stands in Europe; hence, there is a need to monitor the health of these forests. This task is even more important due to the climate changes observed in recent decades. Weather conditions affect both host trees and bark beetle populations. In the current study, the analysed weather data in the period of May-September for 2013-2019 indicated favourable conditions for bark beetle development in the Białowieża Forest District. At the same time, the relatively low observed precipitation under high air temperatures was considered a stress factor in the studied spruce stand (Figure 5b). This supports previous findings that the probability of infestation is high due to drought even when regional bark beetle population pressure is low, and that variation in water availability is more important than variation in temperature [25].
To assess the tree damage manifested as the dehydration of tissues, shortwave infrared light (SWIR) is particularly useful, since the maximum differences in reflectance attributed to plant water content occur in the 1500-1750 nm spectral region. It is well established that the colour of spruce changes over time, beginning from the tree's colonisation by bark beetles. First, very recently infested trees remain green (green-attack). After the beetles leave the host tree, the spruce crowns become reddish-brown (red-attack). Only 36 months after the bark beetle attack, the spruces start to appear grey in colour (dead trees without needles: grey-attack) [26]. Hart and Veblen [28] observed that the reflectance of NIR and SWIR light, rather than just visible wavelengths, improved estimates of grey-stage spruce beetle outbreaks. Consequently, our SWIR-NIR-G false colour composites showed the effects of bark beetle infestation and allowed us to estimate the damage range ( Figure 2). Changes in the spruce health were reflected by a modification of the spectral responses of the spruces (Figure 3). Healthy spruce stands revealed high reflectance in the NIR and, at the same time, lower reflectance in the SWIR due to the absorption of the latter by the water contained in the plants. Spruce stands in a worse state of health, i.e., those with a loss of photosynthetic organs, the destruction of their mesophyll, and dehydration of their tissues, showed a decrease in the NIR accompanied by a considerable increase in the SWIR. Our results were in line with those of Rautiainen et al. [45] concerning the spectral properties of coniferous forests. Moreover, similar reflectance characteristics in the range 450-2480 nm were obtained by [26] for spruce stands located in the Bavarian Forest National Park.
In the current study, the spatial distributions of NDVI also clearly showed the changes in spruce health (Figure 4). In 2013, the NDVI revealed a mean close to 0.8 (Figure 5a, Tables S1 and S2), with the distribution skewed towards high values above ca. 0.75. This confirmed our earlier observations (Figure 2) that the health of the stand at this time was good. In 2015, the NDVI decreased, which indicated deterioration of the spruce. The 2017 and 2019 maps revealed a further decline in the health of the studied spruce stands (Figure 4). The NDVI distributions differed significantly over the studied years: the NDVI range extended throughout the studied dates, the mean and median decreased, and the shape of the NDVI distribution changed (Figure 5a, Table S1).
Previous studies also indicated intra-and inter-annual NDVI variability following diverse external conditions and canopy phenological stages. Matiu et al. [46]  The above literature review indicates that various Norway spruce disturbances resulted in a decrease of NDVI values in a range of 0.20-0.25. Our research showed similar NDVI values and trends for the quantitative NDVI distributions over the studied years. Additionally, the observed alterations in spatial NDVI distributions suggested that the threshold NDVI value separating the studied spruce stands in good and poor health was ca. 0.6. This value was considered when creating maps of areas with a low NDVI in relation to the spruce stand age (Figure 7). These maps revealed that the oldest spruce stands-over 140 years old and located mainly in the northern part of the study area were the most susceptible to damage. The first 0.2-0.6 NDVI spots were visible in that area during 2013. In the following years, the zones of a low NDVI covered other spruce stands, especially those over 80 years old. Generally, stands older than 80 years were in worse health than the younger ones. This supports the earlier findings of Ksepko et al. [49] that the bark beetle outbreaks in the studied area started in tree stands that were more than 100 years old and which were excluded from protective management and sanitary cuttings. Other studies also revealed that spruces >70-100 years old are at the highest risk of infestation [7,41,50].
Based on Landsat 8 imagery, we observed ongoing deterioration of spruce health in 2013, 2015, 2017, and 2019 and an expansion in the area of damaged spruce stands. The results for 2015 agreed well with the map of the field inspection carried out in September 2015 [49], which was required for creating the documentation to elaborate an annex to the 2012-2021 Forest Management Plan for the Białowieża Forest District. The field survey included 424 forest compartments with stands containing spruce with a total area of 2055.17 ha, which accounted for 29.8% of such stands in the Białowieża Forest District outside the nature reserves.
One of the methods for assessing the dynamics of the bark beetle population is to track the volume of infested deadwood through spruce stand screening or recording the amount of sanitation cuts, which is part of the statutory activities of the forest protection services and forest administration. From the point of view of forest management, continuous ecosystem monitoring is very important, as it serves the current needs of forest protection and facilitates modelling the developments in the number of bark beetles and other threats to spruce stands. Such activities can, and should, be supported by remote sensing techniques [51]. The spatial NDVI distributions presented in this paper could serve as valuable support for monitoring at a landscape scale in the Białowieża Forest and in other forest ecosystems. The recorded distributions are suitable for the routine assessment of stand health since they are based on open access products from the Landsat archive. The ready-to-use NDVI images derived from Level 2 surface reflectance data are readily available upon request. An additional advantage of Landsat imagery is its large scene size, covering ca. 185 × 185 km. Therefore, for example, to monitor the Białowieża Forest, only one scene was sufficient; for the Polish part-scene 187/23 (WRS2 path/row), and for the whole Białowieża Forest within Poland and Belarus-scene 186/23. This allows one to process a limited number of scenes for the monitored forest ecosystem, and thereby, to obtain homogenous cartographic results.

Conclusions
A combination of Landsat 8 medium-resolution satellite imagery, thematic and vector data, supplemented by weather data proved to be an effective tool to objectively assess the current state of the studied spruce stand within the Białowieża Forest District under a bark-beetle infestation.
The analysis of temperature and precipitation for the period of May-September in the years 2013-2019 revealed that the weather conditions were beneficial for bark beetle development and detrimental for the spruce stand, particularly in 2015, 2018, and 2019, when extremely low precipitation accompanied by high mean air temperatures was recorded.
Since the visible (G-green, R-red), near-infrared (NIR), and short-wave infrared (SWIR) wavelengths provide information about, respectively, leaf pigments, cell structure, and leaf water content, SWIR-NIR-G false colour composites and NDVI images clearly showed the effects of bark beetle infestation. We observed ongoing deterioration in spruce health during 2013, 2015, 2017, and 2019. The stands older than 80 years were in a worse health than the younger ones. The NDVI maps showed these changes more precisely, which confirmed the applicability of the index for monitoring alterations in spruce stands. We also found statistically significant differences between the distributions of the NDVI values. The minimum, mean, and median decreased, and the distribution shape of the index values changed, in the years of 2013, 2015, 2017, and 2019. Analysis of the spatial NDVI distributions suggested that the threshold NDVI value separating the spruce stands in good and poor health was ca. 0.6. Our study indicated that the spatial NDVI distributions can provide valuable support for forest monitoring at a landscape scale, especially due to the fact that medium-resolution ready-to-use NDVI images are easily available from the Landsat archives, thereby facilitating the routine assessment of stand health.
Supplementary Materials: The following are available online at https://www.mdpi.com/1999-490 7/12/1/34/s1, Table S1: Statistical parameters of the quantitative NDVI distributions of the studied spruce stand, Table S2: Absolute differences in the area of the NDVI classes of the studied spruce stand, Table S3: Relative differences in the area of the NDVI classes of the studied spruce stand.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.