Using Intra-Annual Landsat Time Series for Attributing Forest Disturbance Agents in Central Europe

The attribution of forest disturbances to disturbance agents is a critical challenge for remote sensing-based forest monitoring, promising important insights into drivers and impacts of forest disturbances. Previous studies have used spectral-temporal metrics derived from annual Landsat time series to identify disturbance agents. Here, we extend this approach to new predictors derived from intra-annual time series and test it at three sites in Central Europe, including managed and protected forests. The two newly tested predictors are: (1) intra-annual timing of disturbance events and (2) temporal proximity to windstorms based on prior knowledge. We estimated the intra-annual timing of disturbances using a breakpoint detection algorithm and all available Landsat observations between 1984 and 2016. Using spectral, temporal, and topography-related metrics, we then mapped four disturbance classes: windthrow, cleared windthrow, bark beetles, and other harvest. Disturbance agents were identified with overall accuracies of 76–86%. Temporal proximity to storm events was among the most important predictors, while intra-annual timing itself was less important. Moreover, elevation information was very effective for discriminating disturbance agents. Our results demonstrate the potential of incorporating dense, intra-annual Landsat time series information and prior knowledge of disturbance events for monitoring forest ecosystem change at the disturbance agent level.


Introduction
Forest disturbances are central drivers of ecosystem dynamics, impacting stand structure and composition, carbon flux, biodiversity, and economic value of forests [1][2][3].Ecological and economic impacts may vary strongly between different types of forest disturbance [4][5][6].Thus, information on the occurrence and distribution of disturbance agents is important to assess the carbon balance of forests [5][6][7], impacts on biodiversity [8][9][10], and economic consequences [7,11].Despite the ecological and economic importance of forest disturbances, there is a lack of spatially explicit information on forest disturbances at the agent level, especially covering larger geographic extents [12].This lack of Forests 2017, 8, 251 2 of 24 knowledge hampers the development of process-based ecological models, which might help with predicting changes in disturbance dynamics under future climate change, as well as the development of adequate management strategies [12,13].
Recent advances in satellite-based monitoring of forest ecosystem dynamics have demonstrated the potential for remote sensing applications to map forest disturbances at the agent level.In this regard, the Landsat satellite program constitutes a unique data source, allowing mapping forest dynamics at a stand scale over large geographic extents and long historic periods (>30 years) with relatively high temporal resolutions (annual or higher) [14].Following the opening of the USGS Landsat Global Archive in 2008 [15], new methods have been proposed for exploiting the full potential of the high repetition rate offered by Landsat.That way, several automated change detection algorithms have been developed, which operate by fitting a simplified model to a time series of Landsat observations [16][17][18][19][20]. Approaches based on Landsat time series have been shown to be very effective for detecting and characterizing different processes of forest ecosystem dynamics [21,22].
While the detection of many forest disturbance processes using Landsat data may be considered operational given the diversity of successful applications in forest ecosystems around the world [14,23,24], the discrimination of different disturbance causes (i.e., the attribution to disturbance agents) remains a prevailing challenge for remote sensing research [25].First studies have used the outputs of time series based change detection algorithms to attribute disturbance agents through classification models [26,27].Such approaches allow for the in-depth study of agent-level forest disturbance dynamics at regional to continental scales [28,29].However, current studies mostly mapped forest disturbances in North-America, where clear-cut harvest, fire, and large-scale insect outbreaks dominate as disturbance agents [26,27,30,31].Approaches developed in those forest ecosystems might be difficult to transfer to regions of lower disturbance intensity and smaller disturbance patch sizes, such as the forest ecosystems of Central Europe [32].
In Central Europe, windthrow and bark beetles are the major natural forest disturbance agents [2,33], which have been shown to interact in complex ways [34,35].Both levels of windthrow and bark beetle disturbances have increased in European forests over the last decades [2,36], and are expected to further increase in response to climate change [37].Given the wide range of important ecological and economic impacts associated with both disturbance agents [5, 7,11], the monitoring of windthrow and bark beetle disturbances is critical to assist the sustainable stewardship of forest ecosystems in Europe [12].Thus, remote sensing-based monitoring of disturbance agents in Europe requires methods that can robustly identify and discriminate these two most influential natural disturbance agents.While there is a range of studies on Landsat-based characterization and identification of insect-related forest disturbances [27,[38][39][40], windthrow disturbances have received considerably less research attention [41].Further, most of Central European forests are under active management, where harvest-in particular selective logging (i.e., thinning or shelterwood cutting)-is an important disturbance agent [42].Also, stands disturbed by natural disturbance are commonly salvaged and sanitation felling is routinely applied to dampen the spread of bark beetles [43][44][45][46].These complex interactions between natural and human disturbances make the attribution of disturbance agents challenging.
Spectral-temporal metrics derived from Landsat time series have been successfully used to identify and discriminate different disturbance agents [27,29,40,47,48].Such metrics can be obtained from automatic change detection algorithms and provide key characteristics of spectral-temporal trajectories such as the spectral change magnitude or the duration of spectral change.While spectral-temporal metrics can deliver important insights into the characteristics of forest disturbances [49], recent studies have demonstrated the potential of improving the identification of disturbance agents by incorporating additional sources of information [26,30,31].These additional information streams include context-based metrics, which relate the location of detected disturbances to landscape characteristics influencing the predisposition of forest stands towards different disturbance types Forests 2017, 8, 251 3 of 24 (e.g., topography or forest type [26,31]), as well as spatial metrics describing the spatial configuration of disturbance patches [26,30,50].
In this study, we investigate the potential of detailed temporal information derived from intra-annual Landsat time series to improve the mapping of disturbance agents.Previous studies on disturbance agent attribution have used time series of annual best observation images [26,29,30].At the same time, great advances have been made regarding change detection methods including intra-annual Landsat data [20,51,52].For instance, Zhu et al. [20] have developed a change detection algorithm working on all available Landsat observations that detects forest disturbances by comparing new observations against values predicted by a model fitted to stable forest pixels.Hilker et al. [51] used a fusion model incorporating Landsat and MODIS data for detecting forest disturbances with high temporal and spatial resolution.As different types of forest disturbances often exhibit a distinct seasonality, such as windthrows [53,54], insect outbreaks [55,56] or fires [57,58], information on the timing of disturbance events within a year may potentially help to infer the underlying disturbance agent [13].
Moreover, combining Landsat time series data with additional information on natural disturbance events and/or their underlying causes may improve the attribution of forest disturbances.Previous studies combined climate data [59] or independent disturbance observations from aerial surveys [60][61][62] to predict insect outbreak dynamics.Here, we test whether including external information on the occurrence of storm events will help improve the identification of storm-related disturbances.
The aim of this study was to test Landsat time series-based predictors, context-based topography-related metrics, and prior information on wind disturbances for mapping forest disturbance agents in Central Europe.We tested our approach at three study sites featuring varying disturbances patterns and forest management intensities.Specifically, we wanted to address the following research questions: 1.
How well can the major disturbance agents of Central European forests (i.e., windthrow, bark beetles, and harvest) be discriminated using Landsat-based agent attribution models?2.
Does incorporating information on the intra-annual timing of disturbances improve the attribution to change agents?3.
Does including prior knowledge about the occurrence of major storm events improve the attribution to disturbance agents?4.
Are the disturbance attribution models generalizable across study sites, i.e., are there differences in model performance and relative importance of predictor variables?

Study Area
Since forest management strongly impacts the manifestation of local disturbance patterns, we selected our study sites to include different levels of management intensity.Overall, our selected sites include three types of forest management zones: (1) strictly protected areas where natural disturbances were allowed to progress without human intervention, (2) protected areas with human intervention, and (3) forests actively managed for timber production.We selected three protected forest areas as study sites: The Bohemian Forest Ecosystem consisting of the Bavarian and Šumava National Parks in Germany and the Czech Republic, the Kalkalpen National Park in Austria and the core zone of the Tatra National Park in Slovakia (Figure 1).While the core zones of national park areas display non-intervention zones, low levels of human intervention are present in the surrounding management zones.Around each of the protected area sites, we constructed a 10-km buffer to include actively managed forests.Within the recent past, all study sites were affected by large-scale storm events (Table 1) and outbreaks of European spruce bark beetle (Ips typographus L.) [32,34,63].Disturbance regimes governed by these two natural agents are typical for Central European mountain forests [44,64].However, the manifestation of local disturbance patterns varies between study sites.Following multiple storm events in the 1980s and a series of hot summers in the 1990s, large-scale bark beetle outbreaks have led to extensive forest die-off in the Bohemian Forest Ecosystem [32,65,66].Moreover, the Kyrill storm in January 2007, as well as a thunderstorm event in July 2011, have caused substantial windthrows in high elevation forest stands [32].At the Tatra site, a storm event in November 2004 has caused massive windthrows on the southern flank of the mountain chain [67].The blowdown of large forest areas was followed by an outbreak of bark beetles, which primarily affected forest stands in the core zone of the national park [67,68].A second storm event in 2014 again led to large windthrows [69].Compared to the large-scale disturbances that occurred at the Bohemian Forest and Tatra sites, disturbances at the Kalkalpen site were much smaller in size.In recent years, blowdowns caused by the Kyrill storm in 2007 and two storm events in 2008 triggered a major bark beetle outbreak inside the national park area [34].Elevation ranges are ca.750-1450 m at the Bohemian Forest site, ca.400-2000 m at the Kalkalpen site, and ca.900-2650 m at the Tatra site.At all study sites, coniferous forests dominated by Norway spruce (Picea abies (L.) H. Karst) are the most prominent vegetation type in higher elevation zones (subalpine belt).At the Bohemian Forest and Kalkalpen site, lower elevations are dominated by European beech (Fagus sylvatica L.) forests, while the montane belt is formed by mixed forests of European beech, Norway spruce and Silver fir (Abies alba Mill.).At the Tatra and Kalkalpen sites, higher-elevation forests also feature European larch (Larix decidua Mill.), stone pine (Pinus cembra L.) and mountain pine (Pinus mugo Turra) as admixture species.
Within the recent past, all study sites were affected by large-scale storm events (Table 1) and outbreaks of European spruce bark beetle (Ips typographus L.) [32,34,63].Disturbance regimes governed by these two natural agents are typical for Central European mountain forests [44,64].However, the manifestation of local disturbance patterns varies between study sites.Following multiple storm events in the 1980s and a series of hot summers in the 1990s, large-scale bark beetle outbreaks have led to extensive forest die-off in the Bohemian Forest Ecosystem [32,65,66].Moreover, the Kyrill storm in January 2007, as well as a thunderstorm event in July 2011, have caused substantial windthrows in high elevation forest stands [32].At the Tatra site, a storm event in November 2004 has caused massive windthrows on the southern flank of the mountain chain [67].The blowdown of large forest areas was followed by an outbreak of bark beetles, which primarily affected forest stands in the core zone of the national park [67,68].A second storm event in 2014 again led to large windthrows [69].Compared to the large-scale disturbances that occurred at the Bohemian Forest and Tatra sites, disturbances at the Kalkalpen site were much smaller in size.In recent years, blowdowns caused by the Kyrill storm in 2007 and two storm events in 2008 triggered a major bark beetle outbreak inside the national park area [34].At all study sites, some of the disturbed forest stands located inside the national parks, particularly within the core zones, were left to progress without human intervention [34,45,67,70].Within the management zones of the national parks, most of the blowdown areas and bark-beetle-infested stands located at the park borders were salvage-logged in order to avoid an uncontrolled spread of bark beetles outside the park boundaries [34,45,67 [71].The images were converted to surface reflectance using the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) [72], except for Landsat OLI data, which was processed using the methods described in Vermote et al. [73].Clouds and cloud shadows were masked out using the Fmask algorithm [74].In a final preprocessing step, the Tasseled Cap (TC) indices were computed using the surface reflectance-derived transformation coefficients described in Crist [75] for all Landsat sensors.The different sensor characteristics of the Landsat 8 OLI satellite can lead to differences in derived vegetation indices (VIs) compared to previous Landsat sensors [76], thus potentially affecting the comparability of VI values when using both sensors in a time series framework.To test the impact of sensor differences on data continuity, we compared TC wetness values derived from OLI data with simulated ETM+ data, which were generated by applying the band equations from Roy et al. [76] to OLI reflectance data.We calculated TC wetness for 60,622 samples from the European Land Use and Cover Area-Frame Statistical Survey (LUCAS) [77] comprising different land cover types from across all of Europe.TC wetness calculated from OLI data tended to be lower than from simulated ETM+ data.However, we found that differences were considerably smaller in forested areas compared to other land cover types, which is in agreement with the findings of She et al. [78].As the observed differences generally were small (average difference in TC wetness <0.01 for forested areas) compared to changes associated with forest disturbance and OLI data made up only a small part of our analyzed time series, we did not apply any normalization to the OLI imagery.
The availability of cloud-free observations varied considerably between study sites.The high-elevation Tatra site was the cloudiest with only 3.8 cloud-free observations per year on average.The Bohemian Forest site had 9.3 cloud-free observations, and the Kalkalpen site had 7.8 cloud-free observations per year on average.

Reference Data Collection
For all study sites we collected two sets of training data: (1) a stratified random sample of Landsat pixels that was used for forest masking and disturbance detection; and (2) a sample of disturbance agent polygons that was used for creating classification models discriminating forest disturbance agents.
The sample of Landsat pixels consisted of 500 pixels at each study site.To include enough disturbed forest pixels, the sampling was stratified based on a map reflecting change magnitudes within the pixel time series [79], quantified as the negative change in the TC wetness index [75] associated with the greatest disturbance event in a time series.For every Landsat pixel in the sample, we determined whether it was forested at the beginning of the time series (in 1985) and whether (and when) a disturbance event occurred.The manual labeling process was based on visual interpretation of Landsat image chips and associated spectral-temporal trajectories, and high resolution imagery accessed through Google Earth (Google, Menlo ParK, CA, USA).Such an approach allows creating reliable reference datasets on forest dynamics in the absence of other ground truth data [80] and has been used by other studies [27,81,82].
The disturbance agent polygons were manually digitized on high-resolution imagery available in Google Earth and Landsat data.For every polygon, the disturbance agent was recorded using four disturbance agent classes: windthrow, cleared windthrow, bark beetles, and other harvest.The cleared windthrow class was assigned when wind-damaged stands were harvested (salvaged) within a year following the storm event.All other loggings of forest stands were ascribed the other harvest class.Thus, the class comprises regular harvests for timber production (clear-cut and partial harvests) as well as sanitation and preventative logging used to control bark beetle outbreaks.As sanitation logging in response to bark beetle disturbance is usually applied immediately after an attack is identified by ground monitoring and infested trees are removed, discriminating it from regular clear-cut harvests based on remote sensing data is extremely difficult [83].
In addition to the Landsat time series and high-resolution imagery, we used reference datasets available for the national park areas for interpreting disturbance agents at the Bohemian Forest and Kalkalpen sites.At the Bohemian Forest site, we used annual aerial survey data from the Bavarian Forest National Park going back to 1989.These surveys included mappings of disturbed forest stands with information on disturbance agent and post-disturbance management practices (i.e., whether a stand was salvage logged or not) [84].For the Kalkalpen National Park, we used a reference polygon dataset indicating which forest stands were disturbed since 1992, including information on the disturbance agent and year.
The sampling of disturbance agent polygons was done in two phases similar to the approach by Kennedy et al. [26].We used disturbance maps derived from the Landsat time series and high-resolution imagery to identify disturbed forest areas, which were then assigned with a label of the disturbance agent using all available reference data.In the first phase, we selected a purposive sample of disturbance areas which were then used to train preliminary agent attribution models.In the second phase, we extended our polygon sample by adding additional samples guided by the spatial predictions of the preliminary attribution models.The second sampling phase ensured that classes underrepresented in the preliminary sample were adequately represented in the final sample of disturbance agent polygons.Using model predictions to identify potential additional training samples can be more effective than manual search, even when some of the predictions turn out to be false [26].In total, we collected 478, 258, and 306 reference polygons for the Bohemian Forest, Kalkalpen, and Tatra sites, respectively (Table 2).The size of the digitized polygons varied between 1 and 1300 Landsat pixels, with an average polygon size of 23 pixels (ca.2.1 ha).

Disturbance Detection and Agent Attribution Approach
Our approach to mapping forest disturbances and disturbance agents consisted of three main steps: (1) forest mapping, (2) disturbance detection, and (3) disturbance agent attribution (Figure 2).

Disturbance Detection and Agent Attribution Approach
Our approach to mapping forest disturbances and disturbance agents consisted of three main steps: (1) forest mapping, (2) disturbance detection, and (3) disturbance agent attribution (Figure 2).First, areas forested at the beginning of the time series were delineated by creating a forest mask.During initial tests, we tried all three TC components (greenness, brightness, wetness) in change detection and agent attribution models, but found wetness to be most informative for detecting and separating disturbance agents at our study sites.Thus, we applied a breakpoint detection algorithm to TC wetness time series of all available, cloud-free Landsat observations.TC wetness has been shown to correlate with important tree and stand attributes [85,86] and was applied successfully in several studies on forest disturbance mapping and agent attribution [27,39,41,87].
From the detected breakpoints and trends, we extracted spectral and temporal metrics describing the greatest disturbance event, and then used these metrics as input to classify disturbed and undisturbed forest area.Finally, we classified the disturbed forest area according to one of four disturbance agent classes: windthrow, cleared windthrow, bark beetles, and other harvest.

Breakpoint Detection and Time Series Segmentation
To detect and characterize forest disturbances, we applied the breakpoint detection method employed by the Breaks For Additive Season and Trend (BFAST) algorithm [19,88] to the TC wetness time series.Our approach for detecting breakpoints is similar to the one used by de Jong et al. [89], but allows for detecting multiple breaks in a time series.
First, a linear model including harmonic terms was fitted to the pixel time series to account for effects of seasonality [52,90].Similar to the approach by Zhu et al. [20] for fitting models capturing seasonality and Bidirectional Reflectance Distribution Function (BRDF) effects in time series of all available Landsat data, we used a two-term harmonic (sine and cosine) model since higher-order functions tended to overfit.By accounting for the underlying inter-annual variation, seasonal patterns in the residuals unrelated to potential disturbances are removed [89].The resulting linear season-trend model is then partitioned into segments using the method described in Zeileis et al. [91].
A statistical test is applied to check whether structural breaks occur in the time series.We used the test based on moving sums (MOSUMs) of the residuals in ordinary least squares (OLS) regression models [88].If the test indicates evidence for a structural break (p < 0.05), the number and position of breakpoints are determined based on the method by Bai and Perron [92].Breakpoint First, areas forested at the beginning of the time series were delineated by creating a forest mask.During initial tests, we tried all three TC components (greenness, brightness, wetness) in change detection and agent attribution models, but found wetness to be most informative for detecting and separating disturbance agents at our study sites.Thus, we applied a breakpoint detection algorithm to TC wetness time series of all available, cloud-free Landsat observations.TC wetness has been shown to correlate with important tree and stand attributes [85,86] and was applied successfully in several studies on forest disturbance mapping and agent attribution [27,39,41,87].
From the detected breakpoints and trends, we extracted spectral and temporal metrics describing the greatest disturbance event, and then used these metrics as input to classify disturbed and undisturbed forest area.Finally, we classified the disturbed forest area according to one of four disturbance agent classes: windthrow, cleared windthrow, bark beetles, and other harvest.

Breakpoint Detection and Time Series Segmentation
To detect and characterize forest disturbances, we applied the breakpoint detection method employed by the Breaks For Additive Season and Trend (BFAST) algorithm [19,88] to the TC wetness time series.Our approach for detecting breakpoints is similar to the one used by de Jong et al. [89], but allows for detecting multiple breaks in a time series.
First, a linear model including harmonic terms was fitted to the pixel time series to account for effects of seasonality [52,90].Similar to the approach by Zhu et al. [20] for fitting models capturing seasonality and Bidirectional Reflectance Distribution Function (BRDF) effects in time series of all available Landsat data, we used a two-term harmonic (sine and cosine) model since higher-order functions tended to overfit.By accounting for the underlying inter-annual variation, seasonal patterns in the residuals unrelated to potential disturbances are removed [89].The resulting linear season-trend model is then partitioned into segments using the method described in Zeileis et al. [91].
A statistical test is applied to check whether structural breaks occur in the time series.We used the test based on moving sums (MOSUMs) of the residuals in ordinary least squares (OLS) regression models [88].If the test indicates evidence for a structural break (p < 0.05), the number and position of breakpoints are determined based on the method by Bai and Perron [92].Breakpoint positions Forests 2017, 8, 251 8 of 24 are determined by minimizing the residual sum of squares of linear segments, while the Bayesian Information Criterion (BIC) is used to select the optimal number of breaks [89].
As one parameter of the breakpoint algorithm, the bandwidth of the moving window used in the OLS-MOSUM test must be specified by the user (the bandwidth also determines the minimum length of resulting segments; h-parameter).Based on test runs, we set h to be 5% of the total data, reflecting a minimum segment length of ca.1.5 years.To exclude periods in which no breakpoints could be detected because of the minimum segment size defined by h, we removed the first and the last two years of the time series from further analysis, such that only events occurring between 1986 and 2014 were considered for forest disturbance detection and change agent attribution.The applied method is implemented in the strucchange package in R [93].

Computation of Metrics for Greatest Disturbance
The breakpoint detection algorithm results in a segmentation of the time series, in which OLS segments including linear trend terms and harmonic season terms are separated by breakpoints (i.e., discontinuities) [52,89].As harmonic terms were used to account for seasonality in the time series during the detection of breakpoints (i.e., to avoid the detection of spurious breakpoints due to seasonal effects), linear models without harmonic season terms were fitted to the time series segments to describe trend changes in the Landsat time series.We used these linear trends and the breakpoints separating segments to identify and characterize potential forest disturbances.To determine the major change event (i.e., greatest disturbance) within a pixel time series, the breakpoint or linear segment accounting for the largest negative change in TC wetness was identified.Thus, a disturbance was marked either by a breakpoint (abrupt change) or a linear segment (gradual change).
We estimated the onset date for each greatest disturbance event.Since a breakpoint detected by the time series segmentation method represents the last observation of a segment before the structural break occurs, the mean value between the date of the breakpoint and the next clear observation was taken as a naïve best guess of the actual disturbance date.By incorporating all clear observations in the breakpoint detection model, our approach allowed to derive within-year estimates of the disturbance date.Finally, a suite of metrics describing key characteristics of the spectral-temporal trajectory was extracted for the greatest disturbance (Figure 3, Table 3).positions are determined by minimizing the residual sum of squares of linear segments, while the Bayesian Information Criterion (BIC) is used to select the optimal number of breaks [89].
As one parameter of the breakpoint algorithm, the bandwidth of the moving window used in the OLS-MOSUM test must be specified by the user (the bandwidth also determines the minimum length of resulting segments; h-parameter).Based on test runs, we set h to be 5% of the total data, reflecting a minimum segment length of ca.1.5 years.To exclude periods in which no breakpoints could be detected because of the minimum segment size defined by h, we removed the first and the last two years of the time series from further analysis, such that only events occurring between 1986 and 2014 were considered for forest disturbance detection and change agent attribution.The applied method is implemented in the strucchange package in R [93].

Computation of Metrics for Greatest Disturbance
The breakpoint detection algorithm results in a segmentation of the time series, in which OLS segments including linear trend terms and harmonic season terms are separated by breakpoints (i.e., discontinuities) [52,89].As harmonic terms were used to account for seasonality in the time series during the detection of breakpoints (i.e., to avoid the detection of spurious breakpoints due to seasonal effects), linear models without harmonic season terms were fitted to the time series segments to describe trend changes in the Landsat time series.We used these linear trends and the breakpoints separating segments to identify and characterize potential forest disturbances.To determine the major change event (i.e., greatest disturbance) within a pixel time series, the breakpoint or linear segment accounting for the largest negative change in TC wetness was identified.Thus, a disturbance was marked either by a breakpoint (abrupt change) or a linear segment (gradual change).
We estimated the onset date for each greatest disturbance event.Since a breakpoint detected by the time series segmentation method represents the last observation of a segment before the structural break occurs, the mean value between the date of the breakpoint and the next clear observation was taken as a naïve best guess of the actual disturbance date.By incorporating all clear observations in the breakpoint detection model, our approach allowed to derive within-year estimates of the disturbance date.Finally, a suite of metrics describing key characteristics of the spectral-temporal trajectory was extracted for the greatest disturbance (Figure 3, Table 3).The extracted metrics can be divided into two groups: Spectral metrics describing the spectral change in the pixel time series associated with disturbances, and temporal metrics relating to the onset date and duration of the greatest disturbance (Table 3).While similar spectral metrics and disturbance duration have previously been used to attribute disturbance agents using Landsat time series [26,27,29,30], our approach introduces two additional temporal metrics which are based on the within-year estimates of disturbance dates.The seasonality metric (SEAS) reflects the within-year timing of a disturbance event, while the storm metric (STORM) relates the estimated disturbance date to prior knowledge about the occurrence of major storm events.To quantify the temporal proximity of disturbance date estimates and actual storm events, we used the number of clear observations between detected disturbances and known storm dates, assuming that most disturbances related to a storm event are detected within a few observations of the actual storm date.By using the number of clear Landsat observations instead of an absolute time difference, we wanted to ensure comparability between pixel locations with different numbers of clear observations.
In addition to the spectral and temporal metrics derived from the change detection algorithm, we also calculated context-based topography-related metrics for all study sites (Table 3).Such metrics, (elevation, slope, aspect and terrain roughness) were derived from the Shuttle Radar Topography Mission (SRTM) digital elevation model (DEM) with the same spatial resolution as the Landsat data (30 m).

Classification Models
We used random forest classifications [94] for forest masking, disturbance detection, and disturbance agent attribution.The parametrization of the random forest models was based on the recommended default values of the randomForest package for R [95].To create the forest mask, the sample of 500 interpreted Landsat pixels was used as training data.The classification was applied to a stack of pixel-wise metrics comprised of mean, median, minimum, maximum, 25% and 75% quantile of all TC indices (brightness, greenness, wetness) between 1984 and 1986.The classification model discriminating disturbed forest stands was trained using all interpreted Landsat pixels that fell inside forested areas according to the forest mask.That way, 376, 390, and 335 training pixels were used to train the classification models in the Bohemian Forest, Kalkalpen, and Tatra sites, respectively.We used the set of spectral metrics derived from the Landsat time series as input features for the disturbance classification (Table 3).
The training data used for building the disturbance agent classification models was derived by taking a sample of 250 pixels for every agent class (N = 1000) located inside the collected disturbance agent polygons.To account for the different sizes of the manually digitized training polygons, we used a stratified sampling approach that ensured the allocation of at least one sample in every polygon.The random forest models used for attributing disturbance agents was trained using spectral, temporal and topography-related metrics (Table 3) computed for disturbed forest pixels.
To evaluate our approach, we built a range of different attribution models.To assess the impact of adding a group of metrics on model performance, we fitted models including only the spectral metrics, models including spectral and temporal metrics, models including spectral and topography-related metrics, and full models using all metrics.Models were trained site-wise; however, to investigate whether predictor variables had similar effects across study sites, we also fitted a global model using the training samples from all study sites as input data.To assess the relative importance of predictor variables, we computed the mean decrease in the Gini coefficient [94] and normalized the values to a 0-100 scale of relative importance scores.The performance of the classification models was assessed using the out-of-bag (OOB) classification errors returned by the random forest algorithm [96].
We used the spatial predictions of the attribution models to derive trends of forest disturbance agents for each study site.Agent-specific forest disturbance rates were calculated using the attribution results and the forest masks.

Forest Masking and Disturbance Detection
According to the OOB estimates, the forest masks were highly reliable, with overall accuracies of 95%, 90%, and 88% for the Bohemian Forest, Kalkalpen, and Tatra sites, respectively.Also, the accuracies of the disturbance detection models were high, with overall accuracies of 97%, 91%, and 92% at the Bohemian Forest, Kalkalpen, and Tatra sites, respectively (Table 4).We assessed the temporal accuracy of disturbance detection by comparing the disturbance year estimates obtained from visually interpreting the sample of Landsat pixels with disturbance dates assigned by the change detection algorithm (Figure 4).Overall, we found good agreement between the Landsat-based and interactively derived disturbance year estimates, with 81% of disturbances in the sample pixels being detected within +/−1 year of the interpreted disturbance year.We assessed the temporal accuracy of disturbance detection by comparing the disturbance year estimates obtained from visually interpreting the sample of Landsat pixels with disturbance dates assigned by the change detection algorithm (Figure 4).Overall, we found good agreement between the Landsat-based and interactively derived disturbance year estimates, with 81% of disturbances in the sample pixels being detected within +/−1 year of the interpreted disturbance year.

Disturbance Agent Attribution
The OOB error indicated high accuracies of the disturbance agent models.The full models including all predictor metrics had overall accuracies of 78%, 82%, and 91% for the Bohemian Forest, Kalkalpen, and Tatra sites, respectively (Figure 5).Class errors were relatively well balanced, with producer's accuracies (1-omission error) ranging from 71-80%, 70-87%, and 84-96% and Users' accuracies (1-commission error) ranging from 75-81%, 79-83%, and 88-94% at the Bavarian Forest, Kalkalpen, and Tatra sites, respectively.Models built with only the spectral metrics showed considerably lower class and overall accuracies than models incorporating the temporal and/or

Disturbance Agent Attribution
The OOB error indicated high accuracies of the disturbance agent models.The full models including all predictor metrics had overall accuracies of 78%, 82%, and 91% for the Bohemian Forest, Kalkalpen, and Tatra sites, respectively (Figure 5).Class errors were relatively well balanced, with producer's accuracies (1-omission error) ranging from 71-80%, 70-87%, and 84-96% and Users' accuracies (1-commission error) ranging from 75-81%, 79-83%, and 88-94% at the Bavarian Forest, Kalkalpen, and Tatra sites, respectively.Models built with only the spectral metrics showed considerably lower class and overall accuracies than models incorporating the temporal and/or topography-related metrics.When compared to the model only using the spectral metrics, including the temporal metrics resulted in increases of overall accuracies of 10%, 10%, and 7%.Including the topography-related metrics increased overall accuracies by 8%, 18%, and 10%.Including both the sets of temporal and topography-related metrics increased overall accuracies by 12%, 21%, and 15%, respectively.topography-related metrics.When compared to the model only using the spectral metrics, including the temporal metrics resulted in increases of overall accuracies of 10%, 10%, and 7%.Including the topography-related metrics increased overall accuracies by 8%, 18%, and 10%.Including both the sets of temporal and topography-related metrics increased overall accuracies by 12%, 21%, and 15%, respectively.

Importance and Distributions of Predictor Metrics
The variable importance scores based on the Gini-index were used to assess the importance of the individual metrics for model performance.According to the importance scores, the elevation variable (ELE) was the most important predictor of the disturbance agent in all site-specific models as well as the global model (Figure 6).The temporal metric relating detected disturbance to prior knowledge about the occurrence of major storm events also showed some of the highest variable importance across models, being ranked second most important variable in the global model.Only in the model for the Bohemian Forest Ecosystem, a spectral metric (disturbance magnitude) was among the two most important predictors.Disturbance magnitude (MAG) and post-disturbance

Importance and Distributions of Predictor Metrics
The variable importance scores based on the Gini-index were used to assess the importance of the individual metrics for model performance.According to the importance scores, the elevation variable (ELE) was the most important predictor of the disturbance agent in all site-specific models as well as the global model (Figure 6).The temporal metric relating detected disturbance to prior knowledge about the occurrence of major storm events also showed some of the highest variable importance across models, being ranked second most important variable in the global model.Only in the model for the Bohemian Forest Ecosystem, a spectral metric (disturbance magnitude) was among the two most important predictors.Disturbance magnitude (MAG) and post-disturbance wetness change (POSTCH) showed the highest variable importance scores for spectral metrics across all models.The temporal seasonality metric, quantifying the intra-annual timing of disturbances, was a less important predictor, being ranked seventh most important in the global model.
wetness change (POSTCH) showed the highest variable importance scores for spectral metrics across all models.The temporal seasonality metric, quantifying the intra-annual timing of disturbances, was a less important predictor, being ranked seventh most important in the global model.To investigate differences in the metrics between disturbance agent classes, we analyzed the distributions for the disturbance agent samples.Bark beetle disturbances had the lowest change magnitudes from all disturbance agents (Figure 7).In contrast, other harvests and cleared windthrows showed the highest disturbance magnitudes across all study sites.Moreover, both anthropogenic disturbance agent classes (cleared windthrow and other harvest) were associated with higher increases in TC wetness following the disturbance event than unmanaged natural disturbances by windthrow and bark beetles.To investigate differences in the metrics between disturbance agent classes, we analyzed the distributions for the disturbance agent samples.Bark beetle disturbances had the lowest change magnitudes from all disturbance agents (Figure 7).In contrast, other harvests and cleared windthrows showed the highest disturbance magnitudes across all study sites.Moreover, both anthropogenic disturbance agent classes (cleared windthrow and other harvest) were associated with higher increases in TC wetness following the disturbance event than unmanaged natural disturbances by windthrow and bark beetles.
Thus, windthrows cleared by salvage logging featured spectral change signals very different from unmanaged blowdown areas, which in direct comparison had lower disturbance magnitudes and post-disturbance wetness increases.These strong differences between spectral signatures of unmanaged windthrows and cleared windthrows demonstrate the impact of post-disturbance forest-management on the change signals recorded by satellite sensors.At the same time, cleared windthrows were generally more similar to regular harvest disturbances in their spectral behavior.In contrast to the similar spectral signatures observed for cleared windthrows and other harvest, we found significant differences in the distributions of topography-related and temporal metrics for these agents.Generally, disturbance agents featured distinct elevational distributions at all study sites, explaining the variables' effectiveness for change agent attribution (Figure 8).In particular, windthrow and bark beetle disturbances occurred more frequently at higher elevations than harvests.
Disturbances by windthrow and cleared windthrow were detected very close to dates of known storm events, while bark beetles and other harvest disturbances were much less clustered in their distributions of the metric (Figure 8).Integrating data from all study sites, the majority (75%) of windthrow samples featured disturbances within +/−3 clear observations of actual storm event dates.However, the detection offset for windthrow disturbances varied between study sites.Whereas 91% of the windthrow samples showed disturbances within +/−3 clear observations at the Bohemian Forest Ecosystem, this number was considerably lower at the Kalkalpen (72%) and Tatra (64%) sites.Thus, comparing the study sites, windthrow disturbances were detected closer to the actual storm dates with an increasing number of available clear observations.In contrast to the similar spectral signatures observed for cleared windthrows and other harvest, we found significant differences in the distributions of topography-related and temporal metrics for these agents.Generally, disturbance agents featured distinct elevational distributions at all study sites, explaining the variables' effectiveness for change agent attribution (Figure 8).In particular, windthrow and bark beetle disturbances occurred more frequently at higher elevations than harvests.

Spatial Model Predictions and Disturbance Trends
We created disturbance agent maps showing the spatial distribution of the four disturbance agent classes at the three study sites and compared them with available reference datasets (Figure 9).The spatial patterns of the Landsat-based disturbance agent maps were in good agreement with aerial survey mappings based on visual interpretation of orthophotos.Compared to aerial survey interpretations, the pixel-wise classification of disturbance agents based on Landsat data resulted in more heterogeneous agent patterns.Disturbances by windthrow and cleared windthrow were detected very close to dates of known storm events, while bark beetles and other harvest disturbances were much less clustered in their distributions of the metric (Figure 8).Integrating data from all study sites, the majority (75%) of windthrow samples featured disturbances within +/−3 clear observations of actual storm event dates.However, the detection offset for windthrow disturbances varied between study sites.Whereas 91% of the windthrow samples showed disturbances within +/−3 clear observations at the Bohemian Forest Ecosystem, this number was considerably lower at the Kalkalpen (72%) and Tatra (64%) sites.Thus, comparing the study sites, windthrow disturbances were detected closer to the actual storm dates with an increasing number of available clear observations.

Spatial Model Predictions and Disturbance Trends
We created disturbance agent maps showing the spatial distribution of the four disturbance agent classes at the three study sites and compared them with available reference datasets (Figure 9).The spatial patterns of the Landsat-based disturbance agent maps were in good agreement with aerial survey mappings based on visual interpretation of orthophotos.Compared to aerial survey interpretations, the pixel-wise classification of disturbance agents based on Landsat data resulted in more heterogeneous agent patterns.The disturbance rates derived for the agent classes were highly variable in space and time (Figure 10).In particular, disturbances associated with windstorm events (windthrows and cleared windthrows) were temporally clustered, appearing as spikes in the annual disturbance rates.At the Kalkalpen and Tatra sites, disturbance rates by bark beetles increased following the respective storm events, while at the Bohemian Forest Ecosystem the earlier outbreak cycle in the 1990s preceded the major storm events of 2007 and 2011.Other harvests were the most constant cause of forest disturbance at all study sites.Disturbance rates by harvesting activities increased with higher levels of natural disturbances by bark beetles and windthrow.The disturbance rates derived for the agent classes were highly variable in space and time (Figure 10).In particular, disturbances associated with windstorm events (windthrows and cleared windthrows) were temporally clustered, appearing as spikes in the annual disturbance rates.At the Kalkalpen and Tatra sites, disturbance rates by bark beetles increased following the respective storm events, while at the Bohemian Forest Ecosystem the earlier outbreak cycle in the 1990s preceded the major storm events of 2007 and 2011.Other harvests were the most constant cause of forest disturbance at all study sites.Disturbance rates by harvesting activities increased with higher levels of natural disturbances by bark beetles and windthrow.

Disturbance Agent Attribution
The high accuracies of the disturbance agent models across all test sites indicate that disturbances by windthrow, bark beetles, and harvest can be discriminated using classification models integrating data from dense Landsat time series, topography-related metrics, and prior information about the occurrence of disturbance events.Moreover, we were able to separate cleared windthrows from unmanaged windthrows in our classification models, thus providing valuable information on the forest management response to storm-related disturbances.In general, thanks to the possibility to investigate large areas and the temporal depth provided by the Landsat archive, the integrated Landsat-based mapping of windthrow and bark beetle disturbances offers great potential to improve the understanding of interactions between these two disturbance agents [12,34].
Comparing agent attribution models among study sites revealed several metrics that were consistently important across sites: elevation, the temporal proximity to storm events, the sum of TC wetness change after a disturbance, and TC wetness change magnitude.The disturbance agent maps and derived forest disturbance rates resembled spatial and temporal patterns known from previous studies and aerial overview surveys [32,34,67].For instance, the model predictions reproduced the different outbreak patterns of I. typographus at the study locations.At the Kalkalpen and Tatra sites, bark beetle outbreaks followed as a direct consequence of storm events [27,53], while in the Bohemian Forest Ecosystem, the largest bark beetle outbreak during the 1990s and early 2000s was largely climate-driven and not directly preceded by large-scale windthrows [66].The synchronicity of increased disturbance by other harvests and natural agents (i.e., bark beetles and windthrow)

Disturbance Agent Attribution
The high accuracies of the disturbance agent models across all test sites indicate that disturbances by windthrow, bark beetles, and harvest can be discriminated using classification models integrating data from dense Landsat time series, topography-related metrics, and prior information about the occurrence of disturbance events.Moreover, we were able to separate cleared windthrows from unmanaged windthrows in our classification models, thus providing valuable information on the forest management response to storm-related disturbances.In general, thanks to the possibility to investigate large areas and the temporal depth provided by the Landsat archive, the integrated Landsat-based mapping of windthrow and bark beetle disturbances offers great potential to improve the understanding of interactions between these two disturbance agents [12,34].
Comparing agent attribution models among study sites revealed several metrics that were consistently important across sites: elevation, the temporal proximity to storm events, the sum of TC wetness change after a disturbance, and TC wetness change magnitude.The disturbance agent maps and derived forest disturbance rates resembled spatial and temporal patterns known from previous studies and aerial overview surveys [32,34,67].For instance, the model predictions reproduced the different outbreak patterns of I. typographus at the study locations.At the Kalkalpen and Tatra sites, bark beetle outbreaks followed as a direct consequence of storm events [27,53], while in the Bohemian Forest Ecosystem, the largest bark beetle outbreak during the 1990s and early 2000s was largely climate-driven and not directly preceded by large-scale windthrows [66].The synchronicity of increased disturbance by other harvests and natural agents (i.e., bark beetles and windthrow) observable in the agent-specific disturbance rates can likely be attributed to increased levels of forest management in the form of preventative felling, sanitation and salvage logging in reaction to natural disturbances.
An alternative to attributing disturbance agents on the level of individual pixels is identifying forest disturbance causes on the level of disturbance patches.Previous studies have demonstrated that spatial information can help distinguishing between disturbance agents, as different disturbance types have been shown to result in characteristic spatial patterns [26,30,50].Generally, with forest disturbances typically appearing as disturbance patches, patch-based attribution approaches can be considered to have conceptual advantages over pixel-based approaches [26].However, the definition of appropriate decision heuristics on which pixels are aggregated is critical to derive meaningful patches reflecting underlying disturbance events [81].The fine scale and the strong spatial and temporal interactions of natural and anthropogenic disturbances in European forests [35] pose a challenge for the spatial segmentation of disturbance patches.

Usefulness of Temporal Metrics for Disturbance Agent Attribution
Applying the change detection algorithm to all clear Landsat observations allowed determining the date of forest disturbances on an intra-annual temporal level.Disturbance date estimates for most windthrow disturbances, for which the actual disturbance dates (i.e., the dates of major windstorm events) were known, showed margins of error of only a few observations.This demonstrates the capability of the change detection algorithm to accurately identify disturbance dates in intra-annual time series comprising all available Landsat observations.The accuracy of the estimated disturbance date increased with data availability, indicating that the approach benefits from frequent observations.
Including prior knowledge about windstorm events substantially improved the discrimination of forest disturbance agents, with the temporal proximity to storm events being one of the most informative predictors in the agent attribution models across all test sites.This demonstrates the usefulness of including external information about storm event dates to discriminate wind-related forest disturbances when mapping disturbance agents based on Landsat time series.In contrast to the importance of the proximity to storm events, we found that the within-year timing of disturbances (seasonality metric) was less informative for predicting disturbance agents.The fact that storm events have occurred in both summer and winter months at two test sites (the Bohemian Forest and Tatra sites) likely contributed to the relatively low variable importance of the seasonality metric.While on a continental scale in Europe a very distinct seasonality of windstorm events can be observed [53], differing temporal patterns on smaller spatial scales can make the within-year timing of disturbances less meaningful as a predictor.The number of available Landsat observations likely also compromises the effectiveness of the seasonality metric.With only a few clear observations per year available at many sites for large parts of the Landsat data record, the potential to differentiate between seasonal patterns of forest disturbance agents will therefore vary a lot depending on the region of interest.The effectiveness of this metric may increase when data fusion approaches, integrating data from other sensors such as MODIS, are used [51].With increasing availability of data from Landsat and Landsat-like sensors, such as the recently launched Sentinel 2A/B satellites, the potential for monitoring forest change on an intra-annual temporal scale will improve [20].

Importance of Other Metrics and Transferability to Other Study Regions
The usefulness of spectral change magnitudes to discriminate between disturbance agents observed in our models is in line with the results of other studies [27,30,31].The disturbance magnitudes associated with natural agents (bark beetles and windthrow) were lower than for cleared windthrow and other harvest.These differences are likely caused by varying mortality rates during bark beetle attacks (i.e., non-host tress or trees with high defensive capacity), and residual trees that were not uprooted during storms.Further, changes in the canopy caused by natural disturbance agents exposes understory vegetation and natural regeneration, which further can dampen the disturbance magnitude [39].These differences in change magnitudes are also mirrored in the different post-disturbance spectral changes of disturbance agents, with higher increases of TC wetness following disturbances by forest management (cleared windthrow and other harvest classes), thus making this metric effective for distinguishing different types of forest disturbance.
Moreover, we found topography-related information, particularly elevation, to be very helpful for discriminating disturbance agents.The effectiveness of the elevation variable in our models can be related to a number of factors.First, topography-related site conditions can impact the predisposition of a forest stand towards natural disturbances.For example, wind speed tends to increase with altitude, which in turn increases the likelihood of windthrow disturbances [97].Moreover, at the mountain forests sites investigated in this study, elevation is strongly correlated with the distribution of different forest types and species, which can show varying levels of susceptibility towards different natural disturbance agents.Spruce-dominated forests are prevailing in the higher altitude zones of the study sites.Norway spruce is both the most important host tree for I. typographus [35,44,98] and has been shown to be particularly vulnerable towards windthrow [99][100][101].
However, the elevation variable also reflects a clear elevational gradient of different forest management zones present at all study sites.This effect associated with differences in forest management was probably even more important for the effectiveness of the elevation metric in our agent attribution models.Since protected forest areas with no or limited human intervention are confined to high-elevation areas and natural disturbances are actively managed in the surrounding forests, natural disturbances tended to occur at higher elevations than harvests.
Outside the strictly protected zones of the national parks, the active management of natural disturbances through salvage and sanitation logging decreased the frequency of bark beetle disturbances, and blowdown areas were mostly cleared.Since cleared windthrows exhibited spectral signatures very similar to other harvests, the separation of these two agent classes strongly relied on topography-related and temporal information at our study sites, particularly the elevation and the temporal proximity to storm events.Combining spectral with temporal information derived from intra-annual time series as well as prior information on the occurrence of storm events thus holds great potential for mapping wind disturbances and related salvage logging activities based on Landsat data.This result has important implications for large-area applications of remote sensing-based disturbance agent monitoring: As salvage logging of windthrow disturbances is the norm in most Central European forests [45,102], the discrimination of salvage logging from other harvesting activities is important to accurately reflect the impact of windthrow disturbances on forest ecosystems.Given the lack of standardized reporting systems for forest storm damage and its management in Europe [103], remote sensing-based maps of storm-related salvage logging can provide important insights into associated spatiotemporal patterns and help assessing the effectiveness of management activities [43,67].Moreover, as salvage logging has a range of important ecological consequences [102] and its use is critically discussed [45,104], remote-sensing based datasets on storm-related salvage operations can help more accurately assessing impacts on forest ecosystems.
The importance of topographic information and prior knowledge about windstorm events demonstrates that remote sensing-based mapping of forest disturbance agents can strongly benefit from including additional sources of information that complement the spectral-temporal information contained in satellite imagery [26].Including datasets on forest attributes, such as forest type and management, likely could further improve the discrimination of disturbance agents.Previous studies already have demonstrated the potential to incorporate prior information about disturbance events by including ancillary spatial datasets into remote sensing-based analyses of forest disturbance [31,[60][61][62].However, the transferability of approaches relying on external sources of information strongly depends on the availability of appropriate reference datasets.While detailed spatial datasets on forest disturbance agents are often not available over larger geographic extents for European forests, information on the timing of major natural disturbances events, such as windstorms or fires, can typically be obtained at regional levels from databases or information systems (e.g., the FORESTORM database: http://www.efiatlantic.efi.int/portal/databases/forestorms/).Because of the general beetles, and harvest disturbances showed similar spectral, topographic, and temporal characteristics across study sites, our approach is likely transferable to other European forests with similar species composition and disturbance regimes.

Conclusions
Thanks to recent methodological advances, Landsat-based forest monitoring can provide otherwise often unavailable spatially explicit information on forest disturbance agents that offers important insights for research and forest management.With the increasing availability of analysis-ready data from sensors such as Landsat, MODIS and Sentinel-2, as well as improved methods and computation power to handle dense time series of satellite data over large areas, monitoring forest disturbances on an intra-annual, temporal scale becomes feasible.As the results of this study demonstrated, monitoring forest disturbances on an intra-annual level not only allows for accurate delineation of disturbance dates and trends, but can also improve the identification of underlying disturbance agents.In this way, we have been able to reliably discriminate the major forest disturbance agents of Central Europe.
Extending the methods for attributing disturbance agents based on Landsat data to larger geographic extents could allow for mapping disturbance agents of European forests on a regional to continental scale.Such maps would greatly improve our understanding of spatiotemporal trends in forest disturbances.Moreover, such large-scale spatially explicit datasets on forest disturbance agents would help with gaining new insights into drivers and interactions of disturbance agents, thereby assisting with the sustainable management of European forest ecosystems.

Figure 1 .
Figure 1.Location and extent of the three study sites.Dark gray areas represent initial protected area polygons and light gray areas the constructed buffer zones.Elevation ranges are ca.750-1450 m at the Bohemian Forest site, ca.400-2000 m at the Kalkalpen site, and ca.900-2650 m at the Tatra site.At all study sites, coniferous forests dominated by Norway spruce (Picea abies (L.) H. Karst) are the most prominent vegetation type in higher elevation zones (subalpine belt).At the Bohemian Forest and Kalkalpen site, lower elevations are dominated by European beech (Fagus sylvatica L.) forests, while the montane belt is formed by mixed forests of European beech, Norway spruce and Silver fir (Abies alba Mill.).At the Tatra and Kalkalpen sites, higher-elevation forests also feature European larch (Larix decidua Mill.), stone pine (Pinus cembra L.) and mountain pine (Pinus mugo Turra) as admixture species.Within the recent past, all study sites were affected by large-scale storm events (Table1) and outbreaks of European spruce bark beetle (Ips typographus L.)[32,34,63].Disturbance regimes governed by these two natural agents are typical for Central European mountain forests[44,64].However, the manifestation of local disturbance patterns varies between study sites.Following multiple storm events in the 1980s and a series of hot summers in the 1990s, large-scale bark beetle outbreaks have led to extensive forest die-off in the Bohemian Forest Ecosystem[32,65,66].Moreover, the Kyrill storm in January 2007, as well as a thunderstorm event in July 2011, have caused substantial windthrows in high elevation forest stands[32].At the Tatra site, a storm event in November 2004 has caused massive windthrows on the southern flank of the mountain chain[67].The blowdown of large forest areas was followed by an outbreak of bark beetles, which primarily affected forest stands in the core zone of the national park[67,68].A second storm event in 2014 again led to large windthrows[69].Compared to the large-scale disturbances that occurred at the Bohemian Forest and Tatra sites, disturbances at the Kalkalpen site were much smaller in size.In recent years, blowdowns caused by the Kyrill storm in 2007 and two storm events in 2008 triggered a major bark beetle outbreak inside the national park area[34].

Figure 1 .
Figure 1.Location and extent of the three study sites.Dark gray areas represent initial protected area polygons and light gray areas the constructed buffer zones.

Figure 2 .
Figure 2. Workflow used for detecting and attributing forest disturbances.

Figure 2 .
Figure 2. Workflow used for detecting and attributing forest disturbances.

Figure 3 .
Figure 3. Schematic illustration of metric extraction from linear model segments with the greatest disturbances highlighted in red.Two cases are shown: (A) an abrupt disturbance process detected as a breakpoint and (B) a gradual disturbance detected as a linear segment.

Figure 3 .Table 3 .
Figure 3. Schematic illustration of metric extraction from linear model segments with the greatest disturbances highlighted in red.Two cases are shown: (A) an abrupt disturbance process detected as a breakpoint and (B) a gradual disturbance detected as a linear segment.

Figure 4 .
Figure 4. Scatterplots of disturbance year estimates obtained from visual interpretation (interpreted) and disturbance years obtained from the breakpoint detection algorithm (mapped).Diagonal lines represent complete agreement.Data points are plotted with 75% transparency.

Figure 4 .
Figure 4. Scatterplots of disturbance year estimates obtained from visual interpretation (interpreted) and disturbance years obtained from the breakpoint detection algorithm (mapped).Diagonal lines represent complete agreement.Data points are plotted with 75% transparency.

Figure 5 .
Figure 5. Heatmap of model accuracies including the different sets of metrics for the three site-specific models and the global model.PA = Producer's accuracy, UA = User's accuracy, OA = Overall accuracy.Spec.= Spectral, Temp.= Temporal, Topo.= Topography-related.

Figure 5 .
Figure 5. Heatmap of model accuracies including the different sets of metrics for the three site-specific models and the global model.PA = Producer's accuracy, UA = User's accuracy, OA = Overall accuracy.Spec.= Spectral, Temp.= Temporal, Topo.= Topography-related.

Figure 6 .
Figure 6.Relative variable importance for site-specific models and the combined global model as normalized importance scores based on the Gini index.Colors represent different metric groups (spectral, topography-related, temporal).

Figure 6 .
Figure 6.Relative variable importance for site-specific models and the combined global model as normalized importance scores based on the Gini index.Colors represent different metric groups (spectral, topography-related, temporal).

Figure 7 .
Figure 7. Distributions of disturbance magnitude and post-disturbance change for disturbance agent classes based on visually interpreted agent polygons.TC = Tasseled Cap.Thus, windthrows cleared by salvage logging featured spectral change signals very different from unmanaged blowdown areas, which in direct comparison had lower disturbance magnitudes and post-disturbance wetness increases.These strong differences between spectral signatures of unmanaged windthrows and cleared windthrows demonstrate the impact of post-disturbance forest-management on the change signals recorded by satellite sensors.At the same time, cleared windthrows were generally more similar to regular harvest disturbances in their spectral behavior.In contrast to the similar spectral signatures observed for cleared windthrows and other harvest, we found significant differences in the distributions of topography-related and temporal metrics for these agents.Generally, disturbance agents featured distinct elevational distributions at all study sites, explaining the variables' effectiveness for change agent attribution (Figure8).In particular, windthrow and bark beetle disturbances occurred more frequently at higher elevations than harvests.Disturbances by windthrow and cleared windthrow were detected very close to dates of known storm events, while bark beetles and other harvest disturbances were much less clustered in their distributions of the metric (Figure8).Integrating data from all study sites, the majority (75%) of windthrow samples featured disturbances within +/−3 clear observations of actual storm event dates.However, the detection offset for windthrow disturbances varied between study sites.Whereas 91% of the windthrow samples showed disturbances within +/−3 clear observations at the Bohemian Forest Ecosystem, this number was considerably lower at the Kalkalpen (72%) and Tatra (64%) sites.Thus, comparing the study sites, windthrow disturbances were detected closer to the actual storm dates with an increasing number of available clear observations.

Figure 7 .
Figure 7. Distributions of disturbance magnitude and post-disturbance change for disturbance agent classes based on visually interpreted agent polygons.TC = Tasseled Cap.

Forests 2017, 8 , 251 15 of 25 Figure 8 .
Figure 8. Distributions of elevation and proximity to storm events (storm metric) for disturbance agent classes based on visually interpreted agent polygons.

Figure 8 .
Figure 8. Distributions of elevation and proximity to storm events (storm metric) for disturbance agent classes based on visually interpreted agent polygons.

Figure 9 .
Figure 9.Comparison of (A) Landsat-based disturbance agent map and (B) aerial survey data for a section of the Bavarian Forest National Park.Aerial survey polygons were rasterized at 30 m resolution for comparison with Landsat-based predictions.The disturbances from other harvesting inside the national park were sanitation logging to control bark beetles.

Figure 9 .
Figure 9.Comparison of (A) Landsat-based disturbance agent map and (B) aerial survey data for a section of the Bavarian Forest National Park.Aerial survey polygons were rasterized at 30 m resolution for comparison with Landsat-based predictions.The disturbances from other harvesting inside the national park were sanitation logging to control bark beetles.

Figure 10 .
Figure 10.Annual disturbance rates for the four disturbance agent classes at the study sites.

Figure 10 .
Figure 10.Annual disturbance rates for the four disturbance agent classes at the study sites.

Table 1 .
Overview of major storm events that affected the study sites within the observation period in which we mapped forest disturbances.
]. ESA; https://earth.esa.int/)archives.The ESA images were geometrically aligned with the USGS imagery using the automatic georectification algorithm developed by Gao et al.

Table 2 .
Number of polygon samples collected for the disturbance agent classes at the study sites.

Table 4 .
Overview of out-of-bag (OOB) model accuracies of forest masks and disturbance maps created for the study sites.