Combining an Occurrence Model and a Quantitative Model for the Prediction of the Sanitary Felling of Norway Spruce Because of Bark Beetles

Abstract

The European spruce bark beetle (Ips typographus L.) is an eruptive forest pest that has caused a great deal of damage in the last decades because of increasing climatic extremes. In order to effectively manage outbreaks of this pest, it is important to predict where they will occur in the future. In this study we developed a predictive model of the sanitary felling of Norway spruce (Picea abies (L.) H. Karst.) because of bark beetles. We used a time series of sanitary felling because of bark beetles from 1996 to 2020 in Slovenia. For the explanatory variables, we used soil, site, climate, geographic, and tree damage data from the previous year. The model showed that sanitary felling is negatively correlated with slope, soil depth, soil cation exchange capacity, and Standard Precipitation Index (less sanitary felling in wet years). On the other hand, soil base saturation percentage, temperature, sanitary felling because of bark beetles from the previous year, sanitary felling because of other abiotic factors from the previous year, and the amount of spruce were positively correlated with the sanitary felling of Norway spruce due to bark beetles. The model had an R² of 0.38. A prediction was performed for 2021 combining an occurrence model and a quantitative model. The model can be used to predict the amount of sanitary felling of Norway spruce due to bark beetles and to refine the risk map for the next year, which can be used for forest management planning and economic loss predictions.

1. Introduction

Bark beetles pose a significant threat to forest health in the light of climate change. Many large-scale outbreaks have been recorded in Europe and North America [1,2], with dramatic economic and carbon storage losses [3]. The main reasons for these outbreaks are the type of forest management, invasion in new regions, and large-scale catastrophic climatic events [1,4,5,6,7,8]. It is therefore important to develop an integrated system that takes all phases of the population dynamics of bark beetles into account in order to combat outbreaks [9].

The European spruce bark beetle (Ips typographus L. 1758) causes more damage to European temperate forests than any other forest pest [2,10]. The main host plant is Norway spruce (Picea abies (L.) H. Karst.), which is the most common tree species in Europe and has important ecological and economic functions. Due to planting regimes in secondary habitat and increasing temperature associated with many weather extremes, trees are weakened and therefore susceptible to large-scale spruce bark beetle outbreaks. An integrated management approach is therefore suggested for this species [2,9]. Damage intensity can be reduced by a monitoring system for the early detection and survey of spruce bark beetle outbreaks. When an outbreak is found, direct methods should be used to reduce populations with the help of sanitation felling, trap trees, and pheromone traps. It is important to take into account that these methods should be combined and used in the correct stage of the phenology of the spruce bark beetle. In addition to direct methods, silvicultural practices that increase the forest’s resistance to bark beetles should be used [6,8,11,12]. One of the crucial aspects of direct methods is to detect an outbreak as quickly as possible. Therefore, detection methods should be further developed.

Predictive models can be helpful for determining where to focus surveys of spruce bark beetles [13]. There are phenological models that predict bark beetle phenology, such as PHENIPS [14] and RITY [15]. On the other hand, there are ecological models which take different variables into account [13,16,17,18]. These models have found that outbreak potential is most affected by drought and high temperatures, the amount of spruce and storm felled trees [5,19,20]. However, the majority of these models do not predict future outbreaks in the short term or long term (for example, see [16,17,20,21]). Such models can focus on predicting the occurrence of a potential outbreak or predicting the amount of sanitary felling. Both outcomes are important in their own way. The prediction of occurrence gives foresters the opportunity to locate areas which are under threat, while the prediction of sanitary felling indicates how severe an outbreak will be.

Slovenia has a long history of bark beetle outbreaks. Several methods exist for the short- and long-term predictive modeling of spruce bark beetles [15,16,17,22] and have been implemented for several years [22,23,24]. However, a model that predicts short-term sanitary felling is still lacking.

The aim of this study was to predict the amount (in m³) of sanitary felling of spruce in Slovenia because of spruce bark beetles. This involved three steps: first, we created a model of sanitary felling because of bark beetles on a test dataset based on different soil, geographic, climate, and stand variables. Next, we validated the model on a validation dataset. Finally, we created a risk map of sanitary felling because of bark beetles for the year 2021.

