Continuous Detection of Small-Scale Changes in Scots Pine Dominated Stands Using Dense Sentinel-2 Time Series

: Climate change and severe extreme events, i.e., changes in precipitation and higher drought frequency, have a large impact on forests. In Poland, particularly Norway spruce and Scots pine forest stands are exposed to disturbances and have, thus experienced changes in recent years. Considering that Scots pine stands cover approximately 58% of forests in Poland, mapping these areas with an early and timely detection of forest cover changes is important, e.g., for forest management decisions. A cost-e ﬃ cient way of monitoring forest changes is the use of remote sensing data from the Sentinel-2 satellites. They monitor the Earth’s surface with a high temporal (2–3 days), spatial (10–20 m), and spectral resolution, and thus, enable e ﬀ ective monitoring of vegetation. In this study, we used the dense time series of Sentinel-2 data from the years 2015–2019, (49 images in total), to detect changes in coniferous forest stands dominated by Scots pine. The simple approach was developed to analyze the spectral trajectories of all pixels, which were previously assigned to the probable forest change mask between 2015 and 2019. The spectral trajectories were calculated using the selected Sentinel-2 bands (visible red, red-edge 1–3, near-infrared 1, and short-wave infrared 1–2) and selected vegetation indices (Normalized Di ﬀ erence Moisture Index, Tasseled Cap Wetness, Moisture Stress Index, and Normalized Burn Ratio). Based on these, we calculated the breakpoints to determine when the forest change occurred. Then, a map of forest changes was created, based on the breakpoint dates. An accuracy assessment was performed for each detected date class using 861 points for 46 classes (45 dates and one class representing no changes detected). The results of our study showed that the short-wave infrared 1 band was the most useful for discriminating Scots pine forest stand changes, with the best overall accuracy of 75%. The evaluated vegetation indices underperformed single bands in detecting forest change dates. The presented approach is straightforward and might be useful in operational forest monitoring.


Introduction
Currently, forest ecosystems are exposed to many disturbances related to climate change and extreme weather events. Forests are especially sensitive to climate change, as they cannot rapidly adapt to environmental changes [1]. Changes in the precipitation and a higher drought frequency have an impact on the growth and condition of trees, the occurrence of insects or fungal diseases, and, finally, the increase of mortality risk [1,2]. In temperate zones, forests have been exposed to severe droughts in recent years, which affect forests directly and indirectly [3]. There is high abundance approaches, the integration of MODIS and Landsat data for forest cover change monitoring was used [41,42]. On the other hand, the use of very-high-resolution (VHR) imagery provides a detailed and precise detection of smaller changes [43]. When VHR data are used for forest change mapping, the object-based approach is often employed [44,45]. However, the operational use of VHR imagery may be hampered by the high cost of its acquisition [46]. The methods for mapping forest cover changes include algorithms designed for analyzing the temporal trends of forest disturbances, such as LandTrendr [33], bfastSpatial [20], the Continuous Change Detection and Classification (CCDC) algorithm [47], and the TIMESAT program [48]. These algorithms are specifically designed to analyze long time series of Landsat, NOAA Advanced Very-High Resolution Radiometer (AVHRR), and MODIS imagery.
Data from the Sentinel-2 mission are characterized by a high temporal (up to 5 days), spatial (10, 20 m) and spectral resolution (13 bands), allowing for the detection of small forest cover changes with higher accuracy than that achieved using medium-or coarse-resolution data [24]. Furthermore, the Sentinel-2 imagery is available for free, which, in combination with the variety of open-source software, creates opportunities relating to the successful application of this approach in operational forest health monitoring. So far, the Sentinel-2 time series have been used for the assessment of forest cuttings [49] and treefall gaps [50], as well as the monitoring of cork oak decline [51]. Sentinel-2 was also used in the mapping of bark beetle green attack stages [30], Scots pine stands defoliation [31], and red-edge bands were assessed in forest decline [26]. In general, the Sentinel-2 data were found to be very successful in forest disturbance mapping, and their spatial resolution enables the mapping of changes at a detailed level.
Accordingly, the aim of this study is to develop a simple approach for the continuous monitoring of forest changes using the Sentinel-2 dense time series at a local scale. We focused on changes in Scots pine forest stands, which, since 2015, show symptoms of dieback in Poland. We analyzed all the pixels within the study area, which were previously assigned to the probable change area, based on classifications. We used open-source software and the free Sentinel-2 data to provide an approach, which could later be applied in operational forest management. In particular, this study had the following objectives: To evaluate the dense Sentinel-2 time series for the continuous monitoring of Scots pine stands changes; To compare the usefulness of the Sentinel-2 bands and vegetation indices in detecting forest disturbances; and To assess the temporal accuracy of the detected changes.

