The Significance of Tree Height as a Predictor of Tree Mortality during Bark Beetle Outbreaks in a Small Catchment

Abstract

Bark beetle outbreaks damage forests and kill trees worldwide, but many aspects of their dynamics remain unexplained. Our aim was to identify predictors for individual tree deaths within the small (0.7 km²) Plešné Lake catchment in the Šumava National Park in southwestern Czechia. Within this area, >60,000 trees were geo-referenced and categorized from ten aerial images (20 cm spatial resolution) between 2000 and 2015. For each year for which aerial images were available, we calculated tree densities of different categories and diameters. Tree height was evaluated by means of LiDAR in two terrestrial campaigns (2010 and 2011). A machine learning technique was then used to evaluate the most important variables. The resulting relationships were largely nonlinear and differed among years; however, individual trait tree height proved to be the most influential variable in each year. Higher trees were more likely to have died during either the undisturbed phase (2000 and 2003), the disturbed phase (2005–2011), or the recovery phase (2013). Our results indicate that salvage logging may not be the most effective measure for protecting trees in small catchments.

1. Introduction

Bark beetle outbreaks damage forests worldwide, but recent large-scale bark beetle outbreaks have shown patterns that are heterogeneous on several scales. Although recent empirical [1] and modelling studies [2,3,4] have significantly advanced our understanding, knowledge gaps remain. For example, processes leading to the collapse of bark beetle outbreaks occur on scales from less than three meters to several kilometers, but whether the same factors are decisive on all scales is not yet known. Effects from these different scales might act in the same or opposite directions [5,6]. Past studies either focused on small patches where all trees were characterized, or on larger areas where trees were only roughly characterized. Our objective was to provide a more detailed picture and follow forest damage on an individual tree basis throughout an entire catchment.

1.1. Bark Beetle Behavior—Endemic and Epidemic Phases

Depending on the bark beetle species, bark beetles behave differently in endemic and epidemic phases. In endemic phases, most bark beetles only colonize dead trees [7,8] and are attracted through the volatiles such trees exude [9]. These dead trees then become nuclei for further infestations. In parallel, storm events lead to trees first being uprooted, then dying, and exuding volatiles. In this way, storm events often trigger bark beetle attacks of damaged forests. Most bark beetles’ behavior in the endemic phases is better described by primary attraction by such dead trees than by random orientation [10]. In contrast, in epidemic phases, when the density of dead and damaged trees increases, e.g., after wind fells, bark beetles may enter epidemic phases and start “mass attacks” [11]: they attract congeners to the same tree as well as to healthy trees by emitting aggregation pheromones [7,8]. A positive relationship between the volume of dead broken trees, e.g., by wind and storms, and the proportion of standing trees being infested has been described [12]. While healthy trees with low susceptibility are generally not infested in endemic phases, they may be infested during epidemic phases, when there are mass attacks by large populations [13]. Ips typographus, the species responsible for most spruce tree damage in Central Europe [14], is most often univoltine (i.e., the beetles have only one generation, and fly to the next host only once per year) and thus, each female infests only one tree per year. Thus, the propagation of a bark beetle attack is spatially limited in time in both endemic and epidemic phases. Due to climate change, however, an increasing number of multivoltine populations have been observed [15].

1.2. Spatial Spread of Bark Beetles

The spatial spread of bark beetles has been studied in different areas and with differing focus. For example, univariate statistical models have been constructed based on the mapping of storm damage and subsequent bark beetle damage, assessing damage in a small area over time. One study found a clear relationship between tree damage and infested trees, but the smallest distance tested was 100 m [15]. In another study [7], the colonization of bark-beetle-damaged trees was correlated with storm gaps within 2000 m (distances between 500 m and 2000 m were tested). Yet another study found that 66% of attacks occurred within 100 m of a source, and 100% within 500 m [16]. These results are in line with a summary of earlier studies [17], showing that only about a third of a population is caught by pheromone traps while the other two thirds travel farther. Wichmann and Ravn [18] found that in an endemic situation, 50% of a local population seemed to propagate as far as 500 m. In contrast, in an epidemic situation, the same study found 80% of new attacks occurring within 250 m of an older attack [18]. However, distances smaller than 100 m have rarely been tested, and it is unknown whether there are also differences between epidemic and endemic phases at scales < 100 m. Remote sensing has been employed to derive data on individual trees to help address this issue, but databases are generally limited to just a couple thousand trees [19]. In the catchment studied here, a data set of >60,000 trees, observed over 10 distinct years [20], allowed us to tackle the question of the importance of small scale while comparing a larger number of trees at a wide range of distances.

1.3. Scales in Bark Beetle Attacks

When studying insect outbreaks, approximately four scales have been considered important. Conventionally, the regional or catchment scale encompasses one to several square kilometers [21]. The stand scale usually is formulated to act within about a 200 m radius, representing a homogeneous tree composition [22,23]. The plot or site scale is smaller, and characterized by small-scale elevation, slope, and microclimate conditions, and encompasses less than a 50 m radius [24]. Lastly, some individual properties such as diameter at breast height (DBH), tree health, and tree height, are also relevant. While some individual tree properties can be derived via remote sensing, there are several limitations to the various methods. Very high-resolution (VHR) remote sensing data such as IKONOS (or aerial cameras with very high-resolution sensors) can be used in classification tools such as random forests, and can be trained to recognize individual trees and even species. However, their calibration requires extensive reference data that are rarely available. Thus, they are not generally calibrated [25], and there are only few published accounts of using VHR for evaluating tree state. Thus, individual tree health is usually not studied beyond the scale of plots of a few square meters, because larger scales would encompass immense time requirements.