2. Materials and Methods

2.1. Area Description

Slovenia is a small country in Central Europe. It lies at the intersection of the sub-Mediterranean, Dinaric, sub-Pannonian, and Alpine biogeographical regions. More than 60% of Slovenia is covered with forests, which in 2020 consisted mainly of beech (32.9%), Norway spruce (30.2%), and silver fir (7.4%) [25]. The country has a continental climate with hot summers and cold winters, but the coastal area experiences milder winters. In the last few decades there have been several large-scale climatic catastrophes. In 2003, there was a severe drought which resulted in a long period of bark beetle outbreaks. In 2014, a catastrophic ice storm catalyzed severe bark beetle outbreaks across the entire country. Despite many small-scale windthrows in the last few years, bark beetle abundance dropped in 2020 [25,26] (Figure 1).

2.2. Databases

For the model development we used the same database as in de Groot and Ogris [17] but instead of the relative amount (cubic meters per hectare), the absolute amount (cubic meters) of Norway spruce and sanitary felling was used.

For the dependent variable sanitary felling of spruce because of bark beetle outbreaks in the current year (cubic meters) was used. Independent variables were as follows: amount of spruce (cubic meters), altitude (meters above sea level), relief aspect (degrees), slope (degrees), dominant geology type, exchangeable phosphorus (mg/100 g), soil depth (centimeters), soil cation exchange capacity (mmolc/100 g), soil base saturation (percent) [16,27], Standardized Precipitation Index (SPI—a proxy for drought) during the previous growing season, average temperature during the previous year (degrees Celsius), cumulative amount of precipitation in the previous year (millimeters), amount of sanitary felling of spruce because of bark beetle outbreaks in the previous year (cubic meters), amount of weakened trees because of bark beetles in the previous year (cubic meters), amount of sanitary felling because of abiotic factors in the previous year (cubic meters), and the location (X and Y) of the model grid cell.

The spatial resolution of the data was 1 × 1 km (21,001 model cells). Sanitary felling variables were extracted from the Timber database for the period 1996–2020, which was provided by the Slovenia Forest Service [28]. The Timber database distinguishes between sanitary felled trees due to bark beetle outbreaks and felled trees due to regular bark beetle attacks of weakened trees. Sanitary felling is the felling of trees in the outbreak area, where trees are damaged primarily because of bark beetle attack. It is used to prevent the further spread of the outbreak to surrounding areas. Felling of weakened trees is used when trees are weakened primarily because of other factors such as drought and other weather extremes which make them susceptible to secondary bark beetle attacks. Every felled tree is recorded at the forest sub-compartment level, which measures 22 ha on average (N > 53,000). Additionally, the number of trees and their volume in m³ in different size classes are also recorded, along with the reason (62 possible values) for felling and the causative agent, e.g., I. typographus.

Altitude data were provided by the Surveying and Mapping Authority of the Republic of Slovenia [29]. Data on the amount of Norway spruce were also provided by the Slovenia Forest Service in the Forest Funds database [30]. Exposition and slope were calculated with ESRI ArcGIS Desktop 10.6.1 on the basis of a digital elevation model with a spatial resolution of 12.5 m [29]. Air temperature and precipitation were acquired from the Slovenian Environment Agency [31,32]. Soil data on exchangeable phosphorus (mg/100 g), soil depth (centimeters), soil cation exchange capacity (mmolc/100 g) and soil base saturation (percent) were prepared on the basis of a soil map [33] by Ogris [27].

2.3. Analysis