Study Area
The study area is the Milicz Forest District (Figure 1) located in South-western Poland on the border of Milicz-Głogów Depression, Trzebnica Range, and Southern Wielkopolska Lowland macroregions [52]. It covers an area of 656 km 2 , and the forests cover approximately 42% of it. Scots pine (Pinus sylvestris) stands cover approximately 70% of the forests in the Milicz Forest District [53]. This site was selected due to reported changes related to beetle disturbances.

Satellite Imagery
The Sentinel-2 imagery was downloaded from the Copernicus Open Access Hub repository using the sen2r package in R [54]. All of the available Sentinel-2 data from the years 2015-2019 (tile number: 33UXT), with a cloud cover below 10%, were downloaded and further processed, excluding winter imagery. In total, we used 49 images, acquired between 11 July 2015, and 23 September 2019 ( Figure 2).

Satellite Imagery
The Sentinel-2 imagery was downloaded from the Copernicus Open Access Hub repository using the sen2r package in R [54]. All of the available Sentinel-2 data from the years 2015-2019 (tile number: 33UXT), with a cloud cover below 10%, were downloaded and further processed, excluding winter imagery. In total, we used 49 images, acquired between 11 July 2015, and 23 September 2019 ( Figure 2).

Satellite Imagery
The Sentinel-2 imagery was downloaded from the Copernicus Open Access Hub repository using the sen2r package in R [54]. All of the available Sentinel-2 data from the years 2015-2019 (tile number: 33UXT), with a cloud cover below 10%, were downloaded and further processed, excluding winter imagery. In total, we used 49 images, acquired between 11 July 2015, and 23 September 2019 ( Figure 2).

Methodology
To achieve the objectives of the study, the following steps of analysis were performed: (i) The Sentinel-2 data were downloaded and pre-processed, including sen2cor correction, vegetation indices (VI) calculation, and cloud masking; (ii) "base" and "target" classifications were created, and their accuracy was assessed; based on these classifications, (iii) mask of probable change was generated; (iv) the spectral trajectories of the Sentinel-2 bands and indices were calculated for all the pixels in the area of probable change; (v) the accuracy of the detected changes was assessed using 861 validation points; and (vi) a forest change map was created, based on the highest accuracy change detection ( Figure 3). All of the abovementioned steps were performed in R [55].
Remote Sens. 2020, 12, x FOR PEER REVIEW 5 of 20 To achieve the objectives of the study, the following steps of analysis were performed: (i) The Sentinel-2 data were downloaded and pre-processed, including sen2cor correction, vegetation indices (VI) calculation, and cloud masking; (ii) "base" and "target" classifications were created, and their accuracy was assessed; based on these classifications, (iii) mask of probable change was generated; (iv) the spectral trajectories of the Sentinel-2 bands and indices were calculated for all the pixels in the area of probable change; (v) the accuracy of the detected changes was assessed using 861 validation points; and (vi) a forest change map was created, based on the highest accuracy change detection ( Figure 3). All of the abovementioned steps were performed in R [55].