Many factors on the large, i.e., regional or catchment, scale rendering trees vulnerable to insect attacks have been discussed in the literature. The most obvious ones are the distance to the closest damage, i.e., to breeding grounds for bark beetles, and the biogeographic and climatic conditions, due to geography and large-scale elevation conditions. On the stand scale, the density of conifers in mixed forests, i.e., the density of potential stepping stones, even if not attacked themselves [26], the average age of stands, and general stand health, have all been shown to be critical. Changes to microbial communities around and in roots [27] and mycorrhiza that might have been shared prior to the attack might be also play a role, altering nutrient flow among and to trees [28]. On the site scale, elevation, the exposure of trees in terms of slope and thus wind [15], or solar radiation ([29]; see also below), are important. On the individual scale, factors such as tree age (bark beetles attack rather mature than younger trees; [30]) and tree stress, e.g., from high elevation [31] or drought [32,33], have been found to be key. Also, trees with intermediate crown density have been found to be most affected by bark beetles [30].

The site scale factor of elevation affects tree vulnerability in several additional important ways. First, elevation structures plant communities because of the exposure to the wind and precipitation. At the highest elevations of a catchment, exposure to wind, temperature, and subsequently evaporation [34] can be assumed to be highest, posing the highest stress to vegetation within the catchment [15]. Norway spruce is adapted to such conditions, and low temperatures limit bark beetle developments. However, forests in the surroundings of our site are still in the process of recovering from the major disturbance of acid rain that culminated in the late 1980s [35,36,37]. Before the most recent bark beetle outbreaks, many mountain peaks had already been bare for decades, partly due to the fact that acid rain particularly impacted the highest elevations, and partly as a result of ensuing management. The combination of these factors at the highest elevations in the catchment studied here may have rendered the trees particularly prone to repeated insect infestations.

1.4. Bark Beetles in the Bohemian Forest

The Bohemian Forest, which lies across the Czech–German–Austrian border region and includes the Šumava National Park on the Czech side, has experienced several severe windthrows and bark beetle outbreaks in its recent history. These have affected both actively managed forests and those strictly protected in non-intervention zones [37], but not at the same time and not to the same severity. The most recent large outbreak was triggered by windstorm Kyrill in 2007, affecting ca. six thousand hectares of mature mountain spruce forest [38]. This event initiated heated discussions about the management of wind-felled areas, which may act as a source of bark beetles, putting surrounding forests at risk [1], but also representing important biodiversity hot-spots with positive effects on landscape heterogeneity and recovery [39]. This event also initiated a new wave of interest in bark beetle ecology and population dynamics, which was also fueled by the increasing outbreak risk due to climate change [14,40].

1.5. Hypotheses

We hypothesized that there would be differences in the predictors that would best determine the likelihood of a tree dying during various bark beetle outbreak phases. More precisely, we hypothesized that (1) the importance of factors at different scales (individual site, stand, or regional) would change during the outbreak in the Plešné Lake catchment. We tested whether (2) the site scale was more important than the other scales during the transition from the endemic to epidemic phase in the area, and thus whether the site scale could pre-determine tree survival or dieback. In addition, we tested whether (3) there were important site-scale effects, especially effects from the density of healthy trees in close proximity of an attacked tree. Our expectation was that neighboring healthy trees might potentially protect an attacked tree due to, e.g., shared mycorrhiza before the epidemic phase, while dense stands of healthy trees might attract bark beetles during the epidemic phase. Thus, site-scale effects might play different roles in different outbreak phases.

2. Materials and Methods

2.1. Study Site

The Plešné Lake catchment (0.71 km²) is situated between 13.854 °E and 13.868 °E, and 48.770° N and 48.778° N, at altitudes of 1086 to 1378 m above sea level (Figure 1). The catchment was covered by a primeval conifer forest, strictly protected as a part of the Šumava National Park. The land use classification is 100% forest [41]. About 90% of the catchment was covered by predominantly mature Norway spruce (Picea abies (L.) Karst.) until a bark beetle (Ips typographus L.) outbreak from 2004 to 2008, which led to the death of 75% of the trees [20]. Since then, the forest has increasingly consisted of birch (Betula pubescens Ehrh. and B. pendula L.), rowan (Sorbus aucuparia L.), and European beech (Fagus sylvatica L.), in addition to a large number of spruce seedlings [36]. Apart from the decades-lasting acidic deposition, anthropogenic impacts have been negligible in the Plešné catchment.

The catchment bedrock is formed by granite [41]. The soils are mostly sandy (75%), low in clay content (2%), and include shallow leptosols (38%), podsols (29%), and dystric cambisols (27%), with an average soil depth of 33 cm; the remaining ground surface is bare rock (5%) and wetland (1%) [42]. Annual average air temperature in the Plešné catchment from 2004 to 2017 averaged 3–6 °C [36]. Precipitation averages ~1300 mm annually [43].

2.2. Elevation, Aspect, and Slope Data