For the analysis only the data in which Norway spruce was present were taken into account. The dataset was divided into a training dataset (80% of the data) and a validation dataset (20% of the data). First, the data were explored for outliers, multicollinearity, types of relationships and spatial autocorrelation, using plots to determine outliers and relationship types, the variance inflation factor for multicollinearity and Moran’s I for spatial autocorrelation [34]. A linear mixed model was used for the analysis. The dependent variable was the amount of sanitary felling because of bark beetles in the current year. The independent variables were the amount of spruce, altitude, exposition, slope, amount of exchangeable phosphorus, soil depth, soil cation exchange capacity, soil base saturation percentage, SPI during the previous growing season, average temperature during the previous year, amount of sanitary felling of spruce because of bark beetles in the previous year, amount of weakened trees of spruce because of bark beetles in the previous year, amount of sanitary felling of spruce because of abiotic factors in the previous year, and geographical coordinates. Due to outliers, sanitary felling, slope, sanitary felling because of bark beetles, sanitary felling because of abiotic factors, and amount of spruce were log + 1 transformed. Soil depth was square root transformed. The year and grid ID were used as random effects. As temperature and altitude were multicollinear, altitude was dropped from the analysis. Other variables were not multicollinear. There was no spatial autocorrelation observed.

Before starting model selection, a complete model was built including all of the variables. Model selection was based on the Akaike Information Criterion (AIC). The best model was the model with the lowest AIC. For the validation, the marginal R² and conditional R², RMSE, MSE, and MAE were calculated for the validation dataset in order to estimate how much of the variability was covered by the model. Lastly, the residuals were checked for normality.

Predictions for 2021 were based on the final quantitative model and the probabilistic model developed by de Groot and Ogris [17] and applied already for 2021 [23] (Figure 2). We used 2020 data for all the variables which were included in the model. The prediction was performed on the basis of fixed effects, and random effects were omitted. The final map was prepared in two steps. First, the probabilistic predication was prepared to show the potential for sanitary felling above the threshold. Validation of this model showed that the optimal threshold for the prediction of actual sanitary felling occurrence was 0.55 [35]. For the cells with a threshold of more than 0.55, the amount of sanitary felling was calculated with the model described here. For the prediction, the minimum, maximum, mean, and sum were calculated.

For the analysis, the statistical program R [36] and the libraries “LME4” [37], “car” [38], “spdep” [39,40], “adespatial” [41], “spatialreg” [39,40], and “MuMIn” [42] were used.

3. Results

The best model included the amount of Norway spruce, slope, amount of sanitary felling of spruce in the previous year, amount of sanitary felling because of abiotic factors, SPI, temperature, soil depth, soil cation exchange capacity, and soil base saturation percentage (

Table 1. Model of the prediction of the amount of sanitary felling because of bark beetles.

Parameters	Estimate	Std Error	df	T Value	p
(Intercept)	−2.2044	0.1219	53.0850	−18.0831	2.38 × 10−24
log(1 + slope)	−0.0789	0.0087	12,250.3280	−9.1134	9.22 × 10−20
sqrt(soil depth)	−0.0229	0.0036	11,199.1970	−6.3855	1.78 × 10−10
Soil cation exchange capacity	−0.0104	0.0010	12,178.9620	−10.6647	1.95 × 10−26
Soil base saturation percentage	0.0067	0.0004	12,110.8650	15.6774	7.47 × 10−55
Temperature	0.0852	0.0043	15,433.4990	19.6619	5.01 × 10−85
SPI	−0.2479	0.0088	110,276.0280	−28.0173	4.04 × 10−172
log(1 + sanitary felling because of bark beetles)	0.4052	0.0026	149,932.8310	153.9234	0.00
log(1 + sanitary felling because of abiotic factors)	0.2092	0.0026	159,012.0310	79.8858	0.00
log(1 + amount of spruce)	0.2996	0.0038	42,934.0320	77.8883	0.00

). The amount of spruce, temperature, amount of spruce in the previous year, amount of sanitary felling because of abiotic factors, and soil base saturation percentage were positively correlated with sanitary felling because of bark beetles, while slope, soil depth, soil cation exchange capacity, and SPI were negatively correlated with sanitary felling because of bark beetles. The marginal R² was 0.34 and the conditional R² was 0.46 (

Table 2. Validation measures of the prediction model.

Validation Measures	Marginal Model	Conditional Model
RMSE	1.50	1.39
MSE	2.25	1.94
R2	0.38	0.46
MAE	1.23	1.12

).

