Nonparametric upscaling of bark beetle infestations and management from plot to landscape level by combining individual-based with Markov chain models

Linked to climate change, drivers such as increased temperatures and decreased water availability affect forest health in complex ways by simultaneously weakening tree vitality and promoting insect pest activity. One major beneficiary of climate-induced changes is the European spruce bark beetle (Ips typographus). To improve the mechanistic understanding of climate change impacts on long-term beetle infestation risks, individual-based simulation models (IBM) such as the bark beetle dispersion model IPS-SPREADS have been proven as effective tools. However, the computational costs of IBMs limit their spatial scale of application. While these tools are best suitable to simulate bark beetle dynamics on the plot level, upscaling the process to larger areas is challenging. The larger spatial scale is, nevertheless, often required to support the selection of adequate management intervention. Here, we introduce a novel two-step approach to address this challenge: (1) we use the IPS-SPREADS model to simulate the bark beetle dispersal at a local scale by dividing the research area into 250 × 250 m grid cells; and (2) we then apply a metamodel framework to upscale the results to the landscape level. The metamodel is based on Markov chains derived from the infestation probabilities of IPS-SPREADS results and extended by considering neighbor interaction and spruce dieback of each focal cell. We validated the metamodel by comparing its predictions with infestations observed in 2017 and 2018 in the Saxon Switzerland national park, Germany, and tested sanitation felling as a measure to prevent potential further outbreaks in the region. Validation showed an improvement in predictions by introducing the model extension of beetle spreading from one cell to another. The metamodel forecasts indicated an increase in the risk of infestation for adjacent forest areas. In case of a beetle mass outbreak, sanitation felling intensities of 80 percent and above seem to mitigate further outbreak progression.