Data Pre-Processing
We downloaded all of the available Sentinel-2 imagery, with a cloud cover lower than 10%. For the imagery from after 2017, Bottom-of-Atmosphere (BOA) products were downloaded (if available), and for the older imagery, the sen2cor algorithm from the sen2r package was applied to transform the Top-of-Atmosphere (TOA) into BOA products. Then, four indices were calculated using s2_calcindices: NDMI, MSI, NBR, and TCW (Table 1). In the next step, we applied cloud masking to the Sentinel-2 BOA imagery and indices. The cloud mask was derived from the Sentinel Land Cover classification product (clouds, clouds shadows, and snow have been masked) [56]. The resulting products had a resolution of 10 meters.

Classification and Probable Change Mask
Post-classification change is one of the most common methods of change detection, and it is simple and suitable for many applications. However, its accuracy is dependent on the accuracy of the

Data Pre-Processing
We downloaded all of the available Sentinel-2 imagery, with a cloud cover lower than 10%. For the imagery from after 2017, Bottom-of-Atmosphere (BOA) products were downloaded (if available), and for the older imagery, the sen2cor algorithm from the sen2r package was applied to transform the Top-of-Atmosphere (TOA) into BOA products. Then, four indices were calculated using s2_calcindices: NDMI, MSI, NBR, and TCW (Table 1). In the next step, we applied cloud masking to the Sentinel-2 BOA imagery and indices. The cloud mask was derived from the Sentinel Land Cover classification product (clouds, clouds shadows, and snow have been masked) [56]. The resulting products had a resolution of 10 m.

Classification and Probable Change Mask
Post-classification change is one of the most common methods of change detection, and it is simple and suitable for many applications. However, its accuracy is dependent on the accuracy of the input maps [22,57,58]. In some studies on forest disturbances, a forest mask was created for the beginning of the time series [59]. In our study, besides using this "base" forest map, we also created a "target" classification map for the end of the studied period. Both of the maps were classified into two classes: Coniferous forests (pine stands) and all other land cover types (broad-leaved forests, low vegetation, bare land, built-up areas, and water bodies). Reference samples for classification and validation were visually delineated for areas where no changes occurred during the studied period, and they include 112 polygons for conifer forest stands and 120 polygons for non-conifer forests (broad-leaved, water, agricultural, areas, and built-up). These reference polygons were randomly split into training and validation polygons in a 70:30 ratio. The base classification was calculated for the first available Sentinel-2 image for the study area, i.e., 11 July 2015. The target classification was calculated for the image from 31 August 2019. Based on the differences between the maps for 2015 and 2019, we produced the change map.
The classification was performed using the Random Forest (RF) algorithm, which is currently one of the most broadly-used classifiers in remote sensing land cover studies, mainly because of its ease of use and high accuracy [60]. Random Forest is an ensemble classifier that consists of a combination of decision trees [61]. The hyperparameters of the RF classifier consist of the number of trees that were set to 500 and a number of predictors at each split, which was tested from values of 1,2,3,5, and 10. The classification was performed using the superClass function from the RStoolbox package [62], and the rf classification algorithm from the randomForest package [63]. After the classification of single images, the accuracy assessment was performed using a validation dataset, and the achieved accuracies for both of the classifications exceeded 98% of the overall accuracy (Table 2).

Pixel Trajectory Assessment
For each pixel assigned to the probable change mask, we assessed the spectral trajectories, i.e., the values from the bands and indices for the whole studied period of time. In the case of spectral bands and TCW, the "raw' values were used, while for NBR, NDMI and MSI indices the values were rescaled into the range from 0 to 1, which, based on tests, provided the best results regarding changes detection. To determine the date of a change (i.e., disturbance), we used the breakpoints function from the strucchange package [64]. This function computes the optimal breakpoints by minimizing the residual sum of squares in the regression model [65]. The parameters of the breakpoints function were set based on test runs, and the minimum length of a segment (h) was equal to two. Additionally, the confidence intervals for the model parameters in the breakpoints were calculated at a confidence level of 0.95, and we rejected breakpoints with intervals higher than 2. As the breakpoints function allow to detect multiple breaks, the first detected one was considered as a change in the forest. We tested breakpoints for reflectance values extracted from the visible red band (VIS R), red-edge bands (RE1, RE2, RE3), the near-infrared band (NIR), and two short-wave infrared (SWIR1, SWIR2) bands. Additionally, we tested four vegetation indices: NDMI, MSI, NBR, and TCW.