The predicted sanitary felling for 2021 showed that the minimum sanitary felling was 1 m³ and the maximum was 780 m³ per 1 × 1 km cell, with a median of 3 m³. It was predicted that there would be 243,655 m³ of sanitary felling because of bark beetles for the whole of Slovenia in 2021. The combined prediction with the probabilistic model where only cells with the probability of sanitary felling of Norway spruce due to bark beetles were taken into account was 198,721 m³. The highest amount of sanitary felling was predicted to be in the northwest, north, and south of the country. The actual amount of sanitary felling of spruce due to bark beetles was 261,288 m³ (as of 31 December 2021). Differences at the model cell level were also low, at only 0.02 m³ on average (RMSE = 57.7 m³). Differences at the model cell level for the combined model were moderate, at −10.5 m³ on average (RMSE = 114.8 m³) (Figure 3). Overall, there were no large spatial differences between predicted and actual sanitary felling. Large differences between the predicted and actual amount of sanitary felling were observed only in the west and north of Slovenia.

4. Discussion

We created a model to predict the sanitary felling of spruce because of bark beetles. The model included slope, soil depth, soil cation exchange capacity, soil base saturation percentage, temperature, SPI, sanitary felling in the previous year, and the amount of spruce. The model performed very well considering that it was prepared in an ecological system with a large amount of background noise. Combining the probability and amount model predicted that there would be a relatively low amount of sanitary felling in 2021 in most parts of Slovenia.

The model previously developed for Slovenia for the occurrence of sanitary felling is very similar to our new model of the amount of felling because of bark beetles ([17]; this study,

Table 3. Estimates of the occurrence model and amount of sanitary felling model. * The variable was not transformed. Occurrence model variables were taken from de Groot and Ogris [17].

Variable	Estimates
	Occurrence Model	Amount Model
Intercept	−2.41	−2.2044
log(1 + slope)	−8.21 × 10−3	−0.0789
sqrt(soil depth)		−0.0229
Soil cation exchange capacity	−1.77 × 10−2	−0.0104
Soil base saturation percentage	1.14 × 10−2	0.0067
Temperature	2.66 × 10−1	0.0852
SPI	−3.88 × 10−2	−0.2479
log(1 + sanitary felling because of bark beetles)	1.38	0.4052
log(1 + sanitary felling because of other abiotic factors)	5.50 × 10−1	0.2092
log(1 + amount of spruce)	2.37 × 10−2 *	0.2996
cbind(x, y)x	−4.43 × 10−6
cbind(x, y)y	4.00 × 10−6
Phosphorus	−5.39 × 10−2
log(1 + amount of weakened spruce in the previous year)	5.04 × 10−2

). There are some differences in the estimate values of the parameters, but this has to do with the dependent variables: 0–1 and the amount of felling. However, for some data there was a difference in the transformation of the variables. This was probably due to additional data used in the amount model: the occurrence model was based on data up to 2016, and the amount model was based on data up to 2020. Interestingly, the models also differed with respect to the variables. For instance, the occurrence model included as variables phosphorous in the soil and sanitary felling of weakened trees, while the amount model included soil depth. The occurrence model is a probabilistic model which shows the probability of the occurrence of sanitary felling, while the amount model predicts how much sanitary felling will be performed in the next year. This difference in models also explains why certain variables are included in the occurrence model or in the amount model. For instance, there is a high chance that an attack will still occur in weakened trees as they attract bark beetles. For the model of the amount of sanitary felling, these weakened trees are less important. The amount of sanitary felling of weakened trees is not very high and therefore does not strongly affect timber volume, while it still indicates the occurrence of attacks. The influence of soil variables such as soil depth and phosphorus are different for both models. The fact that phosphorus was only found in the occurrence model could imply a negative association between phosphorus (representative also for other nutrients such as nitrogen and potassium) and spruce bark beetle outbreaks and that phosphorus potentially increases tree vigor [43]. This results in higher demand for water and other nutrients, which increases tree stress and consequently increases the likelihood (but not necessarily the amount) of sanitary felling. On the other hand, soil depth was an important factor for the amount of sanitary felling. As soil depth reflects the availability of nutrients, water storage and potential for root development [44,45], it seems that a lack of nutrients, less water storage, and poor root development are associated with an increase in the amount of sanitary felling, especially in extreme weather events such as drought and windthrows.