Introduction
European forests are under accelerated effects of climate change, showing decreasing intervals between extreme events, as well as increased drought risks (Lindner et al. 2010). Interactions among such change drivers might potentialize negative impacts (Seidl et al. 2017). Droughtinduced mortality (Allen et al. 2010), for example, may not only lead to changes in tree growth due to altered intra-and inter-specific competition (Piao et al. 2011;Sánchez-Salguero et al. 2012) but also to changes in species distribution (Delzon et al. 2013). Rainfall deficits in combination with extreme summer heat, as observed in 2003, lead to a decrease in primary production and could convert temperate ecosystems into carbon sources (Ciais et al. 2005). The effects of extreme events and increased drought risk on current and future biotic and abiotic disturbances are also not expected to subside with time (Allen et al. 2010;Dale et al. 2001;Seidl et al. 2020;Senf et al. 2019;Thiele et al. 2017). Warmer and drier conditions, as well as increased disturbances such as windthrow also favor pests such as the European spruce bark beetle Ips typographus L. (Seidl et al. 2016(Seidl et al. , 2017Temperli et al. 2013). In recent years, beetle mass outbreaks have been observed in both managed and Communicated by Claus Bässler. unmanaged forests. Although outbreaks in large protected areas without management interventions, such as national parks, might provide valuable insights into its starting conditions and development (Lausch et al. 2011), these sites can provide an infestation risk for managed surrounding forests, requiring an adaptive management to avoid negative impacts on non-protected areas (Mezei et al. 2017;Zolubas and Dagilius 2012). On the other hand, while felling infested trees is inconsistent with most conservation goals of national parks (Kulakowski 2016), it is suggested as one of the most effective measures against Ips typographus outbreaks (Jönsson et al. 2012;Kulakowski 2016).
We apply the novel individual-based model IPS-SPREADS (Pietzsch et al. 2021) to investigate the impacts of climate change and forest protection management on the long-term development of bark beetle outbreaks and spreading risk from the Saxon Switzerland national park, Germany, to adjacent managed forests. A mass outbreak of Ips typographus was observed from 2015 to 2018 in this area ( Fig. 1). IPS-SPREADS was developed from the IPS model (Kautz et al. 2014 by implementing GIS import capabilities to promote its application to real forests. This model describes beetle spreading mechanistically, considering explicitly the energy of each beetle, the energy costs during their flights, and the attractivity of trees among others, which is an advantage over statistical models. Like most individualbased models (IBMs), emerging model results and patterns arise from the simulated individuals' traits (Railsback and Grimm 2011). Therefore, the forecasted dispersal and outbreak patterns in the IPS-SPREADS model depend on the individual traits of both beetles and trees. The model can be applied, for instance, to test the effects of sanitation felling on the outbreak, beetle population fitness and the number of damaged trees. As most IBMs, however, the IPS-SPREADS model is computationally expensive and applicable only for smaller forest plots. Upscaling from plot to landscape level, as well as forecasting spatio-temporal dynamics of bark beetle outbreaks 20 years into the future requires a metamodel (MM) approach, based on a general nonparametric framework (Cipriotti et al. 2016).
Such MM describes the behavior of an underlying model on a higher hierarchical level (Gore et al. 2017;Moorcroft et al. 2001;Urban 2005), providing an efficient way to facilitate the predictions of the behavior of the individual-based model (IBM, Fig. 2) over a wide range of parameter combinations and scenarios (Pietzsch et al. 2020). The use of MMs reduces simulation costs and promotes communication and understanding of simulation models' behavior (Kleijnen and Sargent 2000;Mertens et al. 2017). Regression, Bayesian emulators, machine learning, differential equation and Markov chain (Pietzsch et al. 2020) are the most common types of MMs used in an IBM context, being used as prediction tools, as well as are applied for sensitivity analysis, calibration and upscaling (Pietzsch et al. 2020). We used the Markov chain MM developed by Cipriotti et al. (2016) as a template for our MM, as this model type has proven effective for IBM upscaling (Pietzsch et al. 2020). This applies in particular when considering land cover transitions and neighbor interactions, which are of importance for the development and spread of a bark beetle mass outbreak (e.g., Seidl et al. 2016;Wildemeersch et al. 2019).
We hypothesize that (1) increases in number of bark beetle generations per year due to forecasted increases in mean temperature leads to faster outbreaks and more frequent complete dieback of Norway spruces (Picea abies) inside and outside the national park. We also propose that (2) sanitation felling (i.e., removal or debarking of infested trees) at plot level inside the national park can lower the risk of infestation spread from the protected to adjacent forested area in the next 20 years, but requiring a minimum threshold value to be effective.

Data processing
First, a representative selection of grid cells from the research area was chosen using a shapefile of forest stand polygons (Saxon Forest State enterprise; Staatsbetrieb Sachsenforst). We then selected 7524 grid cells of 250 × 250 m cells to model beetle outbreak within the approximately 14,100 ha of forested area. The relatively small size of the chosen area was representative of heterogeneity of the landscape, while limiting the maximum number of beetles to be simulated during possible mass outbreaks to a tractable number. Of the selected area, 2257 cells had spruce trees and 170 of them exhibited an initial infestation in 2015, being used to simulate the outbreak development in the following year using the IBM IPS-SPREADS. To choose the grids used for beetle outbreak simulation, a classification of the three main features of all grid cells was performed and allowed the selection of a reasonable number of representative areas. Therefore, we calculated the mean kairomoneinduced primary attractiveness of trees to beetles, and both the proportion of spruces and the number of trees already infested by beetles in each site. The frequency distributions of these three features (Fig. 3) were used as classification criteria, so that the number of grid cells within each class was approximately the same (ranging from 53 to 61 cells for each class). This ensured we selected at least one grid cell The classification of features resulted in 36 possible combinations (3 classes for primary attractiveness, 3 classes for spruce proportion and 4 classes for infestation level), with at least one grid cell per combination (see supplementary material, Table S1). For each grid cell of the research area, the mean deviation of its features was calculated and used to identify which combination of grid cells was closest to the respective classification level mid. This was done to ensure that effects associated with each feature were as close as possible to a given classification level.

Individual-based simulations
We followed the ODD protocol (Grimm et al. 2020) to describe the applied individual-based model IPS-SPREADS. The IPS-SPREADS model has been described by Pietzsch et al. (2021) and was used here with small modifications. The random placement of spruce trees was done beforehand during the data processing with R, instead of during the setup of the IPS-SPREADS model and incorporated into the imported raster file PAS.asc to ensure that the amount and placement of spruces was equal for all model runs. We also introduced a new global result in the form of number of trees cut after each generation: ncut1, ncut2 and ncut3. The full ODD protocol with adjustments to the original from Pietzsch et al. (2021) can be found in GitHub (linked in the data availability statement) section.
The overall purpose of using IPS-SPREADS was to understand and investigate the impact of individual traits of beetles and trees as well as measures of bark beetle management (e.g., sanitation felling) on the spatial pattern of outbreaks of the spruce bark beetle. We used patterns of location, number and progression from 2015, 2016, and 2017 of infested spruce trees (Picea abies) within the Saxon Switzerland national park to model beetle outbreaks from real data.
The model includes patches, bark beetles and volatiles as entities. The global environment (a forested area) is characterized by wind speed and direction. Patches represent 5 × 5 m areas containing one spruce or are empty. Characteristics of these patches include their primary and secondary attractiveness for bark beetles, a capacity for infesting and emerging beetles and, one of three levels of infestation: not infested, infested or, fully occupied. Patches are squares in the horizontal plane with data about tree height, tree species and capacity for infesting bark beetles. Beetles have unique values of energy level and energy efficiency. Volatiles are emitted by beetles infesting trees and carry information about the location and the total attractiveness of their source at the time of emission. The spatial extent is constrained to the imported GIS data. If a mass outbreak should be simulated, a maximum extent of a quarter square kilometer is advised. The temporal resolution of IPS-SPREADS is defined by time steps, where 200 steps are contained in one day. The overall length of the simulation depends on the chosen parameter settings (e.g., if one, two, or three generations are to be simulated). The simulation ends as soon as all beetles emerged from the beetle source trees and no more dispersal flight or tree attacks are taking place. level ranged from 53 to 61. Primary attractiveness relates to the kairomone-induced attractiveness of each tree to bark beetles (Kautz et al. 2014 The most important processes of the model, which are repeated every time step, are the (1) emission of volatiles by infested trees, (2) beetle flight, (3) beetle's decision to attack a patch, and (4) defense of attacked trees. Each simulation begins with the model setup during which GIS data like tree height or the begin of beetle flight is imported and used to create the model world by distributing spruce trees and assigning trees as beetle sources. After the artificial environment is set up, the first day of the simulation starts with the generation of beetles on all beetle source trees and the updating of all plots in the model interface. Then, on each time step, volatiles are emitted by infested trees and the beetles move either to the source tree of the most attractive volatile within the beetle's perception range, or to a neighboring patch using a correlated random walk. After all beetles have moved, it is checked if the beetle's energy reserve is higher than the total attractiveness of the beetle's current patch. If the target is an already infested tree, the beetle is added to the infesting beetles on this patch. Otherwise, the beetle will wait on the given tree for one day. If enough beetles attack the same tree within one day, the tree defense is overcome and the tree becomes infested. Otherwise, all beetles waiting on the patch are killed by the tree defense. At the end of each day, it is checked if all beetles of the given beetle generation were generated and dispersed. If not all beetles were generated and dispersed, the next day starts with the generation of beetles and the update of all plots. Otherwise, bark beetle management takes place and a predefined proportion (i.e., 50%) of infested trees is removed from the model world by sanitation felling. This is done before brood emergence or sister brood establishment can take place to simulate an optimal timing of management action. If more than one beetle generation was defined, all infested trees that were not removed by sanitation felling are shifted to beetle source trees for the next beetle generation. After this, a new day with the same schedule of actions as described above begins. When the predefined number of generations has been simulated and the sanitation felling for each of those generations has taken place, model results are calculated, the pattern of infested trees is exported as a raster layer (.asc file) and the simulation stops.

Simulation experiments
After processing, the grid cells chosen for simulating the 36 combinations of mean primary attractiveness, spruce proportion and number of infested trees (Table S1) were imported into IPS-SPREADS. During the simulation experiments parameters were set to the expected values, while both wind speed and duration of beetle flight activity were modified: the wind speed was set to zero to avoid beetle drifting out of the relatively small grid cells, and the duration of beetle flight activity was set to 40. The latter was done, as the online version of the phenology model PHENIPS (Baier et al. 2007) recorded for the weather station within the research area (Lichtenhain-Mittelndorf) 80-100 swarming days in years with two beetle generations and 120 possible days in 2018 where three full generations took place. For the sake of simplicity, we therefore used 40 days of swarming for each simulated beetle generation in the IPS-SPREADS model.
The number of infested trees, the number of living spruce trees, the number of trees cut during sanitation felling if applied and the mean primary attractiveness of all living spruce trees were measured at the end of each simulation experiment. The number of infested trees was used as number of beetle source trees for the following year and to calculate the proportion of dead trees on each site. This proportion was used to assess the impact of management and climate-induced increase in number of beetle generations on the development of the long-term outbreak.
To account for a possible impact of climate changes on the bark beetle population and outbreaks, we varied the number of beetle generations per year from two to three, where more generations per year suggested a stronger increase in mean annual temperature (hereby climate change scenarios). Using such approach ensured that climate-induced increases in mean temperatures on beetle development, flight activity duration and tree susceptibility were considered. In addition, forest management scenarios were investigated. We varied the intensity of sanitation felling (removal or debarking of infested trees) from zero to 100 percent in 10% steps after each beetle generation. As both sanitation felling and climate change scenarios per year were applied simultaneously, a total of 22 combinations were simulated on the selected grid cells (n = 27). Performing ten replicates each, the total number of runs with IPS-SPREADS added up to 5940. Only 27 of the selected 36 grid cells were simulated, as nine cells did not exhibit any prior infestation prohibiting beetle emergence and spread in the IPS-SPREADS model.
The results of the simulation experiments with IPS-SPREADS were stored as a data table for each combination of climate change and management scenario. These data tables contained the mean primary attractiveness, spruce proportion and number of infested trees at the end of the simulation with IPS-SPREADS separated by simulated grid cell. If at least one level of the three aforementioned features changed during the simulations, a new combination label was assigned based on the new combination of features and the previously applied labeling of all possible combinations (Table S1). By using the ten replicates for each grid cell and scenario, we determined the probability of each combination of features to either be converted into another combination, or to remain unchanged (Table 1). These determined probabilities were used as the foundation for the Markov chain MM adapted from Cipriotti et al. (2016), described 1 3 in the following section. The resulted data tables are openly accessible from the GitHub repository referenced in the data availability statement section.

Metamodel
For the development, first we imported the data table containing the combination of features (mean primary attractiveness, spruce proportion and number of infested trees) for all 7524 grid cells of the research area into R (Fig. 4). We then imported the processed simulation results of IPS-SPREADS (Table 1), which contain the possible destinations and corresponding transition probabilities for each combination of features.
In the next step, the MM checks if interactions between neighboring cells in the form of beetle spread should be considered or not (N is TRUE in Fig. 4). If this is the case, the model searches for all grid cells with at least one neighbor with an infested tree amount level of 286 (the highest possible amount). All of these cells simultaneously increase their level of infested tree amount by one level for each beetle generation simulated, according to the climate change scenario. With this MM extension, the spread of bark beetles from one grid cell to another is represented and directly correlated with the given amount of beetle generations during one year.
The baseline MM consisting of Markov chains based on (Cipriotti et al. 2016) is then applied (Fig. 5): Depending on the combination of features (primary attractiveness, spruce proportion and number of infested trees), each grid cell is assigned a new combination based on the simulation results with IPS-SPREADS and the chosen scenario of climate change and management (Table 1). If there is more than one possible destination for a combination of features, a random number is drawn to decide the new combination for the given grid cell according to the pre-calculated transition probabilities.
In the next step, the model checks if the grid cells' dieback of spruce trees should be considered or not (C is TRUE in Fig. 4). If this is the case, then the model applies the spruce dieback extension, in which the amount of living spruce trees is updated for all grid cells by removal of infested or felled trees. If the amount of living spruce trees reaches or falls below zero, then the corresponding cell is marked as 'dead' by setting mean primary attractiveness, spruce proportion and number of infested trees to 500 each. Furthermore, the MM calculates the proportion of dead trees over the total number of living spruces at the beginning of the simulation. This relation is later used for analyzing the impact of management and climate change scenarios on the long-term outbreak development in the research area.
The model then exports the data table containing the new combination of features for each grid cell as a CSV file, labeled with the climate change and management scenarios, and the simulated year.
At last, the MM increases the counter of simulated years by one (a + 1) and checks if the total number of years (i.e., 20) were already simulated. If not, the MM begins the next year by checking if the beetle spread extension should be applied (N is TRUE in Fig. 4) and so on. Once the total number of years is reached, the MM stops its calculations.

Model validation: reproduction of outbreak development in 2017 and 2018
To validate the performance of the developed MM and assess if the introduced beetle spread and spruce dieback extensions improve the predictions of outbreaks development, observations from outbreaks dating from 2017 and 2018 in the Saxon Switzerland national park were used as a measure of agreement. Data of both years were processed as described in Data processing and used to compare the MM predicted number of infested trees for each of the 7524 grid cells. To evaluate the goodness of agreement between model prediction and observations, the macro-averaged F1 score was applied. This score is defined as the harmonic mean of precision and recall, and is commonly used for evaluation of classification models or algorithms (Sammut and Webb 2010) as well as an optimization criterion (Yang 1999). As the F1 score is calculated for each class separately, the global mean of the per-class F1 scores (macro-averaged F1 score) was used for comparison of model predictions (Yang 1999).
To also validate the introduced extensions of beetle spread and spruce dieback, four different MM variants were used to predict the infestation development   For each of these variants and for each possible combination of climate change and management scenarios a macroaveraged F1 score was calculated and visualized as tile plot. We then visually inspected the MM variants' results and data, and determined the variant with highest agreement. The intensity of sanitation felling was also varied during the metamodel validation because no information on location, timing or intensity of felling inside and outside the national park was available. Sanitation felling in the MM was applied in the same intensity on all grid cells inside and outside the national park alike.

Long-term forecast of outbreak development
For the long-term forecast of the outbreak development in the Saxon Switzerland national park, the MM variant with the highest overall macro-averaged F1 score during model validation (Section Metamodel Validation) was applied for 20 consecutive years. To investigate possible impacts of climate change and forest protection management on the outbreak development, the amount of yearly beetle generations was varied from two to three and the intensity of sanitation felling from zero to 100% in 10 percent steps. As the interaction between these two scenarios was also from interest, a total of 22 different combinations had to be simulated. A bar plot was generated to visualize the number of sites with 0, 25, 50, 75 or 100% of spruce dieback inside and outside the national park in dependency of the chosen combination of climate change and management scenario. With this bar plot it was possible to detect thresholds where sanitation felling was able to avoid the spread of infestations from the national park to adjacent forest areas or when no mass outbreak happened in the protected area at all. A map was generated for each scenario showcasing the proportion of dead spruces of each grid cell for the research area as an index of outbreak risks.

Model validation: reproduction of outbreak development in 2017 and 2018
The calculated macro-averaged F1 scores of the four different MM variants showed that the highest mean agreement between model prediction and observations was achieved by the model version with both beetle spread and spruce dieback extensions (Fig. 6). While the introduction of spruce dieback only marginally improved the prediction, the introduction of beetle spread raised the agreement substantially. The full model version however was used for the long-term forecast because of its capability to also assess both the impact of climate and management scenarios on the dieback of spruces inside and outside the national park. In general, the agreement between model prediction and observations was higher for 2017 than 2018, regardless of the chosen scenario or MM type. The validation of MM variants using the macro-averaged F1 score also revealed that scenarios with zero to 90% sanitation felling intensity achieved the highest agreement between model prediction and observations regardless of the chosen amount of beetle generations. Lowest agreements were reached when a sanitation felling intensity of 100% was applied.  Table S1). By using the transition probabilities and possible destinations of stages derived from simulation experiments with the IBM IPS-SPREADS (Table 1), new combinations of spruce proportion, mean primary attractiveness and infested trees are assigned to all grid cells of the research area per time step

Long-term forecast of outbreak development
Applying the MM variant with beetle spread and spruce dieback extensions for 20 years using climate and management scenarios revealed the possible risks of complete dieback of spruce forests in nearly all grid cells outside the national park under the warmer climate change scenario (i.e., three beetle generations per year) and without sanitation felling (Fig. 7). This risk was diminished by sanitation felling of 80% and more, leading to nearly no dead trees outside the national park. For the scenario with two beetle generations per year, even without management intervention only a small proportion of the grid cells located outside the national park were affected by infestations. Inside the national park, nearly half of the grid cells populated with spruces died completely when no management intervention took place. With increasing sanitation felling intensity, the proportion of completely dead grid cells inside the national park was reduced and reached zero percent when a felling intensity of 90 or more percent was applied. The proportion of grid cells with complete spruce dieback outside the protected area also decreased with increasing management intervention and reached zero with 80% felling intensity. A direct comparison of both climate scenarios (two and three beetle generations per year) revealed that an outbreak with complete dieback of all spruces in forests adjacent to the national park did not occur regardless of management intervention when only two beetle generations were simulated. In both scenarios, a felling intensity of 80% and more substantially decreases or diminishes the amount of dead spruces inside and outside the national park. While low management intervention activities had no effect on the outbreak with three beetle generations, the amount of dead spruces linearly decreased with sanitation felling intensity when two beetle Fig. 6 Macro-averaged F1 scores (harmonized mean of prediction and recall) for rating the agreement between metamodel predictions and observed European spruce bark beetle (Ips typographus) outbreak development in 2017 and 2018 in Saxon Switzerland national park, Germany. Columns represent the investigated years, rows the applied metamodel variants. The higher the scores are (the darker the blue color), the better the agreement is generations were present. The maps on spruce dieback risk generated from the 20-year simulation show hotspot areas with high risk for infestations spreading from the national park to adjacent forests located in the north of the eastern part as well as in the northeast of the western part of the protected area (Figs. 8, 9).

Discussion
In this study, we tackled a core problem of forest ecology and management: upscaling from tractable small-scale observations to forecasting large-scale patterns. We used the metamodel framework, MM, adapted from Cipriotti et al. (2016) to scale a computationally intensive yet mechanistic precise bark beetle model to a whole landscape encompassing a national park. Our study introduced a novel analysis of forest pests using simulation models. The approach focuses on defining important state variables (i.e., plot) using a small-scale model to design a state space that is partitioned into a finite number of discrete states. The transition probabilities from one state to the next (i.e., from a non-infested state to an infested state) are then determined by monitoring extensive simulation runs of the small model, which covers the entire range of initial conditions and state combinations. We introduced two extensions to the MM approach, which incorporate both the spread of beetles from one grid cell to another, as well as the possible dieback of spruce individuals With this MM we then investigated the impact of climate change and of management on the long-term development of outbreaks on a landscape level through an increase in yearly beetle generation numbers and sanitation felling, respectively.
The validation of four different MM variants revealed that the baseline model version (i.e., without beetle spread or spruce dieback extension) was not sufficient to hindcast the outbreak development of 2017 and 2018 in the national park. Particularly, the introduction of an indirect interaction of neighboring grid cells by the spreading of beetle infestations substantially improved model predictions. Sites with weakened spruces have been reported as step stones for successful infestation of neighboring and more robust areas, eventually leading to mass outbreaks at regional scale (Křivan et al. 2016), which can include highly resilient forests. Furthermore, connecting habitats has a particularly driving role in facilitating and producing mass outbreaks of European spruce bark beetle (Kärvemo et al. 2016;Seidl et al. 2016;Wildemeersch et al. 2019). The higher agreement between model scenarios and observations from 2017 and 2018 are compatible with two to three yearly beetle generations, suggesting that climate-induced changes to temperature and rainfall regimes in the region might have been already benefitting beetle outbreaks. An increased number of yearly beetle generations have the potential to lead to more disturbances and to a two-to threefold increase in damaged timber volume (Dobor et al. 2020;Seidl et al. 2008). This highlights the importance of forecasting climate change scenarios on both short-and long-term bark beetle outbreak development modeling. Even though the combination of beetle spread and spruce dieback extensions did not substantially improve the agreement between model prediction and observations in comparison to the use of the beetle spread extension alone, by combining both extensions we were able to depict and assess the risks of infestation of adjacent forests through the proportion of dead spruce trees per grid cell and the proportion of grid cells with completely dead spruces. Nevertheless, we suggest that the marginal impact of the spruce dieback extension on the agreement between models and observations might be related to the short simulation period and/or, more likely, to the lack of data on sites with 100% spruce dieback; grid cells with most spruce trees killed were still classified as the highest level of infestations and no consideration of complete dieback was done during the rating of agreement between model results and observations. Another possible reason for the discrepancy between model predictions and observations may lie on weather. The year of 2018 was one of the hottest and driest years observed and resulted in three beetle generations and two sister broods by the end of September (IFFF 2022). This might have promoted the ongoing mass-outbreak in the national park, leading to a potential underestimation of infestations by the MM. It is also possible, that the MM extension of beetle spread from grid cell to neighboring cells may not be enough to capture the speed and intensity of the outbreak development sufficiently. Most new infestations occur within 100 m (Kautz et al. 2011;Wichmann and Ravn 2001) to 500 m (Kautz et al. 2011;Stadelmann et al. 2013) around initial infested trees. Therefore, further research could focus on testing whether including beetle spreading from a high infested grid cell to both direct and secondary neighbors could improve agreement between model outputs and observations. While two annual bark beetle generations imposed only a marginal risk for forests adjacent to the national park, an increase to three generations per year due to drier and hotter years as associated with climate change might have shifted the predictions to complete diebacks within the national park (Dobor et al. 2020;Seidl et al. 2008), as well as to high infestation risks for adjacent forests. To address the outbreak risk and decrease potential risks for forests surrounding the national park, a timely removal of at least 80% of infested trees might be needed. Similar magnitudes for a successful reduction in infestation risk were found with ≥ 95% for salvage logging (Dobor et al. 2020) and with 80% for bark beetle mortality in general . Such a high reduction is difficult to achieve and only feasible if damage is limited to small areas (Dobor et al. 2020) as the removal of infested trees might not effectively stop a mass outbreak once it reached a large scale (Kulakowski 2016). Nonetheless, timely removal or debarking of infested trees is still considered to be the most efficient measure against the spread of bark beetles (Jönsson et al. 2012;Stadelmann et al. 2013;Wermelinger 2004). Our study results support this assessment: with 80% or more sanitation felling intensity the infestations are reduced substantially regardless of the amount of yearly beetle generations present. Although further studies could consider other measures used in bark beetle management, such as pheromone traps, we did not include it in our study as these might not be suited for the reduction in infestations, being mostly used for monitoring of seasonal flight activity (Dimitri et al. 1992;Lobinger and Skatulla 1996;Wichmann and Ravn 2001). The almost complete dieback of spruce trees inside the national park, with low sanitation felling and two or three beetle generations might be related to the large bark beetle population within parts of the national park where no management intervention took place. There, population densities are reached, that cause infestation of nearly all spruces of a given dimension regardless of climatic factors or individual tree conditions (Marini et al. 2017). As a consequence of this mass outbreak and the resulting reduction in spruce individuals at a regional scale, we suggest that a collapse of the beetle population could be expected, similarly to what has been observed for spruce beetle (Dendroctonus rufipennis) interacting with Engelmann spruce (Picea engelmannii) across the Southern Rocky Mountains of northwestern Colorado, USA (Hart et al. 2015).
A possible reason for the low infestation risk of adjacent forests with only two yearly beetle generations could lie within the missing connectivity of forested grid cells between the national park and its surrounding forests. Beetle spread was simulated in the MM as interaction between grid cells in direct neighborhood, which neglected the possibility of beetle flight from grid cells farther away. Furthermore, in reality sister broods happen regularly in national park region and lead to higher beetle densities and to more infestations than the IPS-SPREADS model or the MM predicted. In addition, as we used the infestation maps from 2016 as baseline for our model predictions, we had no information on the beetle population present outside the protected area at the beginning. Therefore, infestations predicted by the model could only originate from the infested trees inside the national park, which again leads to an underestimation of infestation risk in forests adjacent to the protected area. While our simulation results indicate possible reductions in infestation spread, sanitation felling in our model was implemented in a near perfect manner: detection and removal of infested trees always took place right after the last beetle of a generation finished its dispersal flight and right before sister brood establishment or brood emergence could take place. As this procedure is considered as using an optimal timing with the highest effectiveness (Kautz et al. 2022), our model application possibly overestimated the effects of sanitation felling. Further challenges influencing the effectiveness of sanitation felling include the delay in visual detection of infestations, limitation of technical and human resources, high costs of debarking or challenging terrain conditions (Wermelinger 2004). With such challenges, sanitation felling intensities practice are unlikely to reach 100% with level of even 75% not being feasible on larger scales (Dobor et al. 2020). Nevertheless, while lower felling intensities did not stop the beetle spread from the national park to adjacent forests entirely, sanitation felling can be a valid option for small high-risk areas that could potentially serve as a bridge for infestation spread to adjacent forests.
The implementation of a buffer zone at the border of the national park as an outbreak contention measure might be a next step to be tested within our proposed MM framework on upscaling an IBM on bark beetle outbreak. In other national parks in Germany it was shown that such a buffer zone between 100 m (Niemeyer et al. 1995) and 1500 m (Heurich 2001) could be enough to substantially reduce the spread of infestations from protected to adjacent forest areas (Wermelinger 2004). This buffer zone could be implemented in the model to test the impact of different sanitation felling intensities in this zone as well as different widths of the buffer on the beetle spread. Implications from such modeling studies could then be used for the national park management to identify areas with high risks and to focus management actions accordingly.