The fourth-generation digital elevation model (DEM) was obtained from the Land Survey Office, Prague, Czech Republic (64110-DMR 4G SM 5 (2.5 × 2 km); accessible at www.geoportal.cuzk.cz (accessed on 15 February 2018); metadata webpage https://geoportal.cuzk.cz/(S(2uc3e2sw25prlc3tbf503ezx))/Default.aspx?lng=EN&mode=TextMeta&side=vyskopis&metadataID=CZ-CUZK-DMR4G-V&mapid=8&menu=301 (accessed on 15 February 2018)). The elevations of trees were extracted from the raster map built based on this text file. Slope and aspect were calculated in ArcGIS Desktop 10.8 based on DEM.

2.3. Tree Mapping and Characteristics

Trees were mapped in a geodatabase, similar to previously published approaches. In the first approach, attacked trees were mapped as a layer of points, and areas with wind-thrown trees, mapped by the forest managers, as single layer polygons [18]. The second approach used a LiDAR-derived digital surface model (DSM) to evaluate the canopy height and to differentiate standing from lying trees [38].

Aerial color images covering the Plešné Lake catchment had been captured by manned aircraft and georeferenced by GEOREAL spol. s r.o., Plzeň, Czech Republic. The aerial imagery had a 20 cm spatial resolution. This imagery covered the visible part of the spectrum (380–750 nm) with a measuring camera. Images were photogrammetrically adjusted for the orthophoto map.

In this study, one aerial image per year from the years 2000, 2003, 2005, 2007, 2008, 2009, 2010, 2011, 2013, and 2015 was manually digitized as described by Fluksová et al. [20], covering the full cycle of an outbreak, from endemic to epidemic and back to endemic phase [20]. Each tree’s position was georeferenced in a GIS layer “tree stand”, and the tree state or condition was individually noted. The category of the tree condition was assigned from changes in crown color. Shades of green in the aerial image were classified as “healthy”, unaffected (

Table 1. Classification of trees in the geodatabase.

Class	Explanation
0	future point—there is nothing there at the time of observation, but there will be a sapling/seedling in a later year
1	tall healthy
2	small healthy
3	sapling, seedling
4	tall dead
5	small dead
6	tree stump

). Dead (or dying) trees (affected) were those that appeared as orange, brown, grey, or white. Dead trees were first “standing”, or later “lying”. Where a dead lying tree was degraded, only a “stump” may have remained. Trees that appeared for the first time were categorized as a “sapling” or “seedling”, and attributed a height of 0 m. When the crown diameter of such a tree exceeded two meters, it was moved to the category “young healthy tree” (Table 1). The year in which a tree was first classified as either class 4 or 5 (i.e., as dead) was attributed to the tree as the year of death. Tree digitization and categorization were performed by one person (H.F.) for all images according to the described methodology. While somewhat subjective, the appraisals were consistent throughout the data set.

Tree heights were extracted from LiDAR data for all trees (except for saplings) using the method of local maxima for height estimation [44] (further details below). The terrestrial LiDAR ground surface data (laser point cloud) were taken in 2010 and 2011 by GEOREAL spol. s r.o. In 2010, the data were taken from the lake dam and in 2011, from the trail leading to the ridge. The main steps of the procedure are shown in Figure 2 and are further explained in Supplementary Information, Part B.1.

A visual correction of the calculated DEM was performed as follows (see Supplementary Information, Part B.1 for details). Tree heights were derived from the first return of the LiDAR point as an elevation (height) for the top of a tree. Circular crown projections for each individual surviving tree were manually derived from the local maximum raster. Centroids of the crown projections were converted to a point layer (“tree height” layer) with local maximum height values [44]. Trees from the LiDAR “tree height” layer were then matched to the closest trees in the “tree stand” layer (see Figure 2) and the height from LiDAR was added to each tree point in the “tree stand” layer. This resulted in 11,711 tree heights matched. This approach thus used the transition parameter crown projection to find the corresponding points from the tree layer (i.e., orthophoto digitization) and from the elevation layer (derived from LiDAR).

Height values for all trees that were not assigned a height from the “tree height” layer due to gaps in the LiDAR data were derived from an interpolated pixel value from these newly created rasters (“tree height-interpolated raster”). Buffers (polygons) were drawn around these non-assigned trees so that the polygons represented the crown of the tree. The maximum tree height from the appropriate “tree height-interpolated raster” within a buffer zone around these polygons was extracted and assigned to the respective tree point shape in the “tree height” geodatabase.

Alive and dead standing trees were each categorized as either tall or small based on the tree height from LiDAR, with trees between 7 and 20 m in height categorized as small, and trees from 20 to 50 m as tall. This led to, e.g., a tall healthy tree being assigned to class 1 (Table 1; [20]). Trees were further distinguished into deciduous or coniferous trees, based on crown appearance in the aerial image. For periods with missing aerial images (Supplementary Information, Table S2), annual changes in tree categories were extrapolated (Supplementary Information, Table S3). We classified the tree data beginning with those of 2005 up to and including 2011, so as to represent an epidemic outbreak phase (Supplementary Information, Figure S3). This was based on the combination of the tree counts of the various categories in Table S1, and the cFII index ([45]; Supplementary Information, Part B.2; Figure S3).

The resulting geodatabase (shapefile) was then used for further calculations. The sum of newly dead trees in each aerial image is noted in Tables S1 and S2. We first created subsets of the complete database for time periods up until the respective year named (i.e., data set 2000 comprised trees present in 2000, to represent the period up until and including 2000; data set 2003 comprised trees present in 2001–2003, etc.). In other words, for each such period, we excluded trees of category 0 (Table 1), which would only appear in later years (i.e., future sapling), and we excluded for each year those trees that were already marked dead in previous years (i.e., contained an entry in the “Year of Death” column in the shapefile table which was before or the same as the year of the respective data set), because we wanted to predict the death in a particular year (or period), and not whether a tree had already died previously. This resulted in data sets ranging from 14,214 trees for 2011, when most old trees were already dead but not much regrowth had yet happened, up to 47,838 trees, mainly saplings, in 2015 ([20]; Table S1).

2.4. Geostatistical Analysis—General approach to Spatial Scale Effects

To address the hypothesis that small-scale factors determined a tree’s fate during the outbreak phase when bark beetles were already predominant in the area, we checked effects on different scales. We regarded the “tree height” (see above), longitude and latitude coordinates (abbreviated as “X” and “Y”), “leaf type” (either coniferous or deciduous), and year of “death” as variables on the individual scale (overview in

Table 2. Explanation of the abbreviations of predictors, here for 2003 as an example.

Predictor	Explanation
Individual scale
leaf type	j = coniferous; l = deciduous
Tre_hgh	Tree height (m), estimated from LiDAR
X	X coordinate
Y	Y coordinate
Death	Year of death
Site scale (including up to 30 m)
pl9_s_1	Slope (°)
plec_4g	Elevation (m)
pl3_a_t	Aspect (°)
FII_00, FII_03	Raster cell value from the Forest Infrared Index (UHUL, pers. comm.) for the year 20xx, i.e., FII_00 for the year 2000
c1_30m_j, c2_30m_j, c3_30m_j, c4_30m_j, c5_30m_j, c6_30m_j, c1_30m_l, c2_30m_l, c3_30m_l, c4_30m_l, c5_30m_l, c6_30m_l, c1_10m_j, c2_10m_j, c3_10m_j, c4_10m_j, c5_10m_j, c6_10m_j, c1_10m_l, c2_10m_l, c3_10m_l, c4_10m_l, c5_10m_l, c6_10m_l, c1_5m_j, c2_5m_j, c3_5m_j, c4_5m_j, c5_5m_j, c6_5m_j, c1_5m_l, c2_5m_l, c3_5m_l, c4_5m_l, c5_3m_l, c6_3m_j, c1_3m_j, c2_3m_j, c3_3m_j, c4_3m_j, c5_3m_j, c1_3m_l, c2_3m_l, c3_3m_l, c4_3m_l, c5_3m_l, c6_3m_l, c1_30m, c2_30m, c3_30m, c4_30m, c5_30m, c6_30m, c1_10m, c2_10m, c3_10m, c4_10m, c5_10m, c6_10m, c1_5m, c2_5m, c3_5m, c4_5m, c5_5m, c6_5m, c1_3m, c2_3m, c3_3m, c4_3m, c5_3m, c6_3m	Count of trees of one of six categories (c1, …, c6) within a buffer zone of xx m; if applicable, belonging to tree “species” (l = deciduous; j = coniferous); e.g., “c1_30m_j” stands for the count of trees of category 1, i.e., tall, within 30 m, of “species” coniferous
Stand scale (50–100 m)
d_50_00, d_50_03, h_50_00, h_50_03, etc.	Kernel Smoothing-estimated density of d = dead, or h = healthy trees within a raster cell of 50 m or 100 m, in the year 20xx; e.g., d_50_00 = estimated density of dead trees within a 50 m raster cell in the year 2000
subbasn	sub basin; geostatistically calculated; apart from the four known inflows (see, e.g., Kopáček et al. [36]), two more inflows and their subcatchments were identified and used as a variable
Regional scale (however, some distances turned out to be less than 30 m, i.e., rather site-scale level)
d_ls_xx	Distance to Czech Šumava National Park-listed damages in year 20xx; only available from year 2006 onwards
dm_20xx	Distance to damage in the five years up to 2000, or in the previous years since the last appraisal in the forest state geodatabase; e.g., dm_2000 = damage in the five previous years, i.e., 1996–2000; dm_2003 = damage in the years 2001–2003 according to Senf and Seidl [46,47] which is on a 30 m grain size scale

).

On the site scale, i.e., within 30 m, slope, aspect (in degrees), and elevation were included. The immediate surroundings of a tree encompass the directly neighboring trees, which we addressed by counting the trees from the Fluksová et al. [20] database within a three-, five-, ten-, or thirty-meter radius, to represent site-scale tree densities. We counted trees per category, to address the hypothesis that on the site scale, the density of trees of a certain category would be important, either protecting neighboring trees before the attack, or attracting bark beetles during the epidemic phase. Apart from site-scale characteristics in the Fluksová database, we also used UHUL’s forest health index FII of the respective year as a predictor.

We calculated rasters from the tree densities per area via Kernel Density Estimators (see below) to be able to build models with a % dead trees per area response. If a forest stand is described by its structure as being somewhat homogeneous in terms of—among others—species composition and tree density, then stand size in a topographically complex area such as a catchment with a wide range of slopes and altitudes will range from a few meters to hundreds of meters in diameter (see Figure 1, and compare Moran’s I which is an index quantifying non-randomness in distribution, and thus spatial autocorrelation, in every year high spatial autocorrelation, Supplementary Information, Part A, Table S4; for the method see below). We decided to set the limit between the site and the stand scales at between 30 and 50 m, with 30 m here representing the upper limit of site-scale evaluations, and 50 m here representing the lower limit of stand-scale evaluations. Therefore, we calculated extrapolating rasters from the tree database by Fluksová et al. [20] for 50 and 100 m. A stand may well be larger than 100 m, but the tree database encompassed only a few meters beyond the catchment border, and larger raster cells at the edges of the catchment would have had lower densities because not all counts were available, and would thus have distorted density values.

For the distance to the closest windthrow, i.e., scales potentially larger than 100 m, representing the regional scale, we used two databases. The first consisted of polygons containing various detailed information, but covered only the Czech side (geodatabase “ZjistovaniStavLesa“ of the Šumava National Park administration, pers. comm.; this is an internally compiled long-term study of the Šumava National Park administration to provide information about forest health for all sites within the National Park area in the aftermath of the devastating Kyrill windstorm and massive bark beetle outbreaks, etc.), and the second was a generalized raster data set based on yearly satellite imagery [46,47]. For an overview of the two regional data sources and comparisons with the Fluksová et al. [20] database, refer to Supplementary Materials, Part B.3, Figure S4, and for further details, see below. As can be gathered from Figure S4, some distances from trees to damaged areas recorded in these two databases were in the end less than 30 m, i.e., rather site-scale level than regional.

In order to extract the tree count of a specific type (i.e., small healthy, dead adult, stump) within a given radius around a given tree, we used the command st_intersects() from the R package sf [48,49] for each year. We used four different radii: 3, 5, 10, and 30 m.

2.5. Geostatistical Analysis—Spatial Kernel Density Smoothing to Calculated Stand Density from the Tree Database by Fluksová et al. [20]

Spatial kernel density smoothing was performed using the command kde() in the R package SpatialKDE version 0.8.1 [50], with the setting kernel = “quartic”, to estimate stand density as a predictor for statistical analyses. Non-coniferous trees were excluded for this estimate, since they are not hosts for bark beetles and their state is thus unrelated at the scale at which we used the kernel density smoothing. For kernel density estimations, the band width is decisive. We calculated the appropriate band width using the command bw.smoothppp() in the R package “spatstat.core”, version 2.4-4 [51]. This is a nonparametric test and performs least-squares cross-validation. The calculated band width was 10.24.

We chose cell sizes of 5, 50, and 100 m. The rasters with 5 m cell sizes were only used for plotting (Figure 3) and for calculating Moran’s I (see below), while the 50 and 100 m cell size rasters were used to represent two spatial resolutions of the stand scale. Kernel densities were then calculated separately for healthy trees, combining healthy tall and healthy small trees from the Fluksová et al. [20] database. For further details on the derived geodatabase columns which are predictors as listed in Table 1, see Supplementary Information, Part B.1. For each year, the respective dead, newly dead, or healthy tree density in the kernel smoothed grid cell were extracted for each tree in the Fluksová et al. [20] database, using the command extract() from the “raster” R package [52]. The extracted values were attributed to the respective tree, and included in statistical treatment (see below).

2.6. Forest Damage on the Regional Sale—Czech Data Set of the Šumava National Park Administration

For evaluating the distance to the nearest forest damage from which bark beetles could have come, we used the Šumava National Park administration polygon geodatabase mentioned above (pers. comm.). In this data set, polygons are distinguished into “sites with dry or drying standing trees”, “sites with dry lying trees”, “sites after incidental tree removal”, and “breaks and upheavals not resolved”. The sites after incidental removal are further distinguished into “sites with the removal of most of the wood”, and “sites leaving most or all of the wood mass to decay” (further explanations in Supplementary Information part C). This geodatabase was only available from 2006 to 2021, thus, information for most years prior to the Kyrill storm (2007) was not available. We assumed that the “sites with dry or drying standing trees”, “sites with dry lying trees”, “sites leaving most or all of the wood mass to decay”, and “breaks and upheavals not resolved” were the most important as potential bark beetle breeding sites. Bark beetles usually only breed in the bark, so lying dead wood from which the bark was removed may not offer breeding grounds. “Not resolved” implies that the bark was left on the trees, and such dead wood might thus become a breeding ground.

For each tree of the Fluksová et al. [20] database, we calculated the smallest distance to the nearest polygon of damage in the respective year, using the command st_distance() from the R package sf [48,49], assuming that each polygon included not only fresh disturbances, but also legacy disturbances for the respective year. These data are only available from 2006 onwards.

2.7. Forest Damage on the Regional Sale—European Scale

Since the detailed data set from the Šumava National Park only covered the Czech area, and was only available from 2006 onwards, we also used the European scale database by Senf and Seidl [46,47]. For more explanation, refer to Supplementary Information, Part B.3.

In the final database, the smallest distance of each tree to the nearest damage central point for the respective year was calculated using distanceFromPoints() from the raster R package [52], resulting in the shapefile column “dm_year”, where “year” was the respective year for which the database provided an entry and which overlapped with the time period in the tree database.

Both for the 50 and 100 m raster cell kernel smoothed rasters that had been based on the Fluksová et al. [20] database, and the rasters from the Senf and Seidl [46,47] data set which are on a 30 m grain, we also used rasters from the respective previous years or periods as predictors, because a disturbance might not act immediately. In order to keep the set of predictors to a minimum, and since the site scale should operate at shorter time periods, we used only the tree class counts (for the definition of the six tree state classes, refer to Table 1) within 3, 5, 10, and 30 m from the same year for the site-scale predictors.

We compared the distances to the closest disturbance from the Šumava polygon geodatabase with those from the database by Senf and Seidl [46,47] in Figure S5. The distances from the Šumava geodatabase were consistently smaller than those from Senf and Seidl [46,47], almost certainly due to the smaller size of some of the polygons in the Czech database. Both databases, however, consistently show that the distances to disturbances were the smallest in 2006–2008 (Figure S5).

2.8. Statistical Treatment of the Results from Geostatistical Analyses

2.8.1. Spatial Distribution

To test the hypothesis on the random distribution of dead and living trees, we tested the autocorrelation of the year in which a tree died (“Year of death”, or “death” in short in the following tables and figures) with Moran’s I.

2.8.2. Gradient Boosting Machine (GBM) Predicting Tree Death

We used the machine learning tool Gradient Boosting Machines (GBMs; [53]) for classifying the yearly outcome of an individual tree as either alive or dead, based on the variables that were named or derived as above. Like random forests, GBMs build ensembles of trees, but have been found to be superior in some applications [54]. While random forests build deep independent trees, GBMs build shallow trees and incrementally minimize the error term [55]. The trees are built sequentially, and each tree learns and improves on the previous one. Shallow trees by themselves are weak predictive models. However, in GBMs, they are “boosted” to produce a powerful “committee” [56]. GBMs, when appropriately tuned, are thus often hard to beat with other algorithms [57]. They deal with multicollinearity in the data [58], and when there are collinear variables, the tree takes one of them and the accuracy of the model stays the same no matter which of the collinear variables are used. In the variable importance calculations, they may achieve the same or similar importance, and thus, it is important to not just focus on the most important predictor. This was demonstrated by using GAMs in the next step with all variables that had had a variable importance of at least four in any year (see Section 2.8.3 below). For calculating the GBMs, we used the command gbm() in the R package “gbm” v. 2.1.8.1 [59], setting the binomial distribution “huberized”, and narrowed down the respective settings with a grid search, as recommended in [60].

On this fine spatial scale, we focused on the interaction between coniferous and deciduous trees and therefore included not only counts of coniferous trees, but also counts of deciduous trees and total counts of both deciduous and coniferous trees as predictors. All trees, including deciduous, are part of the ecosystem and stabilize the soil. A mixture of trees might prevent bark beetle attacks, which focus on dense stands of conifers [61]. To cross-check the general patterns of the GBM results, we also calculated random forest models, using the command “randomForest” of the R package “randomForest” v. 4.7-1.1, [62].

Following the GBM yearly analyses, we plotted the three most important variables that predicted tree deaths in variable importance plots as well as line plots to visualize the importance of the variance of all variables over time and to check whether the variable importance changed throughout the epidemic phases. Because the total number of predictors (111, see Table S5) is too high to visualize, we concentrated on those predictors that achieved a variable importance threshold of at least 4% in at least one year. Some collinearity between variables of similar character, i.e., between c1_5m and c1_3m, was expected. GBMs deal well with collinearity, as explained above, but it is possible that of two collinear variables, one is attributed a high variable importance in one year, and the other in another year. Plotting all variables that passed our 4% threshold at least once allowed us to check whether such similar variables became interchanged among years. In such cases, both had to be regarded as important and the specific value of variable importance was less decisive. We thus used the individual-scale GBMs to narrow down the list of variables to be used for general additive models (GAMs) in the subsequent step.

2.8.3. General Additive Models (GAMs) Predicting the Percent of Conifers Killed

We chose the nonparametric spline fitting GAM method [63] for fitting models to show which variables predicted the percent kill rate (number of newly dead trees per sum of newly dead and healthy trees × 100) in the kernel-smoothing-derived rasters. While the GBMs were calculated on individual trees and models ran separately for each year, the GAMs were calculated on raster data, like in previous studies, and the response was the % of trees killed, not the individual tree’s survival. On this coarser spatial scale, we focused on coniferous trees only, since they are the ones succumbing to the bark beetle attack. As described in the previous section, the variables used here as predictors were those that showed a variable importance of at least 4% in at least one year. We used the command gam() from the R package “mgcv” [64], with the formula response ~ s(predictor) and the settings bs = “cs”, method = “REML”.

3. Results

3.1. Spatial-Temporal Dynamics

In 2000, just a few dead trees were scattered across the catchment (Figure 3). By 2003, a few hot-spots of tree death appeared in the northwest, and in the following period, 2004–2005, more trees died in a strip along the northwestern edge of the catchment. From 2006 to 2007, tree death dominated in the southeast of the catchment, while from 2008 to 2009, trees died in a few of the areas that had been spared earlier. Only a few additional trees died between 2010 and 2015 (Figure 3).

3.2. Spatial Autocorrelation: Moran’s I

Although Figure 3 shows clear patterns in the temporal progress of tree death, there was significant spatial autocorrelation in every observed year, regardless of the phase of bark beetle attack (Supplementary Information, Table S3). Figure 4 shows how the spatial autocorrelation, estimated using the observed Moran’s I, coincided with new tree deaths, particularly in the first phase. After the climax of tree death, new mortality largely declined, but the spatial autocorrelation between dead trees remained high.

The spatial correlation of the year of a tree’s death according to Moran’s I was highly significant (I between 0.002 and 0.1; p < 0.0001) for all years, contradicting our first hypothesis of neutral dispersal. In more detail, the distribution of dead tall and dead small trees per year was also non-neutral (Supplementary Information, Table S4). This can be seen in Figure 3, with patches of early death in the south of the catchment. In contrast, in the north, where tree death was earliest and most extensive, a few islands of healthy spruces remained even until 2015 (Supplementary Information, Part B.4, Figure S6). The earliest tree death in the north might have resulted from windthrows, as suggested by Figure S7. It is, however, unclear where the bark beetle invasion came from. The windthrow damage within the Plešné catchment was largely in line with damage throughout the Šumava area (Figure S8). Only beginning with 2010 were there any spatial patterns discernible, with the north and south border of the catchment being closest to the nearest damage (Figure S9).

We also checked Moran’s I on the rasters that were derived from the individual tree positions. The spatial resolutions are thus consecutively coarser, from 5 m cell size, over 50 to 100 m cell size. Moran’s I was highest for the 5 m cell size, with maximum values of ~0.8 (Figure 5). The values for the 50 and 100 m resolutions were similar to each other and in the range of Moran’s I for individual tree positions (compare with Figure 4). There were no clear trends with time, particularly for the 100 m resolution.

3.3. The Most Influential Variables Predicting Tree Death

The variables that had at least a 4% influence in at least one year of the observation are shown in Figure 6 and in Supplementary Information Part D. In brief, the major predictors were as follows: tall and small healthy conifers; tall and small dead conifers within distances of between 3 and 30 m; altitude; tree height; distance from the closest damage of the Senf and Seidl [46,47] raster data set for 2003; and coordinates of latitude and longitude. Full details are given in Table S5.

Surprisingly, tree height was always the best predictor, even if its influence varied over the years. It was lowest in 2010, when the variables representing dead tall conifers within 3 m and healthy small conifers within 30 m were almost as important (Figure 6 and Figure S16; note, however, that the predicted tree death was negative, i.e., the model was not reliable). The individual property representing latitude was the third best predictor in 2005 (Table S5). The other most influential predictors were mainly site-based, and were usually the second or third most important predictors (Figure 6 and Table S5). No stand-scale predictor was important. Among the regional-scale predictors, only the distance from the disturbance raster of 2003 according to Senf and Seidl [46,47] was somewhat important, and only in 2008 (Figure 6 and Figure S14), but the actual estimates of a tree to be dead or alive were negative, thus again unreliable. None of the variable influences stayed the same during the endemic phase (2000–2003) or during the epidemic phase (2005–2011), and there were no clear patterns over time.

The regional scale only became important as a predictor at the peak of the epidemic phase, likely because the Czech forest state database did not contain information for the years before 2008. Note, however, that beginning with 2011, the estimated values were calculated to be below 0, and are thus not plotted in the figure—in other words, they were artifacts, and cannot be interpreted. Already in the data for 2007, there were gaps in the range for which predictions could be calculated, and short lines, not spanning the whole range of predictors, are the result. The year 2007 was also when considerable windfall was observed after storm Kyrill [38]. Thus, from 2007 on, dead trees not only reflected the bark beetle attack, but also the windfall.

To exclude the possibility that the GBM method we used influenced the results, we also conducted random forest analyses from the same data sets as were submitted to the GBMs. Figures S19–S27 show that the general patterns were the same: tree height emerged as the most important variable in all years, followed by site-scale variables of tree densities within 30 m distances.

For the variables that passed our 4% influence threshold we calculated GAM models on the 50 m and 100 m cell size rasters to test which of these models were statistically significant. In Figure 7, as an example, the results of the GAM models are given for the predictors shown to be most important by the GBM of 2000 (all other GAMs are in the Supporting Information). Not all predictors produced valid GAM models. An intermediate tree height of just under 20 m was correlated with the highest proportion of dead trees in 2000 (Figure 7). Both very low and very high numbers of tall healthy coniferous within 30 m were linked to the lowest percentage of dead trees, while higher numbers of small healthy coniferous within 10 and 30 m were associated with slightly lower percentages of dead trees. Altitude and northing showed nonlinear trends (Figure 7).

Likewise, in 2003, the proportion of newly dead trees within a 50 m raster cell increased with increasing tree height (Supplementary Information, Figure S28). In other years, patterns were less clear (Supplementary Information, Figures S29–S35). The percentage of dead trees sharply decreased in 2003 with an increasing density of tall healthy coniferous within 30 m (Supplementary Information, Figure S28). The pattern was similar in 2005, but in 2007 and 2008, i.e., during the epidemic phase, the percentage of dead trees sharply increased with an increasing density of tall healthy coniferous within 30 m. The behavior of all other variables also depended on the year of observation and the patterns were largely nonlinear. The same was true for GAM models of the 100 m rasters (Supplementary Information, Figures S36–S44).

4. Discussion

The observation that individual- and small-scale tree properties are the most important factors determining a tree’s fate before, during, and after a bark beetle epidemic is only valid within a catchment. The outbreak/epidemic phase in 2005–2009 in the Plešné catchment itself was most likely triggered by windthrows in the vicinity of the catchment, i.e., on a regional scale. The effective dispersal of bark beetles is several hundred meters [7,15]. There is always some external trigger needed to start the outbreak, and it cannot start without some additional driver (e.g., a windthrow or severe drought). Individual- or site-scale factors shown to be influential within the Plešné catchment cannot by themselves trigger or prevent a bark beetle outbreak and epidemic. However, which tree survives within such an epidemic is decided on the individual or site scale. This may partially explain why different catchments react very differently to the presence of bark beetles.

In this study, individual-scale and site-scale predictors performed better at predicting tree death on the individual scale (with GBM) and on the stand scale (with GAM) than stand-scale or regional-scale predictors. In particular, tree height stood out as an important predictor.

However, the patterns with which tree death could be predicted from predictors at different scales were largely non-linear in the GAM models. This is particularly clear for one of the best predictors, latitude—if latitude predicted tree death linearly, tree death would be most pronounced in either the north or the south. Instead, tree death usually occurred in clusters or nests, some of which were at some point at more northern points than others. Latitude’s prediction value thus occurred with several minima and maxima along this axis, but these “islands of death” did not show any trend in number or density. The reason that some predictors predicted the response in a more linear fashion in some years but not in others may lie in a particular variable’s interactions with other predictors. Depending on the range of an interaction partner in a respective year, their interaction may be stronger or weaker.

The catchment studied here is small and a unilaterally oriented glacial karst, so we assumed that variables that proved most influential in other studies, e.g., mean growing season temperature [17,23] and mean growing season precipitation sum [65], would be less influential in this study. However, these climatic variables were not measured at enough points in the catchment for us to be able to extrapolate them to every tree in the whole catchment. Also, because of the unilateral orientation, they probably did not vary much in this small catchment, and would likely have correlated with the tested variables altitude, slope, and aspect. Soil moisture has been shown to be influential in a larger scale study [66], but like climatic variables, in this study, soil moisture was measured at too few points to extrapolate to every tree in the catchment.

4.1. The Importance of Predictors at Different Scales during the Bark Beetle Attack Phases

In contrast to our hypothesis, the importance of predictors at different scales (individual scale versus site scale, stand scale, or regional scale) changed only slightly from one phase to the other. In the endemic phase in 2000, the highest percentage of dead trees within a 50 m raster cell was linked to an intermediate tree height. However, just before the main phase of tree death, in 2003, the highest trees were found in those raster cells with the highest proportion of tree death (Figure S28). Thus, it seems that the tallest trees were already attacked in 2003, and might actually have already been in the epidemic phase.

During epidemic phases, bark beetles attack healthy trees when the trees are clustered together [8]. This can be seen in 2007. Before 2007, the relationship between tall healthy coniferous trees within 30 m (c2_30m_j) and the proportion of trees dying (c4_xm, with x standing for distances between 5 and 30 m) had been largely negative (e.g., Figure 7). From 2008 on, the relationship between tall healthy coniferous trees within 30 m (c2_30m_j) and the % of dead trees was positive (e.g., Figure S31), and thus a high number of healthy conifers surrounding a coniferous tree did not prevent high numbers of trees from dying. This increase in tree death with an increasing number of healthy trees ceased in 2010, but confidence intervals for that year in most models were wide, i.e., patterns were not strong. The reason probably lies in the small number of healthy trees that had survived until 2010.

Ips typographus, the species responsible for all infestations of Norway spruce stands in the Šumava National Park, has been shown to be deterred at 1.5 m from a source of pheromones and to prefer empty spots without pheromones within this small distance [67]. This is a density-regulating mechanism, since each tree can only support a limited number of beetles [67]. However, at 3 m from the source, only a few beetles preferred empty spots [67] and instead seemed to be attracted by the pheromone [68], but this preference rapidly decreased beyond a distance of 9 m [69]. In another study, within the first year of an outbreak, only 5% of fallen trees were affected, while the proportion was 50% in the subsequent year [7]. The infestation thus was negligible in the first year after the disturbance, but started to be significant in the second year. While our study is not experimental but observational, we confirmed that distances less than 10 m between a dead tree, i.e., a potential bark beetle source, and healthy trees may be important, particularly in the pre-epidemic phase. We found two variables within 3 m and one within 5 m to be among the thirteen most important variables, while stand-scale predictors were never significant, and regional-scale predictors were only significant in 2008. Another study at an adjacent site (and therefore with similar conditions) [70] used very high-resolution satellite and drone imagery to develop a blueprint that could be useful in predicting potential dispersal patterns that the bark beetles might follow. This blueprint should assist management in separating the unhealthy (attacked) trees from the healthy ones in the early epidemic phase, aiding in the protection of healthy trees [70]. Such an approach would likely be more efficient than a detailed manual characterization of individual trees as was done at present.

4.2. Aerial Images Versus Satellite Remote Sensing

In this study, we used a combination of manual appraisals of aerial imagery to derive individual tree properties, and satellite (LiDAR) data, to derive individual tree height. While the manual interpretation of the aerial imagery was particularly time-consuming, we think that this combination was necessary for clarifying which factors most influence the susceptibility to bark beetle attacks at the individual tree level. For other studies, it might be advantageous to combine LiDAR with very high-resolution remote sensing.

4.3. Is the Site Scale the Most Influential Scale?

We hypothesized that, during the epidemic phase, site-scale factors would be most important. We assumed that bark beetles would already be so widespread and abundant in the area that individual features would no longer play a role. However, the GBM variable importance models showed that, in each year, and in each phase, the individual-scale variable tree height was the most important variable in predicting the death of a tree. Another individual-scale variable, latitude, was among the most important predictors, but did not play a discernible role. However, the second and third most important variables were indeed site-scale variables, namely, the densities of healthy or dead trees within varying distances. Interestingly, aspect did not play a role, in contrast to, e.g., a study in the Dinaric mountain forests of Slovenia [71], even though the Plešné catchment has steep slopes and a mosaic of aspects which would be expected to influence tree susceptibility and resistance.

4.4. Does the Density of Healthy Trees on the Small Scale Protect Trees?

Among the influential variables in the GBM model on the site scale were tall healthy coniferous trees within 3 and 30 m, respectively, and small healthy coniferous trees within 5, 10, and 30 m. The hypothesized negative relationship with the percent of dead trees within a 50 or 100 m raster cell in the GAM model was, however, present only in a few years, e.g., on the 50 m scale response. Such negative correlation between the increasing density of healthy trees and the percent of dead trees within a raster cell might indicate, e.g., a diverse root microbial community protecting trees until an epidemic phase. For instance, the roots of healthy trees are surrounded by a diverse network of mycorrhizal fungi that cycle nutrients to and from the tree roots [72,73]. This should be all the more effective within a short range between living trees, i.e., when the density of living trees is high. After bark beetle attacks, symbiotic root fungi and most enzyme activity has been shown to decrease and endophytes increase [27], which may negatively affect spruce [74,75]. The disturbance of trees’ phloem flow to the root, transporting photosynthates, seems to be responsible for the disappearance of ectomycorrhizal fungi [76], disrupting the nutrient transport to trees [77].

Another factor that might have been responsible for the small-scale variables and individual-tree characteristics being more important than larger scale variables, is that trees in the Šumava National Park have already been stressed by increasing temperatures (e.g., [43]). Previous studies have shown that tree health estimated from satellite images (Landsat) was already decreased prior to the bark beetle attacks (Figure S3 and [45,78]), possibly making this area more susceptible to bark beetle attacks than areas with healthier trees.

The fact that we found patterns at the 100 m scale to be less important than at the 50 m scale, in contrast to other studies that have shown stand scales as large as 180 m to be decisive for stand health [23], is probably due to the spatial extent of the data. In our data set, 100 m is probably too large to represent a homogenous stand, because the buffer zone around the catchment, within which trees were counted by Fluksová et al. [20], was less than 100 m wide, and so the densities at the catchment margins may have been underestimated. Therefore, our results including variables and rasters at a 100 m spatial resolution must be evaluated cautiously.

4.5. Are Spatial Scales Useful for Discussing an Ecological Effect?

From an ecological point of view, the differentiation of spatial scales may seem artificial and thus unhelpful. A tree is attacked under two conditions: (1) there is a population of beetles and (2) a tree is susceptible and not able to defend itself. The presence of a beetle population depends on the number of parent beetles and their dispersal, which is a function of climate. The susceptibility of a tree is a function of climate, tree age and size, and local stem density. Thus, climate is an important factor which also varies locally. However, within the small Plešné catchment, climatic conditions are similar, except for air temperature that decreases with increasing elevation. Hence, we used longitude and latitude as a proxy for small-scale effects, and indeed, based on the statistical analyses, they played a minor role in determining a tree’s fate. We had expected slope and aspect to have more of an influence on a tree’s death or survival, since they not only vary more than climate in the catchment, but also predetermine incoming solar radiation to the forest and affect air and soil temperatures. However, we found coordinates to have stronger effects, probably because bark beetles mostly spread from the southwest to the rest of catchment.

Density and tree age are not always homogenous across stands, so we considered density and tree age not as a stand-scale predictor, but as a site-scale or individual-scale predictor, respectively. However, most drivers act on more than one scale. We found tree properties and their location (specifically latitude) to be more important in determining an individual tree’s fate in this small catchment than the distance to the next patch or raster cell of damage. Generalizations are difficult, but for forest management, it is important to know at which spatial scale an intervention has which effects.

4.6. Salvage Logging for Managing Bark Beetle Attacks

Salvage and sanitation logging consist of taking dead trees (a potential infestation nucleus) out of an area and are usually performed at a scale of hectares to remove stepping stones for pest dispersal. Our study, however, shows that on a catchment scale, the distance to the closest disturbance, i.e., to the closest stepping stone, is of less importance than features at the individual- and site-scale, i.e., within 30 m. This has important implications for forest management. In an area where trees are already stressed by high temperatures and where bark beetles are already present, salvage logging may not be an effective measure. Instead, management might focus on improving the condition of trees, particularly the tallest ones, to prevent them from becoming attacked. For instance, the thinning of trees within planted stands might improve individual fitness by reducing competition, which gives individual trees an advantage in withstanding infestation [79]. In contrast, salvage logging has been shown to increase the susceptibility to windthrows at newly created stand borders [80]. In addition, microclimatic conditions at forest margins, e.g., higher bark temperatures, and olfactory signals from standing edge trees have been shown to favor beetle swarming [81]. Thus, new forest margins, as a result of salvage logging, may actually facilitate bark beetle infestation into previously unaffected stands.

5. Conclusions

The spatial patterns of tree death following a bark beetle outbreak in the Plešné catchment showed clear dependences on individual tree properties and on tree densities at the site scale (within 30 m). Stand-scale features (50–100 m) did not predict tree death, and the regional scale (>100 m) was of minor importance. Implications for forest management in similar catchments include that individual tree survival during a bark beetle infestation can be achieved by avoiding hectare-scale salvage logging and by taking care that as many trees as possible, including deciduous trees and saplings, survive.