Precipitation and SPI are important factors influencing the prediction of the sanitary felling because of bark beetles. Previous research has shown that outbreaks are not necessarily dependent on windthrows or the ice storm, but can also be due to acute drought stress which affects host trees [46]. This also happened in the years after the 2003 drought in Slovenia (see also Figure 1).

The slope was negatively associated with the sanitary felling because of bark beetles. There are two possible explanations for this: the poor accessibility of steep, high-altitude stands results in less sanitary felling in these regions. On the other hand there could also be fewer attacks in spruce stands in mountainous/sub-alpine areas because of the less appropriate climate for spruce bark beetles [4].

Other studies have found that forest structure and species composition; forest management; climate; soil; tree vitality; predisposition to storm, ice and snow damage; and species-specific phenology affect bark beetle damage [5,6,11,19,47,48,49]. These results are comparable with our results. Pasztor et al. [11] developed a similar model with a 10-year time series in Austria which also used sanitary felling data. We had the opportunity to use a 25-year time series in neighboring Slovenia. The main difference between the models is that climate data were not included in the model of Pasztor et al. [11]. However, in that model, temperature is indirectly included in the number of generations and the altitude. We decided not to include the number of generations or altitude because the number of generations is a derivate of temperature. Altitude is strongly correlated with temperature but is static. It has been shown that previous year temperature influences bark beetle outbreaks [47]. Therefore, it is better for prediction models to include dynamic variables such as temperature. In our model, forest structure was represented by the amount of spruce, but species composition and stand age are also important [11,49] and could improve the model when these data are available. However, the predictive power of our model is quite high. We found that the model explained 38% of the variability, which is better than the model from Austria (13%). Furthermore, the random effects of site and year accounted for only 8% of the variability, which means that our variables were very well chosen. Given that there was a great deal of ecological noise in the data, the total explained variability of 0.46 is very good.

Implications for Forest Management

The combined model with occurrence and amount of sanitary felling predicts sanitary felling on 1 × 1 km cells over the whole area of Slovenia. As shown, this yields a map of predicted sanitary felling in the next year, which dramatically reduces the area, and therefore a more targeted search for bark beetle infested trees can take place. The prediction of the amount of sanitary felling makes it possible to (1) refine the risk map, (2) estimate the economic cost of the damage in the next year, and (3) estimate the spatial distribution of predicted sanitary felling. As already mentioned above, the occurrence map gives foresters an indication of where to search for attacked trees in general irrespective of the number of attacked trees on an area of 1 by 1 km. This is very important, because just a few attacked trees have the potential to become large outbreaks if they are not found early enough. The map produced based on the model of the amount of sanitary felling shows the areas which are likely to be heavily affected by bark beetles. Foresters could therefore prepare for larger outbreaks well in advance (up to half a year). This is important, as normal forest management is not organized to carry out the intensive sanitary felling required after large-scale outbreaks. Another aspect for which the map is a helpful tool is cost estimation and actual cutting. The amount of sanitary felling reflects the lost yield of Norway spruce and therefore can be taken into account for cost and yield estimation for the next year. On the other hand, the activity of sanitary felling, especially in the case of large-scale outbreaks, might affect normal forest management. Such a map would indicate which areas will not achieve the cutting quota. Therefore, the preparation of a map of the predicted amount of sanitary felling could be beneficial in a number of ways.

5. Conclusions

In conclusion, we developed a model to predict the amount of sanitary felling of spruce because of bark beetles. Many climate, forest structure, site specific, and biotic and abiotic disturbance variables were included. This resulted in a model which, despite the ecological background noise, explained a relatively high proportion of the variability. Developing models and risk maps is a very important aspect of forest pest management as they can enable foresters to respond quickly to potential outbreaks. In the era of climate change and elevated probability of large-scale bark beetle outbreaks, such prediction maps are of vital importance to reduce the risks of bark beetles, especially in large areas where outbreaks might be missed.