Accuracy Assessment and Forest Change Map
For the accuracy assessment of forest change maps we used the following measures of accuracy: An error matrix, an Overall Accuracy (OA), and the producer's/user's accuracies. Additionally, we calculated the macro-averaged F1 score, which is a harmonized mean of the producer's and user's accuracies [66]. To obtain reference samples, we used stratified random sampling, and 861 points were selected for accuracy assessment. Then, they were visually inspected based on the color and texture of the Sentinel-2 time series, and the date of the change (or no change class) was assigned to each of them. Finally, based on the highest accuracy change detection, forest change map was created, where each class represented the date of change (or no change).

Map of Coniferous Forests and Probable Changes
The forest mask for 2015 covered approximately 158.2 km 2 , and that for 2019 covered approximately 158.4 km 2 . The potential change mask, where a change from coniferous forest to any non-coniferous class occurred, covered 4.95 km 2 , which equals 49,511 pixels ( Figure 4). The increase of forest mask by 0.2 km 2 indicates that the forest cover gain is also observed within the study area, and it is higher than forest cover loss. However, the main interest of this study was to analyze the areas where the loss of coniferous forests occurred.

Breakpoints Detection and Accuracy Assessment
The changes in the coniferous forest stands caused an increase in the values in SWIR, RE1, VIS R bands, and MSI, and a decrease in NBR, TCW, and NDMI ( Figures 5 and 6).
The examples of pixel trajectories with both correctly and incorrectly detected breaks are shown in Figures 5 and 6. In some cases, the breaks were detected after the actual break ( Figure 5). There are clearly visible patterns of phenological trajectories of Scots pine stands ( Figure 6); however, in most of the cases, they were not detected as breaks. The variations in reflectance and indices values are significantly larger after the change occurs.

Breakpoints Detection and Accuracy Assessment
The changes in the coniferous forest stands caused an increase in the values in SWIR, RE1, VIS R bands, and MSI, and a decrease in NBR, TCW, and NDMI ( Figures 5 and 6).    The examples of pixel trajectories with both correctly and incorrectly detected breaks are shown in Figures 5 and 6. In some cases, the breaks were detected after the actual break ( Figure 5). There are clearly visible patterns of phenological trajectories of Scots pine stands ( Figure 6); however, in most of the cases, they were not detected as breaks. The variations in reflectance and indices values are significantly larger after the change occurs.
The highest overall accuracy of forest change detection was obtained using the SWIR1 band (75.1% of OA; Table 3). SWIR1 was followed by RE1 (69.5% of OA) and SWIR2 (65.0% of OA). In The highest overall accuracy of forest change detection was obtained using the SWIR1 band (75.1% of OA; Table 3). SWIR1 was followed by RE1 (69.5% of OA) and SWIR2 (65.0% of OA). In general, single Sentinel-2 bands performed better than the vegetation indices. Among the indices, the best performance in detecting the changes was provided by the TCW (59% OA). Moderate accuracies above 60% of OA were also provided by the VIS R bands. The RE2, RE3, and NIR1 bands obtained very low accuracies (OA of 10% or below), and they were, thus not furthered analyzed. In the case of the best detection (SWIR1; Figure 7), the incorrectly detected pixels constituted 24.9% of all reference pixels. Among these misclassified pixels, 61.7% were detected after the actual break, 32.5% were not detected at all, and only 5.8% of all misclassified pixels were detected before the actual break. Furthermore, among the pixels detected after an actual change, 29% were detected within one date, and 50% within two dates after the actual change. Therefore, if also the neighboring detected dates were considered as positive detection, the overall accuracy increased to 80%. The assessment of temporal accuracy of the obtained forest change map, i.e., when the changes were detected with higher or lower accuracy, also showed some regularities (Figure 8). The lowest (zero or near-zero) macro-averaged F1 accuracies characterized early spring but, in all of the cases, not the first available spring date. Stable F1 accuracies generally characterized the summer months.
Remote Sens. 2020, 12, x FOR PEER REVIEW 11 of 20 general, single Sentinel-2 bands performed better than the vegetation indices. Among the indices, the best performance in detecting the changes was provided by the TCW (59% OA). Moderate accuracies above 60% of OA were also provided by the VIS R bands. The RE2, RE3, and NIR1 bands obtained very low accuracies (OA of 10% or below), and they were, thus not furthered analyzed. In the case of the best detection (SWIR1; Figure 7), the incorrectly detected pixels constituted 24.9% of all reference pixels. Among these misclassified pixels, 61.7% were detected after the actual break, 32.5% were not detected at all, and only 5.8% of all misclassified pixels were detected before the actual break. Furthermore, among the pixels detected after an actual change, 29% were detected within one date, and 50% within two dates after the actual change. Therefore, if also the neighboring detected dates were considered as positive detection, the overall accuracy increased to 80%. The assessment of temporal accuracy of the obtained forest change map, i.e., when the changes were detected with higher or lower accuracy, also showed some regularities (Figure 8). The lowest (zero or near-zero) macro-averaged F1 accuracies characterized early spring but, in all of the cases, not the first available spring date. Stable F1 accuracies generally characterized the summer months.

Forest Changes Map
Based on the best detection in terms of accuracy (SWIR 1), we produced changes map in coniferous forest stands for the Milicz Forest District. Out of the 49,511 pixels, the breakpoints were detected for 42,971, which equals an area of the changes of 429.7 hectares (Figure 9). The most numerous changes were detected in the spring months, with the largest change area detected on 2017/04/01, which was equal to nearly 59 hectares, followed by that on 2016/03/17 (  The obtained changes map allowed for the assessment of the spatial patterns and sizes of the coniferous forest changes ( Figure 10). The area of the particular changes in coniferous forests varied from 100 m 2 (i.e., a single Sentinel-2 pixel) to 29,100 m 2 . The largest forest changes were characterized by spatial homogeneity and regular shapes. In many cases, gradual changes were visible on consecutive dates.

Forest Changes Map
Based on the best detection in terms of accuracy (SWIR 1), we produced changes map in coniferous forest stands for the Milicz Forest District. Out of the 49,511 pixels, the breakpoints were detected for 42,971, which equals an area of the changes of 429.7 hectares (Figure 9). The most numerous changes were detected in the spring months, with the largest change area detected on 2017/04/01, which was equal to nearly 59 hectares, followed by that on 2016/03/17 (

Forest Changes Map
Based on the best detection in terms of accuracy (SWIR 1), we produced changes map in coniferous forest stands for the Milicz Forest District. Out of the 49,511 pixels, the breakpoints were detected for 42,971, which equals an area of the changes of 429.7 hectares (Figure 9). The most numerous changes were detected in the spring months, with the largest change area detected on 2017/04/01, which was equal to nearly 59 hectares, followed by that on 2016/03/17 (  The obtained changes map allowed for the assessment of the spatial patterns and sizes of the coniferous forest changes ( Figure 10). The area of the particular changes in coniferous forests varied from 100 m 2 (i.e., a single Sentinel-2 pixel) to 29,100 m 2 . The largest forest changes were characterized by spatial homogeneity and regular shapes. In many cases, gradual changes were visible on consecutive dates. The obtained changes map allowed for the assessment of the spatial patterns and sizes of the coniferous forest changes ( Figure 10). The area of the particular changes in coniferous forests varied from 100 m 2 (i.e., a single Sentinel-2 pixel) to 29,100 m 2 . The largest forest changes were characterized by spatial homogeneity and regular shapes. In many cases, gradual changes were visible on consecutive dates.

Discussion
In this study, we developed a simple approach for detecting changes in Scots pine forest stands in the temperate zone using the Sentinel-2 dense time series. We determined pixel trajectories for the selected Sentinel-2 bands and vegetation indices for all the pixels assigned to the probable change mask between the years 2015 and 2019, with a total of 49,511 pixels (495 hectares). Our approach allowed for the detection of forest disturbances, with a satisfying accuracy of 75% of OA obtained for 46 classes (45 dates and one class representing no changes). We analyzed the sensitivity of different indices and Sentinel-2 bands in the detection of breaks. Interestingly, the indices commonly used in analyzing forest disturbances, such as NDMI or NBR, provided lower accuracies in change detection than the SWIR1 and RE1 Sentinel-2 bands.
Dense time series of Sentinel-2 imagery allows for continuous monitoring of forest ecosystems. We were able to detect Scots pine forest changes with the accuracy of 75%, and it increased to 80% if the neighboring detected dates were considered. Our results are in line with other studies on the Sentinel-2 time series in forest disturbances monitoring, which number, however, is still limited. In forest disturbance studies using the Sentinel-2 data, Puletti and Bascietto [49] used the NDVI series to separate cut from uncut forests, with high accuracy. The potential of the Sentinel-2 series was confirmed in the study of Navarro et al. [51] in the monitoring of cork oak decline. Puletti et al. [67] used Sentinel-2 time series to analyze the vegetation condition following extreme drought events and showed a reduction of trees photosynthetic activity.
In our study, particularly the SWIR1 band performed effectively in detecting changes in Scots pine stands. This part of the spectrum is reported as a good indicator of vegetation water content [68] and is less impacted by soils and atmospheric noise [22]. Infested or dead trees have a significantly

Discussion
In this study, we developed a simple approach for detecting changes in Scots pine forest stands in the temperate zone using the Sentinel-2 dense time series. We determined pixel trajectories for the selected Sentinel-2 bands and vegetation indices for all the pixels assigned to the probable change mask between the years 2015 and 2019, with a total of 49,511 pixels (495 hectares). Our approach allowed for the detection of forest disturbances, with a satisfying accuracy of 75% of OA obtained for 46 classes (45 dates and one class representing no changes). We analyzed the sensitivity of different indices and Sentinel-2 bands in the detection of breaks. Interestingly, the indices commonly used in analyzing forest disturbances, such as NDMI or NBR, provided lower accuracies in change detection than the SWIR1 and RE1 Sentinel-2 bands.
Dense time series of Sentinel-2 imagery allows for continuous monitoring of forest ecosystems. We were able to detect Scots pine forest changes with the accuracy of 75%, and it increased to 80% if the neighboring detected dates were considered. Our results are in line with other studies on the Sentinel-2 time series in forest disturbances monitoring, which number, however, is still limited. In forest disturbance studies using the Sentinel-2 data, Puletti and Bascietto [49] used the NDVI series to separate cut from uncut forests, with high accuracy. The potential of the Sentinel-2 series was confirmed in the study of Navarro et al. [51] in the monitoring of cork oak decline. Puletti et al. [67] used Sentinel-2 time series to analyze the vegetation condition following extreme drought events and showed a reduction of trees photosynthetic activity.
In our study, particularly the SWIR1 band performed effectively in detecting changes in Scots pine stands. This part of the spectrum is reported as a good indicator of vegetation water content [68] and is less impacted by soils and atmospheric noise [22]. Infested or dead trees have a significantly higher reflectance in SWIR than healthy ones at both the leaf and canopy level [30,51]. SWIR was also reported to be insensitive to internal shadows [39]. Schroeder et al. [69] reported that the SWIR part of the spectrum is the most effective in discriminating fires and clear-cuts. Sentinel-2, in comparison to other sensors, incorporates three red-edge spectral bands, which are highly sensitive to chlorophyll and nitrogen content and unaffected by structural properties [27,51]. In our study, the RE1 band provided high accuracy of 69% and, in contrast to RE2 and RE3, enabled the mapping of conifer forest changes. This significant difference in the performance of red-edge bands is due to very significant changes in this part of the spectrum as the wavelength increases ( Figure 11). The red-edge part of the spectrum was also found to accurately retrieve green LAI and canopy chlorophyll content in the study of Delegido et al. [70]. The sensitivity of the red-edge region to the temporal dimension of the forest condition was reported by Zarco-Tejada et al. [26]. The study of Korhonen et al. [71] indicated a higher correlation between single bands, especially RE1, and the canopy cover and leaf area index than indices. Gitelson et al. [72] also reported the high sensitivity of reflectance, near 700 nm, to the chlorophyll concentration.
In other studies, where satellite imagery was used in forest disturbance mapping-mainly from Landsat missions-the indices were found to be an effective method. For example, the evaluation of bands and indices from the Landsat time series, based on a disturbance signal-to-noise metric, showed SWIR reflectance, and SWIR-based indices had the potential for disturbance detection [73]. SWIR-based indices-NDMI and NBR-were also found to be the most sensitive to forest changes among other vegetation indices, such as NDVI or EVI, in selective logging detection [39]. Hughes et al. [74] confirmed the potential of NDMI and NBR indices for forest change detection using the Landsat time series. In separating types of disturbance indices, utilizing SWIR bands was found to be of high importance in the study of Huo et al. [75]. Also, Abdullah et al. [30] confirmed that the Sentinel-2 based, red-edge-and SWIR-dependent indices had a high potential for detecting vegetation changes. In the study of Perbet et al. [76], Change Vector Analysis, based on the NIR, SWIR1, and VIS R bands, was performed in detecting deforestation due to their ability to discriminate forests from bare soil. Many forest disturbances in temperate zones of Europe are small-scale, making the Sentinel-2 data a great tool for their detection. However, there are several sources of errors and uncertainties that should be considered in further studies. Firstly, there may be errors derived from the probable change mask due to errors in the base and target classifications. This is reported in many studies on change detection. Thus, these steps always have to be taken very carefully. Secondly, there are noises (h) NDMI. The example pixel (red dot on true-color composition) was correctly detected using both indices but falsely detected using single bands as a change from 27 March, which was the consecutive date in the Sentinel-2 time series.
Different forest types and parameters may influence the performance of indices and bands in detecting changes. Additionally, different types of disturbances may have an impact on the performance of particular indices [39]. Among the four indices used in this study, the TCW performed better than the other indices in detecting changes. All of the examined indices, by definition, take advantage of the contrast between the SWIR and NIR parts of the spectrum. On the other hand, TCW, utilizes information also from other bands, including the SWIR2 and VIS R bands, which both performed moderately well in this study. Additionally, TCW has been reported to be insensitive to a topographically induced illumination angle [77].
According to our results, 429.7 hectares of coniferous stands were disturbed/deforested since 2015 in the studied area. Particularly in recent years, there is an increasing risk of drought in Poland as one of the effects of climate change. In Poland, a severe drought occurred in 2015 in almost all of the forests, particularly in Central and South-western Poland, where our study site is located [78]. A lack of sufficient amount of water available for the trees caused, among other things, a decrease in the health status of this area [78]. Currently, the climate change in temperate regions is related to higher temperatures, which increase tree water stress, cause physiological effects, and increase the tree vulnerability [3]. In central Europe, particularly conifer species, such as spruce and pine, are exposed to disturbances. Also, in other regions, e.g., in the Mediterranean basin, forests are exposed to drought effects [67]. However, there are different species-and site-specific aspects of growth responses to drought [79]. To evaluate the exact area of forests affected directly by drought, further studies in determining the types of disturbances should be carried out. In particular, planned clear cuts should be separated from disturbances caused by biotic and abiotic factors. Additionally, an explanation of the specific factors influencing the sensitivity to a particular change should be provided. As seen in Figure 10, the changes in forests are sometimes very rapid and occur within the months or even days. Therefore, the use of Sentinel-2, even taking into account the irregularity of acquisition, allows detecting forest disturbances with high temporal accuracy.
Many forest disturbances in temperate zones of Europe are small-scale, making the Sentinel-2 data a great tool for their detection. However, there are several sources of errors and uncertainties that should be considered in further studies. Firstly, there may be errors derived from the probable change mask due to errors in the base and target classifications. This is reported in many studies on change detection. Thus, these steps always have to be taken very carefully. Secondly, there are noises in the time series due to, e.g., undetected clouds or cloud shadows, as well as forest shadows. Forest shadows, which are very problematic in early spring and late autumn imagery, can result in errors, for example, changes detected at later dates. As seen in Figure 11, parts of the non-forested area covered by shadows are characterized by much lower values in single bands, although they are not so evident in the indices. The problem of object shadows is typically known in the processing of very-high-resolution imagery. Thus, the Sentinel-2 data with a 10-and 20-m resolution are now also sensitive to this issue. However, the indices are less sensitive to forest shadows. This problem might be solved in future studies by, for example, applying digital terrain models for shadow detection and then de-shadowing the problematic areas. Furthermore, frequent cloud cover results in the acquisition of irregular time series and sometimes missing data for the key periods of time. The Sentinel-2 time series used in this study are irregularly distributed over the studied years and growing seasons, which results from both the time of launching of Sentinel-2A and -2B missions and the atmospheric conditions, i.e., cloud cover. Thus, the date of the detection may not be an exact timing of the disturbance. This is partially a reason for the most numerous changes detected in the first spring dates, which in fact may also include of the changes from the winter period. Still, this information may be very useful for e.g., forest management, providing precise information where the change in forest stands occurred.
Another factor affecting errors in change detection is the misregistration of the Sentinel-2 imagery, which is visible as the geometric shifts for different acquisitions. In this study, there were several cases of false change detections on the border of forest/non-forest areas due to pixel shifts. Both the internal forest shadows and misregistration of the Sentinel-2 imagery may be reasons for the better performance of bands with a 20-m, rather than a 10-m resolution. The problem of misregistration was reported in other studies using the multi-temporal Sentinel-2 imagery [80]. Yan et al. [81] reported that, after applying a new processing baseline by ESA, the mean misregistration was about 0.4 of 10 m pixels.
The accuracy of the obtained forest change map results, in large part, from the abovementioned factors. Additionally, our approach has some limitations, i.e., it was tested on Scots pine stands, which have recently experienced a lot of disturbances. However, when analyzing broad-leaved forests, seasonal effects caused by phenology influence the result to a larger extent. Still, conifer species also exhibit some seasonal reflectance variations [82], which is seen in Figures 5 and 6. Thus, they, in some cases, may also have an influence on the change detection algorithms.
In detecting forest changes, we used methods from the strucchange R package, designed to detect structural breaks [64]. As a change detection algorithm, it was used in the study of Hislop et al. on the Landsat time series [83], and it performed less accurately than LandTrendR. Also, Oeser et al. [59] used strucchange method to detect disturbance agents with high accuracy. In our study, it performed well and was easily applied to the Sentinel-2 time series. Therefore, it should be tested in other vegetation change studies using Sentinel-2.
While the definitions of forest disturbances differ, most disturbance events are related to the rapid removal or reduction of forest canopy [84]. Still, there are also subtle, non-stand replacements or gradual forest changes, such as thinning or stress damages. In comparison to abrupt changes in forest canopies and forest cover, these are more difficult to discriminate in terms of both their precise detection and the exact timing of their changes. Additionally, it might not be possible to detect changes covering, e.g., single trees, with the Sentinel-2 resolution. In our study, we focused mainly on stand replacement disturbances (i.e., forest stem removal) during this four-year time period. In the future, a longer time series of the Sentinel-2 imagery will be an invaluable source of information for forest and vegetation change analysis, so the automated methods of processing and detecting changes should be developed and used.

Conclusions
This study evaluated the use of the Sentinel-2 time series, consisting of 49 images from 2015 to 2019, in the detection of small-scale changes in coniferous forest stands. The comparison of selected Sentinel-2 bands and vegetation indices, based on an analysis of the breakpoints in pixels trajectories, showed varied usefulness in detecting forest disturbances. Surprisingly, the single Sentinel-2 SWIR1 and RE1 bands performed better than the analyzed vegetation indices, among which the Tasseled Cap Wetness provided the highest accuracy. However, the indices showed higher insensitivity to forest shadows. Given the high accuracy obtained with this simple approach, it might be useful in operational forest management to monitor forest changes.