Bark beetle damage in Norwegian forests: a study of model suitability and projected impact under climate change

ABSTRACT

Bark beetle (Ips typographus) outbreaks have the potential to damage large areas of spruce-dominated forests in Scandinavia. To define forest management strategies that will minimize the risk of bark beetle attacks, we need robust models that link forest structure and composition to the risk and potential damage of bark beetle attacks. Since data on bark beetle infestation rates and corresponding damages does not exist in Norway, we implement a previously published meta-model for estimating I. typographus damage probability and intensity. Using both current and projected climatic conditions we used the model to estimate damage inflicted by I. typographus in Norwegian spruce stands. The model produces feasible results for most of Norway’s climate and forest conditions, but a revised model tailored to Norway should be fitted to a dataset that includes older stands and lower temperatures. Based on current climate and forest conditions, the model predicts that approximately nine percent of productive forests within Norway’s main spruce-growing region will experience a loss ranging from 1.7 to 11 m³/ha of spruce over a span of five years. However, climate change is predicted to exacerbate the annual damage caused by I. typographus, potentially leading to a doubling of its detrimental effects.

KEYWORDS:

• Ips typographus
• Picea abies
• forest management
• model assessment
• boreal forest
• Norway spruce

Introduction

Norway spruce, Picea abies (L.) Karst. (Pinales: Pinaceae), dominates Scandinavian forests (Schlyter et al. 2006) and is a key species for the Norwegian forest sector. This commercially important timber species is often more productive than other native tree species and grows well under a wide range of climatic conditions (Jandl 2020). However, Norway spruce is also the main host of the European spruce bark beetle, Ips typographus (L.) (Coleoptera: Curculionidae: Scolytinae) (Christiansen and Bakke 1988; Jeger et al. 2017), one of the most devastating biotic disturbance agents in European forests. In Norway, Norway spruce is the most economically important tree species, constituting 42% of the growing stock, and approximately 80% of the annual cuts (Breidenbach et al. 2021). Thus, a large portion of Norwegian forests is currently susceptible to spruce bark beetle infestation.

Bark beetles can kill millions of trees during intermittent outbreaks and thus have substantial negative impacts on forests. In addition to increasing tree mortality rates, outbreaks also reduce forest stand densities and disrupt the carbon balance (Abdullah et al. 2018; Kosunen et al. 2020; Nowakowska et al. 2020). Outbreaks may also affect the socioeconomic conditions of forest owners and communities that depend on the forest, e.g. by reducing the timber value of attacked trees and increasing forest management costs (Schowalter 2012), as well as causing a short-term oversupply of wood that reduces timber prices (Stereńczak et al. 2019; Hlásny et al. 2021a).

Most of the time, Ips typographus breeds in dead or severely weakened trees; however, when beetle densities cross a threshold, they can kill healthy spruce trees through pheromone-coordinated mass attacks (Hlásny et al. 2019). A given population of bark beetles can be in either of two states – endemic or epidemic. The endemic state is characterized by population densities so low that the beetles cannot overcome the defenses of healthy spruce trees. In the epidemic state, the beetles can successfully colonize healthy trees and sustain an outbreak. This will further increase bark beetle densities, resulting in a self-sustaining outbreak that can last for years and kill millions of trees. Importantly, the epidemic threshold separating the endemic and epidemic states is not determined solely by beetle density, but rather by a combination of beetle density and tree resistance/vigor. When tree vigor is reduced by, e.g. drought or wind damage, fewer bark beetles are required for successful tree colonization and initiation of an outbreak.

In recent times, there has been a significant rise in bark beetle-induced tree mortality in Europe, particularly in spruce monocultures established outside their natural range (Hlásny et al. 2021a). While Europe saw an annual average of 2.9 million m³ of timber loss to bark beetles between 1950 and 2000 (Schelhaas et al. 2003), the situation intensified from 2018 to 2020: southern Sweden experienced 17 million m³ of infested spruce timber (Schroeder and Fritscher 2020; Müller et al. 2022), while the Czech Republic and Germany faced losses of 50 million m³ and 86 million m³ of spruce, respectively (Fernandez-Carrillo et al. 2020; Statistisches Bundesamt 2021; Hlásny et al. 2021b; VÚLHM 2023). This surge in mortality likely reflects the negative impacts of drought and temperature stress on tree vigor following the exceptionally warm summer of 2018.

In southern Norway, Ips typographus killed around 5 million m³ of spruce trees from 1971 to 1981, facilitated by extensive storms and extreme droughts that lowered trees’ resistance (Bakke 1989). In some municipalities more than 20% of the spruce growing stock was lost (Worrell 1983). Although the outbreak subsided in 1981, Norway’s I. typographus populations persisted endemically (Worrell 1983; Økland et al. 2021), still causing some economic damage by killing weakened trees. This problem is now compounded by climate change, which will likely increase the risk of new bark beetle epidemics in Norway. For instance, in 2018 and 2022, the region in southern Norway that had previously suffered the 1971–1981 outbreak had very low early summer rainfall – only 10% of average precipitation in some places – along with high summer temperatures (Bakke et al. 2020; Meteorologisk Institutt 2022).

The number of completed I. typographus generations per year varies geographically. In lowland continental Europe and southern Scandinavia, warmer temperatures enable up to two or three generations per year, while a single generation still is the norm in the rest of Scandinavia (Lindman et al. 2023). Current temperature conditions allow for two beetle generations (bivoltinism) only during warm summers in coastal areas of southern Norway (Wermelinger and Seifert 1998). However, climate change, with its longer and warmer growing seasons, is expected to facilitate more frequent occurrences of a second generation (Jönsson et al. 2011), leading to more attack periods on trees and faster population growth.

Despite the growing threat, no studies have analyzed nor estimated bark beetle damage risk in Norwegian forests. Recent observations indicate a surge in beetle numbers to levels not seen since 1997 in some areas (Økland et al. 2021). The situation is aggravated by the possible transition from one to two beetle generations per year in southern Norway (Økland et al. 2021). To estimate current and future bark beetle damage potentials in Norway, we need a model linking climatic and forest conditions with bark beetle damage. However, building a new model or even refitting an existing model to Norwegian conditions require data on bark beetle infestation rates or number of trees killed by the insect. Unfortunately, systematic data collection on the volume of trees killed by beetles has never been implemented in Norway.

In such cases, it is common to utilize models developed in countries with similar climatic and biological conditions. Several models exist for estimating bark beetle damage or susceptibility in boreal and temperate forests (Müller et al. 2022; Romeiro et al. 2022; Nordkvist et al. 2023), but most of them require input data that are not currently available in Norway (such as information on previous damage). Among the available models, the Austrian meta-model by Seidl et al. (2009) is calibrated for I. typographus and can be used to estimate bark beetle damage at a whole-country level. The model can be run with inputs that are available for Norway, it can estimate probability and intensity of damage, and it can also address climate change (Romeiro et al. 2022). The model uses a two-stage modelling approach in which the probability of bark beetle damage (pBB; probability that one or more trees in a stand are killed by bark beetles) and damage intensity (iBB; the share of spruce killed by bark beetles) are estimated separately. The meta-model was parametrized in Austria for a range of ecological conditions – from lowland areas to the high Alps – and was validated using Austrian National Forest Inventory data and bark beetle damage information extracted from province-level harvest records. The model has also been applied to the whole of Europe (Hlásny et al. 2021a).

Given this model's applicability and success in varied European conditions, it presents a promising tool for our research. Therefore, our study aims to (1) evaluate if the I. typographus damage model by Seidl et al. (2009) is suited to Norwegian conditions, and, if it is, (2) utilize the model to estimate current and climate change-related bark beetle damage in Norway.

Materials and methods

To test if Seidl’s model is suitable for Norwegian conditions, we first conducted a sensitivity analysis of the key climatic and forest variables included in the model. Subsequently, probabilities of bark beetle damage (pBB) and damage intensity (iBB) were modelled for current climate and climate change conditions. Lastly, model results were compared to trap catches from the Norwegian bark beetle monitoring program and to the outputs of an I. typographus developmental model.

The study area

The present study covered mainland Norway and islands near the coast. Much of Norway has large topographical variation over short distances. Major vegetation zones in Norway include deciduous-dominated nemoral and boreonemoral zones in the south and the boreal zone where conifer woodlands dominate. In addition, alpine vegetation dominates at higher altitudes above the climatic woodland limit (Moen 1999). Historically, average annual temperatures in Norway vary from −2°C in northern regions and the high mountains, to 8 °C in southern and western parts of the country. Annual precipitation ranges from 300 to 3500 mm (Beylich 2021).

We used the circular sample plots (250 m²) measured by the Norwegian National Forest Inventory (NNFI) to estimate annual probability of bark beetle damage and damage intensity. These plots are systematically distributed over the entire landmass of Norway (Breidenbach et al. 2021). We used the most recent NNFI period from 2015 to 2019, where 20% of the total number of plots are measured each year in the 5-year span. Initially, we considered all 9832 plots with productive forest (production >1 m³/ha/year). Plots dominated by forest younger than 40 years and plots with <10% Norway spruce share were excluded from further analysis since bark beetle damage are insignificant under such conditions (Wermelinger 2004). Furthermore, plots that were split into subplots (due to varying forest conditions within the plot) were also excluded from further analysis. After these exclusions, 2924 plots (30% of the total) remained to be used in the analysis.

We divided the country into seven geographical regions (Figure 1): Nord-Norge (182 plots), Trøndelag (595 plots), Innlandet (960 plots), Østviken (292 plots), Vestviken (618 plots), Sørlandet (171 plots), and Vestlandet (107 plots). Since the most important regions for spruce production in Norway are Trøndelag, Innlandet (previously Oppland and Hedmark counties), and Vestviken (Buskerud, Vestfold, and Telemark counties), we focus the results of this study on those three regions (2173 plots in total).

Figure 1. The seven regions considered in this study (1: Nord-Norge; 2: Trøndelag; 3: Innlandet; 4: Østviken; 5: Vestviken; 6: Sørlandet; 7: Vestlandet). Sample plots from the Norwegian National Forest Inventory (n = 2924) are shown as pink points.

Models of bark beetle damage

The two-stage statistical meta-model from Seidl et al. (2009) was fitted to bark beetle damages at stand level as simulated in the hybrid forest patch model PICUS v1.41 (Lexer and Hönninger 2001; Seidl et al. 2005), which includes a sub-model of bark beetle induced mortality of Norway spruce (Seidl et al. 2007). For annual probability of bark beetle damage ($pB{B}_i$, the probability that one or more trees in a stand will be killed by bark beetles in year i), the original version combined two components: the PAS system (Netherer and Nopp-Mayr 2005), which evaluates the predisposition of forest stands to bark beetle outbreaks, and the PHENIPS model for predicting potential bark beetle generations (Baier et al. 2007), and was later integrated into the PICUS hybrid forest patch model. The pBB model used here (Eq. 1(1) $\begin{array}{rl}pB{B}_i& ={\displaystyle \frac{e^{Z_{jklmn}}}{1+e^{Z_{jklmn}}}}\\ Z_{ijklm}& =B0+B1\cdot \mathrm{LogTem}{\mathrm{p}}_j\\ & +B2\cdot \mathrm{LogPreci}{\mathrm{p}}_k+B3\cdot \mathrm{Ag}{\mathrm{e}}_l\\ & +B4\cdot \mathrm{Densit}{\mathrm{y}}_m+B5\cdot \mathrm{Shar}{\mathrm{e}}_n\\ & +B6\cdot (\mathrm{LogTemp}\cdot \mathrm{LogPrecip}{)}_{jk}\\ & +B7\cdot (\mathrm{LogTemp}\cdot \mathrm{Age}{)}_{jl}\\ & +B8\cdot (\mathrm{LogTemp}\cdot \mathrm{Density}{)}_{jm}\\ & +B9\cdot (\mathrm{LogTemp}\cdot \mathrm{Share}{)}_{jn}\\ & +B10\cdot (\mathrm{LogPrecip}\cdot \mathrm{Age}{)}_{kl}\\ & +B11\cdot (\mathrm{LogPrecip}\cdot \mathrm{Density}{)}_{km}\\ & +B12\cdot (\mathrm{LogPrecip}\cdot \mathrm{Share}{)}_{kn}+{\epsilon }_{jklmn}\end{array}$(1) ) is, therefore, a meta-model fitted to PICUS simulations that accounted for various combinations of climate and stand conditions.

For SMI (Soil Moisture Index over the growing season), the initial computation in Seidl et al. (2007) involved using both actual and potential evapotranspi­ration. For the model of annual damage intensity per infestation ($iB{B}_i$; the percentage of spruce killed by bark beetles in year i), SMI was integrated into a Stand Hazard Index (SHI) along with other factors such as a Stand Edge Index (SEE) and Host Tree Share (HTS). Due to the complexity of using evapotranspiration variables for large-scale applications, the initial SMI model formulated in Seidl et al. (2007) was later translated into the final generalized linear model SMI presented here (Eq. 2(2) $\begin{array}{rl}iB{B}_i& ={\displaystyle \frac{1}{1+e^{3.9725-2.9673\cdot \mathrm{SHI}}}}\\ \mathrm{SHI}& =\mathrm{SMI}\cdot \mathrm{HTS}\\ \mathrm{SM}{\mathrm{I}}_{jk}& =B0+B1\cdot \mathrm{Tem}{\mathrm{p}}_j+B2\cdot {\mathrm{Precip}}_k^{-2}\\ & +B3\cdot (\mathrm{Temp}\cdot \mathrm{Preci}{\mathrm{p}}^{-2}{)}_{jk}+{\epsilon }_{jk}\end{array}$(2) ; Seidl et al. 2009), where the main predictors are mean annual temperature and precipitation.

The final annual damage volume ($vB{B}_i$) was calculated as the product of $pB{B}_i$ and $iB{B}_i$ (Eq. 3(3) $vB{B}_i=pB{B}_i\cdot iB{B}_i$(3) ). (1) $\begin{array}{rl}pB{B}_i& ={\displaystyle \frac{e^{Z_{jklmn}}}{1+e^{Z_{jklmn}}}}\\ Z_{ijklm}& =B0+B1\cdot \mathrm{LogTem}{\mathrm{p}}_j\\ & +B2\cdot \mathrm{LogPreci}{\mathrm{p}}_k+B3\cdot \mathrm{Ag}{\mathrm{e}}_l\\ & +B4\cdot \mathrm{Densit}{\mathrm{y}}_m+B5\cdot \mathrm{Shar}{\mathrm{e}}_n\\ & +B6\cdot (\mathrm{LogTemp}\cdot \mathrm{LogPrecip}{)}_{jk}\\ & +B7\cdot (\mathrm{LogTemp}\cdot \mathrm{Age}{)}_{jl}\\ & +B8\cdot (\mathrm{LogTemp}\cdot \mathrm{Density}{)}_{jm}\\ & +B9\cdot (\mathrm{LogTemp}\cdot \mathrm{Share}{)}_{jn}\\ & +B10\cdot (\mathrm{LogPrecip}\cdot \mathrm{Age}{)}_{kl}\\ & +B11\cdot (\mathrm{LogPrecip}\cdot \mathrm{Density}{)}_{km}\\ & +B12\cdot (\mathrm{LogPrecip}\cdot \mathrm{Share}{)}_{kn}+{\epsilon }_{jklmn}\end{array}$(1) where:

• $pB{B}_i$ – annual probability of bark beetle damage

• Z_jklmn – linear combination of predictor variables

• LogTemp_j – logarithmic mean annual temperature

• LogPrecip_k – logarithmic annual precipitation

• Age_l – stand age

• Density_m – stocking density relative to fully stocked yield table stands

• Share_n – host tree share

• ${\epsilon }_{jklmn}$ – error term

• B0 to B12 – see Table 1 (2) $\begin{array}{rl}iB{B}_i& ={\displaystyle \frac{1}{1+e^{3.9725-2.9673\cdot \mathrm{SHI}}}}\\ \mathrm{SHI}& =\mathrm{SMI}\cdot \mathrm{HTS}\\ \mathrm{SM}{\mathrm{I}}_{jk}& =B0+B1\cdot \mathrm{Tem}{\mathrm{p}}_j+B2\cdot {\mathrm{Precip}}_k^{-2}\\ & +B3\cdot (\mathrm{Temp}\cdot \mathrm{Preci}{\mathrm{p}}^{-2}{)}_{jk}+{\epsilon }_{jk}\end{array}$(2) where:

• $iB{B}_i$ – annually damaged relative stem number per stand (0–1)

• SHI – stand hazard index

• SMI – soil moisture index over the growing season

• HTS – host tree share

• Temp_j – mean annual temperature

• Precip⁻²_l – (precipitation)⁻²

• ${\epsilon }_{jk}$ – error term

• B0 to B3 – see Table 2 (3) $vB{B}_i=pB{B}_i\cdot iB{B}_i$(3)

The parameters for the damage probability model and the soil moisture index regression model, as well as their respective standard errors, were retrieved from Tables 1 and 2 in Seidl et al. (2009).

Table 1. Parameter estimates and standard errors of the parameter estimates of the bark beetle probability model (Seidl et al. 2009).

Parameter	Description	Parameter estimate	Standard error
B1	Mean annual temperature (log)	4.234	0.0123
B2	Annual precipitation (log)	−1.196	0.00531
B3	Stand age	−0.0433	0.000142
B4	Stocking density	0.658	0.0249
B5 (dummy 50%)	Host tree share	−1.425	0.0107
B5 (dummy 100%)	Host tree share	−2.606	0.0104
B6	Interaction temp . precip	−0.214	0.00171
B7	Interaction temp . stand age	−0.000383	0.0000239
B8	Interaction temp . density	−0.0527	0.00417
B9 (dummy 50%)	Interaction temp . host tree share	0.0993	0.00183
B9 (dummy 100%)	Interaction temp . host tree share	0.245	0.00178
B10	Interaction precip . stand age	0.00746	0.0000195
B11	Interaction precip . density	−0.104	0.00341
B12 (dummy 50%)	Interaction precip . host tree share	0.278	0.00146
B12 (dummy 100%)	Interaction precip . host tree share	0.457	0.00143
B0	Intercept	−0.301	0.0371

Table 2. Parameter estimates and standard errors of the parameter estimates of the soil moisture index regression model (Seidl et al. 2009).

Parameter	Description	Parameter estimate	Standard error
B1	Mean annual temperature	0.00107	0.0000967
B2	(Mean annual precipitation)^−2	25490	574.4
B3	Interaction temp . precip^−2	8282	58.93
B0	Intercept	−0.0413	0.000943

Seidl et al. (2009) originally fitted their model for Austria using mean annual temperatures (Temp_j), varying from 2 to 15°C, and annual precipitation (Precip_k), varying from 500 to 2000 mm. For stand age (Age_l), the probability model of bark beetle damage was fitted using stands that were 40–120 years old. Stocking density (Density_m) was defined as the basal area of a stand relative to the basal area of a fully stocked stand as given in local yield tables (Seidl et al. 2009) and was originally tested using values from 0.6 to 1.0. Host tree (i.e. P. abies) share (Share_n) was used as a categorical variable with three levels (10%, 50%, and 100% of trees in the stand being P. abies).

Data sources and data analysis

Sensitivity analysis

We ran a sensitivity analysis of pBB using ranges of mean annual temperature and annual precipitation that cover most of the conditions found in Norway (considering the 5th, 50th and 95th percentiles, i.e. 1, 4 and 7°C for annual temperature, and 630, 915 and 1750 mm for annual precipitation). The effect of stand age on pBB was analyzed using the minimum and maximum stand age values found in Norway (40- and 268-year-old). For a sensitivity analysis of the iBB model, we considered the 5th, 50th and 95th precipitation percentiles (630, 915 and 1750 mm) and the full historical range of mean annual temperatures for Norway (from 0.1 to 12°C). For both pBB and iBB sensitivity analyses we also used three levels of host tree share (10%, 50% and 100%).

Current climate and current stand conditions

To obtain a clearer view of the current risk of bark beetle damage in Norway we compounded the annual damage probability $(pB{B}_i)$ and intensity $(iB{B}_i)$ estimations from Seidl’s models to account for the 5-year-period measured in the NNFI data (2015–2019). Damages may happen in any of the 5 years or in none at all. Therefore, we need to calculate the expected value of the combined damage given all possible damage combinations ($2^5=32$). Each distinct damage combination can be represented using binary integers $k\in \{0,1,2,\dots ,31\}$ in the following manner:

• $k=0:00000\to$ no damage

• $k=1:00001\to$ damage in the first period

• $k=2:00010\to$ damage in the second period

• $k=3:00011\to$ damage in the first and second period

• … 

• $k=30:11110\to$ damage in periods 2, 3, 4 and 5

• $k=31:11111\to$ damage in all periods

To calculate the expected value of the combined damage probability given all possible damage combinations $(pB{B}_{5y})$, we first calculated the damage probability for year $i$ $(pB{B}_i)$ (Eq. 1(1) $\begin{array}{rl}pB{B}_i& ={\displaystyle \frac{e^{Z_{jklmn}}}{1+e^{Z_{jklmn}}}}\\ Z_{ijklm}& =B0+B1\cdot \mathrm{LogTem}{\mathrm{p}}_j\\ & +B2\cdot \mathrm{LogPreci}{\mathrm{p}}_k+B3\cdot \mathrm{Ag}{\mathrm{e}}_l\\ & +B4\cdot \mathrm{Densit}{\mathrm{y}}_m+B5\cdot \mathrm{Shar}{\mathrm{e}}_n\\ & +B6\cdot (\mathrm{LogTemp}\cdot \mathrm{LogPrecip}{)}_{jk}\\ & +B7\cdot (\mathrm{LogTemp}\cdot \mathrm{Age}{)}_{jl}\\ & +B8\cdot (\mathrm{LogTemp}\cdot \mathrm{Density}{)}_{jm}\\ & +B9\cdot (\mathrm{LogTemp}\cdot \mathrm{Share}{)}_{jn}\\ & +B10\cdot (\mathrm{LogPrecip}\cdot \mathrm{Age}{)}_{kl}\\ & +B11\cdot (\mathrm{LogPrecip}\cdot \mathrm{Density}{)}_{km}\\ & +B12\cdot (\mathrm{LogPrecip}\cdot \mathrm{Share}{)}_{kn}+{\epsilon }_{jklmn}\end{array}$(1) ). Then we calculated the 31 combinations of damage probabilities (excluding k = 0, which signifies the absence of damage) encoded in the binary representation of $k$ ($pBB(k)$) using Eq. (4(4) $pBB(k)=\prod _{i=1}^{5}bit(i,k)\cdot pB{B}_i+(1-bit(i,k))\cdot (1-pB{B}_i)$(4) ) below. Finally, we summed all the probabilities for the 31 combinations (Eq. 5(5) $pB{B}_{5y}=\sum _{k=1}^{31}pBB(k)$(5) ).

To calculate $pBB(k)$ and $iBB(k)$we defined a helper function $bit(i,k)={\displaystyle \frac{k}{2^i}\text{\%}2}$ that returns the $i$th binary digit of $k$. (4) $pBB(k)=\prod _{i=1}^{5}bit(i,k)\cdot pB{B}_i+(1-bit(i,k))\cdot (1-pB{B}_i)$(4) (5) $pB{B}_{5y}=\sum _{k=1}^{31}pBB(k)$(5) Following a similar approach, we used the damage probability for year $i$ $(iB{B}_i)$ (Eq. 2(2) $\begin{array}{rl}iB{B}_i& ={\displaystyle \frac{1}{1+e^{3.9725-2.9673\cdot \mathrm{SHI}}}}\\ \mathrm{SHI}& =\mathrm{SMI}\cdot \mathrm{HTS}\\ \mathrm{SM}{\mathrm{I}}_{jk}& =B0+B1\cdot \mathrm{Tem}{\mathrm{p}}_j+B2\cdot {\mathrm{Precip}}_k^{-2}\\ & +B3\cdot (\mathrm{Temp}\cdot \mathrm{Preci}{\mathrm{p}}^{-2}{)}_{jk}+{\epsilon }_{jk}\end{array}$(2) ) to calculate the combined intensity of damages $iBB(k)$ in combination $k$ with Eq. (6(6) $iBB(k)=1-\prod _{i=1}^{5}1-bit(i,k)\cdot iB{B}_i$(6) ). (6) $iBB(k)=1-\prod _{i=1}^{5}1-bit(i,k)\cdot iB{B}_i$(6) Damage volume $(vB{B}_{5y})$ was only calculated for stands with a 5-year probability of damage greater than 5%. To calculate 5-year damage volume (%) we multiplied each $pBB(k)$and $iBB(k)$ and then summed all of them (Eq. 7(7) $vB{B}_{5y}(\text{\%})=\sum _{k=0}^{31}pBB(k)\cdot iBB(k)$(7) ). To visualize actual spruce volumes that can be lost due to bark beetles, the damage volume (in m³/ha) was calculated by multiplying the modelled 5-year damage volume (%) by the current NNFI estimates of standing P. abies volume (m³/ha without bark) (Eq. 8(8) $vB{B}_{5y}(m^3/ha)=vB{B}_{5y}(\text{\%})\cdot vol$(8) ). (7) $vB{B}_{5y}(\text{\%})=\sum _{k=0}^{31}pBB(k)\cdot iBB(k)$(7) (8) $vB{B}_{5y}(m^3/ha)=vB{B}_{5y}(\text{\%})\cdot vol$(8) Forest data (stand age, host share, and stocking density) and climatic data (precipitation and temperature) were retrieved from the NNFI plots and the Norwegian Meteorological Institute, respectively, and all variables were measured in the period 2015–2019. While stand age and host share are measured directly in the NNFI plots, relative stocking density is only categorized as either fully stocked or understocked. Therefore, we assumed a relative stocking density of 75% in plots classified as “understocked”, and 100% in plots classified as “fully stocked”.

Calculations of pBB and iBB were performed in RStudio v2022.02.3 (RStudio Team 2022) using R version 4.2.2 (R Core Team (2023)). Probabilities of bark beetle damage, damage volume (%), and damage volume (m³/ha) for each plot were exported to QGIS v3.16.10 and presented on maps.

Time series and climate change conditions

We investigated the effects of future climatic conditions on bark beetle damage by keeping forest variables constant and varying annual temperatures and precipitation according to the projections from different climate models. Stand age was kept as 100 years old, stocking density as 100%, and spruce share as 100%. For the years 1960–2017, we used historical temperature and precipitation data from the Norwegian Meteorological Institute, and from 2018 to 2080 we used model projections. We used the representative concentration pathway (RCP) 4.5 estimated with EURO-CORDEX, a high-resolution regional climate change ensemble that was established for Europe and has been bias-corrected for Norway (Jacob et al. 2014; Wong et al. 2016). We used five EURO-CORDEX simulations resulting from selected combinations of three global climate models (EC-EARTH, CNRM-CM5, and MPI) and three regional climate models (HIRHAM5, RACMO22E, and CCLM4-8-7 – here shortened to CCLM). The five EURO-CORDEX runs used in this study were: CNRM-CM5_CCLM, EC-EARTH_CCLM, EC-EARTH_HIRHAM5, EC-EARTH_RACMO22E, and MPI_CCLM.

Bark beetle damage volume $(vBB)$ under climate change conditions was estimated for all 2924 sample plots each year. We used an annual estimation of damage volume instead of calculating it for a 5-year period, since we wanted to focus on the year-to-year variability in bark beetle damage created by the climatic models.

After having estimated bark beetle damage volume for individual plots, we grouped plots into the seven regions shown in Figure 1 and calculated the average $vBB$ for all plots in each region. Time-series per region were calculated using the averaged $vBB$ values. In the results we focus on the main spruce-growing regions in Norway (Trøndelag, Innlandet and Vestviken), but the results from the other regions are presented in the supplementary material (Figure S11).

Ips typographus development and number of generations in Norway

To evaluate Seidl’s model, we ran a temperature-driven development model to compare if the predicted bark beetle damage coincided with areas where I. typographus is able to complete its development under current climatic conditions. Following the approach of Lange et al. (2006), the development of I. typographus was calculated for the years 2015–2019 using a developmental model parameterized with data from Wermelinger and Seifert (1998). The model uses accumulated daily temperature sums to calculate the time of first beetle flight in spring, the laying of eggs, development of larvae, pupae, and sexual maturation of adults of the first generation. Further, the development of a second generation is modelled, starting with a second flight in summer. The final developmental stage reached at different localities was calculated by running the model until a date after which low temperatures allow no further development (October 31). The median developmental stage reached in the years 2015–2019 is reported for the whole distribution area of I. typographus in Norway and presented at a 1 × 1 km resolution.

Bark beetle monitoring

Trap catches from the Norwegian bark beetle monitoring program for the years 2015–2019 was included as a quantitative foundation for evaluating the modelling results. Every year since 1979, the monitoring program has captured bark beetles in the parts of Norway where I. typographus is present and prevalent, using multifunnel traps and pheromone dispensers releasing methylbutenol, cis-verbenol, and ipsdienol (Økland et al. 2021). Trap counts are a proxy for local population density and is by no means equivalent to bark beetle damage (as the vigor of spruce trees also determines the outcome of bark beetle colonization attempts, especially in an endemic population state). However, trap data does provide a basis for evaluating whether the damage models produce feasible or nonsensical results.

Results

Sensitivity analysis

Our sensitivity analyses showed that temperature had the most impact in the damage probability model (pBB), and it limited the effects of all other variables (Figure 2a). Low temperatures kept the probability of bark beetle damage close to zero even when precipitation and stand age were at their extremes. Stand age had an almost linear effect within the range used to fit Seidl’s model (40–120 years). However, 29% of the 2924 Norwegian plots included in our analyses were older than the fitted range (up to 268 years). In the sensitivity analysis, the probability of bark beetle damage rapidly increased for stands older than 150 years (Figure 2a, and Figure S6 supplementary material), particularly under conditions of elevated temperature and precipitation. This contrasts with empirical evidence suggesting that drier locations are in general more susceptible to bark beetle damage (Wermelinger 2004). After detecting this unrealistic model behavior, all 259 plots (8.8% of the total) older than 150 years were assigned an age of 150 years (see Section 4.1 for a discussion of this).

Figure 2. Sensitivity analysis of mean annual temperature, annual precipitation, stand age, and host share for the bark beetle damage probability model (pBB) (a) and the damage intensity model (iBB) (b). Dashed lines in the upper plot indicate the maximum stand age (120 years-old) used by (Seidl et al. 2009). Stock density is set as 1.

The sensitivity analysis of the damage intensity model revealed a linear increase in damage intensity in response to increasing annual temperature (Figure 2b). An increasing share of P. abies and drier conditions accentuated the impact of temperature. For example, plots with P. abies monocultures, a mean annual temperature of 12°C, and 630 mm annual precipitation are expected to experience 4% annual tree mortality. Under wetter conditions (1750 mm), such plots would experience <2% annual mortality, irrespective of the temperature.

Sixty percent of the NNFI plots included in this study have climate and forest characteristics that fall within the range of the fitting dataset used by Seidl et al. (2009). For annual temperature, 90% of the plots are in the range used for fitting the model (2–15°C). For stand age, 71% of the plots are in the range tested by Seidl et al. (2009) (40–120 years old). Host tree share was a categorical variable (either 10%, 50% or 100%), therefore all NNFI plots were assigned to one of these groups according to their proportion of spruce. For annual precipitation, 92% of the plots fell within the model’s original range (500 to 2000 mm), with the minimum mean annual precipitation being 312 mm and the maximum 4576 mm. Table 3 summarizes the values of current climate and current stand conditions used in the models for all 2924 NNFI plots.

Table 3. Number of sample plots with Norway spruce and average values for stand conditions and climatic variables. Minimum and maximum values are given in parenthesis.

Region	Number of plots	Stand age (years)	Spruce share	Stocking density (0.75; 1)	Annual mean temperature (°C)	Annual precipitation (mm)
Nord-Norge	182	89 (40–150)	0.72 (0.13–0.99)	0.95	3.6 (0.8–5.9)	1722 (443–3080)
Trøndelag	595	101 (40–150)	0.69 (0.11–0.99)	0.96	4.1 (1.3–6.8)	1432 (627–3237)
Innlandet	960	91 (40–150)	0.68 (0.11–0.99)	0.97	3.0 (0.5–6.2)	868 (312–1399)
Østviken	292	84 (40–150)	0.61 (0.11–0.99)	0.98	5.7 (2.5–8.5)	979 (719–1518)
Vestviken	618	96 (40–150)	0.61 (0.12–0.99)	0.97	4.4 (1.0–7.9)	1087 (532–1763)
Sørlandet	171	84 (40–150)	0.54 (0.11–0.99)	0.98	6.4 (2.5–8.3)	1756 (1131–2806)
Vestlandet	107	59 (40–120)	0.71 (0.13–0.99)	0.97	6.3 (3.5–8.6)	2238 (895–4576)
Norway	2924	86 (40–150)	0.65 (0.11–0.99)	0.96	4.8 (0.5–8.6)	1440 (312–4576)

Current (2015–2019) climate and stand conditions

A map of the current 5-year probability of bark beetle damage ($pB{B}_{5y}$), damage volume ($vB{B}_{5y}(\text{\%})$), and damage volume ($vB{B}_{5y}({\mathrm{m}}^3/\mathrm{ha})$), for the regions Vestviken, Innlandet, Trøndelag, and Østviken is shown in Figure 3. Sample plots with less than 5% $pB{B}_{5y}$ (552 blue dots in the map of damage probability in Figure 3) were excluded from the analysis of volume of damage, therefore 1621 plots were used for damage calculations in those regions.

Figure 3. Modelled Ips typographus 5-year damage probability, damage volume (%), and damage volume (m³/ha) in plots from the Norwegian National Forest Inventory (NNFI). Forest variables (stand age, host share, stocking density) and climatic variables (temperature, precipitation) were measured annually from 2015 to 2019. Pie charts show the average shares per damage class for the regions Vestviken (1), Innlandet (2), and Trøndelag (3). Zoomed area in the map shows individual NNFI plots in the Innlandet, Vestviken, Trøndelag, and Østviken regions.

Under current conditions, plots in Vestviken, Innlandet, and Trøndelag have on average an estimated 15% probability of bark beetle damage from 2015 to 2019. Forty-one percent of the plots (891 plots) have between 15 and 80% 5-year probability of bark beetle damage. The current average damage volume is very low (around 0.4%), which translates to a risk of losing on average 0.6 m³/ha of P. abies from 2015 to 2019. Furthermore, 11% of the plots in the three regions are at risk of losing 0.8–3.2% of its P. abies stock, and nine percent of the productive plots are currently at risk of losing 1.7–11 m³/ha over the 5-year period (Figure 3).

Time series and climate change conditions

Bark beetle damage is expected to vary greatly between NNFI plots under climate change (Figures S7, S8, S9, S10, supplementary material), due to differences in climatic conditions and landscape characteristics between regions. Figure 4 shows the estimated annual bark beetle damage volume $(vB{B}_i)$ in Vestviken, Innlandet, and Trøndelag under historical (1960–2017) climatic conditions and under five scenarios from different climate models that simulate climate change conditions (2018–2080). Plots of estimated $vB{B}_i$ for the four remaining regions in Norway are presented in the supplementary material (Figure S11).

Figure 4. Annual modelled damage volume (%) averaged by region for Vestviken, Innlandet, and Trøndelag (from bottom to top). The dotted red line marks the end of the modelling period using historical climate data (2017) and the start of the period using climate data from five EURO-CORDEX simulations with selected combinations of global and regional climate models (see Methods for details).

In Trøndelag, the model estimated an average of 0.06% annual damage volume (%) from 1960 to 2017, with a peak of 0.30% in recent years. Damage intensity increased to an average of 0.11% per year from 2018 to 2080 in this region, with a peak of 0.34% annual damage with the climatic model EC-EARTH_HIRHAM5. Innlandet had an estimated average damage of 0.04% from 1960 to 2017, with occasional small peaks reaching 0.20%. With climate change projections, the estimated average damage volume in Innlandet increased to 0.10% per year, with a peak of 0.40% with the climatic model CNRM-CM5_CCLM. In Vestviken, estimated annual damage for the historical period was 0.09%, with a peak of 0.30%. In this region, climate change projections are estimated to increase the average probability of bark beetle damage to 0.18%, with a peak of 0.70% with climatic variables modeled by EC-EARTH_HIRHAM5.

Ips typographus development and number of generations in Norway

The I. typographus developmental model (Lange et al. 2006) estimated that average temperature conditions in the years 2015–2019 would allow one completed beetle generation in most of the areas estimated by our risk model to have significant bark beetle attack and damage (Figure 5). The potential for bark beetles producing two generations was, in most years, restricted to limited areas in the south of Norway, around the Oslofjord. However, in 2018, which was exceptionally warm, there was potential for two beetle generations over much of SE Norway (Figure S13).

Figure 5. Left: Median developmental potential of local Ips typographus populations in the main Picea abies-growing areas of Norway during the years 2015–2019, modelled using accumulated daily temperature sums (Lange et al. 2006). Right: Average I. typographus trap count per municipality for the years 2015–2019. Trap data is from the Norwegian bark beetle monitoring program. Municipalities with grey color were not sampled during the 5-year period.

Bark beetle monitoring

Data from the Norwegian bark beetle monitoring program for the years 2015–2019 confirmed that I. typographus was present and abundant in the areas where our model estimated at least moderate levels of damage risk (Figure 5).

Discussion

Sensitivity analysis

Most Norwegian National Forest Inventory (NNFI) plots fell within the ranges of climatic and forest conditions used by Seidl et al. (2009) to fit their model. Ninety percent of the NNFI plots were within the mean annual temperature range tested by Seidl et al. (2009) (2–15°C). The remaining 10% of the plots had a mean annual temperature below 2°C. Therefore, no NNFI plots had an annual temperature higher than 15°C (Table 3).

As seen in our sensitivity analysis, the damage intensity model described a clear linear influence of temperature: an increase in annual temperature corresponded directly to an increase in damage intensity levels. Moreover, when conditions were drier, the influence of temperature was more pronounced. For instance, in NNFI plots with Norway spruce monocultures, with an average annual temperature of 12°C, and 630 mm annual precipitation, the model estimated 4% annual tree mortality. In wetter conditions (1750 mm annual precipitation), stands experienced low damage intensities across the whole range of tested temperatures. A noteworthy aspect of the sensitivity analysis is that the minimum damage intensity estimated by the model was 2%, not zero (Figure 2). This might indicate a limitation within the model.

Due to a lack of precise data, we had to assume a stocking density of 75% in understocked NNFI plots and 100% in fully stocked plots. This approach was not fully concordant with the one used by Seidl et al. (2009), as they defined stocking density as “the basal area of a stand relative to the basal area of a fully stocked stand as given in the respective local yield tables”. Our results might have differed if we had been able to use a similar approach as Seidl et al. (2009).

For stand age, Seidl et al. (2009) used a fixed minimum age of 40 years for risk of bark beetle damage when selecting data to calibrate their model. This is a possible limitation of the model, also pointed out by Seidl et al. (2007), since using stand age does not consider that trees in highly productive stands grow more and faster than trees in low productivity stands. A possible alternative approach would be to calibrate the model using a stand-level mean host tree diameter at breast height (DBH). Many studies predicting tree susceptibility to bark beetle damage use a DBH of 15 cm as a minimum threshold for attack (Komonen et al. 2011; Seidl and Rammer 2017). The rationale for this is that the bark of smaller trees is too thin to shelter adult beetles and provide sufficient nutrients for offspring development (Honkaniemi et al. 2018). Also, using actual stand age as an independent variable, as the probability model does, may be problematic because of varying site productivity. Tree vitality, senescence, and resistance to bark beetle attack will vary substantially between high and low productivity areas even for trees of the same age. An alternative to using actual stand age in the model could be to use relative stand age, which, for example, could be the actual stand age divided by the stand’s biologically optimal rotation age, which varies according to site productivity.

The bark beetle damage probability was most influenced by temperature. This is in agreement with several other studies, as temperature is the most influential factor in insect development (Wermelinger 2004; Baier et al. 2007; Seidl et al. 2007). An important limitation of Seidl et al.’s (2009) probability model when applied to Norwegian conditions was the unrealistic behavior of the model in the 259 plots that were older than 150 years (8.8% of the NNFI plots included in this study). Our sensitivity analysis showed that at 4–7°C mean annual temperature the probability of bark beetle damage increases sharply when plots are older than 150 years, especially when precipitation is high (Figure 2a and Figure S6 supplementary material). The effect of precipitation on damage probability shifts from being inversely proportional to stand age in younger stands, to being directly proportional to age in older stands. This is biologically unrealistic, as drier sites tend to have a higher probability of being damaged by bark beetles (Wermelinger 2004). Including the 259 plots with older stands in our model with their original age would produce biologically unrealistic results. Excluding those plots completely, however, would result in a loss of important data on older stands, since the probability of bark beetle damage increases with stand age (Wermelinger 2004). Therefore, as a compromise, we included the 259 plots but assigned all of them with an age of 150 years.

Bark beetle risk in Norway

Under current climate and forest conditions, i.e. considering the 5-years from 2015 to 2019, the estimated damage volume (%) due to bark beetle infestations was around 0.4%. This translates to an average potential loss of 0.6 m³ of Norway spruce per hectare in the NNFI plots. Across the regions Trøndelag, Innlandet, and Vestviken, about 11% of the plots face a risk of losing 0.8–3.2% of its Norway spruce stock due to bark beetle attacks. Furthermore, nine percent of the productive plots are at risk of losing from 1.7 to 11 m³ of spruce per hectare (Figure 3). For comparison, the growth of Norwegian spruce stands in locations with medium site productivities typically varies from 35 to 50 m³/ha over a 5-year span, while standing volumes for 50- to 70-year-old stands at such sites may vary from 100 to 350 m³/ha depending on the management.

Even though the model we used includes important forest variables (stand age, host tree share, and stocking density) and climatic variables (temperature and precipitation), the estimated bark beetle damage volumes still appear to be conservative. This may be because Seidl’s models do not consider other factors that can predispose stands to a greater damage probability and damage intensity, such as wind damage and droughts, bedrock variability, and soil composition. For example, bedrock and soil properties in parts of Vestviken – especially Vestfold – might predispose stands for drought stress (Worrell 1983). This may, in turn, make the trees more susceptible to bark beetle attacks and could partially explain the high spruce mortality in the region in the 1970s.

In an early outbreak phase, bark beetles need weakened trees for reproduction (Lausch et al. 2013) and outbreaks are therefore often associated with past forest disturbances (Annila and Petäistö 1978; Berec et al. 2013; Marini et al. 2017; Fernandez-Carrillo et al. 2020). Most critically, I. typographus developmental rates may be substantially influenced by even minor increases in mean temperatures, which agrees well with what was found for damage volume in the present study (Figure 4). Higher temperatures due to climate change will, therefore, increase the risk of significant timber losses from bark beetle attacks. For example, Norway and Sweden could previously support only one bark beetle generation per year, but now some regions support the development of a second generation because of higher temperatures (Jönsson et al. 2009; Økland et al. 2021).

Under climate change conditions, the annual damage volume is estimated to reach 0.70% in Vestviken, 0.40% in Innlandet, and 0.34% in Trøndelag. The average annual damage volume in those regions will likely double (from 0.06% to 0.13%). These damage volumes are not very large but align well with the results of Zimová et al. (2020) for central-eastern Slovakia, where climate change is predicted to double annual bark beetle-induced damages from 2.2 to 4.5 m³/ha.

Seidl’s model considers each year to be independent from the previous ones, but this is not necessarily a realistic assumption. In practice, if forest managers do not act promptly to e.g. remove weakened trees after a disturbance event, such as a heatwave and/or severe drought, affected trees could be susceptible to bark beetle attack also in the subsequent years, even if the weather conditions return to normal. Ecological disturbances might thus have longer-term consequences for bark beetle damage probabilities and intensities.

Model evaluation

The magnitude of the damage probabilities reported by Seidl’s model (Figure 3) is, in our experience with Norwegian forestry, realistic. In limited areas, particularly in Vestfold county in the Vestviken region, spruce forests are vulnerable to drought stress due to geological conditions and soil types. In these areas, even higher damage volumes than those estimated by the model may be reasonable. High incidences of attack and damage may also be more widespread in areas where local topographical conditions cause high runoff and very dry soil conditions. In fact, due to the high risk of drought and bark beetle damage in such localities, many foresters are moving away from planting Norway spruce in such drought-prone localities. There are no quantitative data on bark beetle damage rates in Norway, but the geographical distribution of damage reported by Seidl’s model is in line with our experience with bark beetle attack rates in the country. However, in the absence of quantitative data this assessment is based solely on our long-term dialog with stakeholders in Norwegian forestry.

While we do not have access to data on damage or attack rates in Norwegian forests, our estimation of the developmental potential of I. typographus shows that areas with moderate to high bark beetle risks largely overlap with areas where temperature sums are predicted to allow the beetles to complete one generation per year (Figure 5). Bivoltinism has been observed in areas around the Oslofjord and this aligns with the predictions of the developmental model, and with the distribution of higher damage predicted by Seidl’s model. While the actual prevalence of bivoltinism is not well known in these areas, it may currently be such a rare and geographically limited phenomenon that the implications for damage rates are small. The higher rates of damage in this area, as predicted by our model and reported by forestry stakeholders, are more likely driven by high temperatures (by Norwegian standards) and the relatively dry climate in the south-eastern part of Norway. Data from the Norwegian bark beetle monitoring program confirms that I. typographus is abundant in the areas where significant risk was predicted by our model (Figure 5). As with the developmental predictions, trap catches also align well with model predictions.

Alternative models and future studies

We see some alternative approaches to build a bark beetle damage model tailored to Norwegian conditions. A more complex model that has been widely used in European conditions for smaller areas is the bark beetle module from the forest landscape model Iland (Seidl and Rammer 2017; Dobor et al. 2020). Unlike the model from Seidl et al. (2009), Iland incorporates the interaction between wind and bark beetle damage. This consideration is vital for more accurate damage prediction, as many bark beetle outbreaks are frequently linked to windthrow events (Wermelinger 2004). However, apart from requiring spatially explicit data, the Iland model needs some data that are not available in Norway, such as recent bark beetle damage records and a verified relationship between temperature and beetle development under Norwegian conditions. Another option for obtaining a model tailored to Norwegian conditions is to build a new model from scratch. This also, however, requires data that is not easily produced, such as bark beetle infestation rates or number of trees killed by the insect. Furthermore, bark beetle damage records are needed to validate the Seidl et al. (2009) model or any other models. Therefore, even with the limitations the Seidl et al. (2009) model presents under Norwegian forest and climate conditions, we believe the results presented here provide a useful first estimate of probability of I. typographus damage, damage intensity, and risk in Norway.

Adaptative management strategies are crucial to avoid major economic losses to the Norwegian forest sector. Some preliminary steps have been taken by the Norwegian Agriculture Agency (Landbruksdirektoratet 2021) in a preparedness plan for how to deal with challenges related to bark beetles. The plan is inspired by the review by Hlásny et al. (2021a), which highlighted the need for a comprehensive management framework with an interdisciplinary social-ecological approach. The Norwegian plan deals with roles, responsibilities and expectations of different authorities and stakeholders, as well as operative measures. Most operative measures are relatively short-term, for situations in which the probability of an outbreak is high, or the outbreak has already started. Such measures include identification and prioritization of potential areas for harvest (areas susceptive to drought, windthrow, and vulnerable stand borders), preparedness measures for harvesters, and identification and removal of already attacked trees. Some long-term measures related to silvicultural practices are also mentioned in the preparedness plan: e.g. thinning to reduce competition between trees and increase vitality, diversifying tree species in regeneration, and facilitate a more diverse stand/forest structure.

However, Norway has seen limited efforts and considerations on long-term measures against bark beetles in a climate change perspective, both by authorities and researchers. Furthermore, research and public debate on economic losses, profitability, and adaptations of forest management strategies is almost absent. Therefore, an important and timely topic for future research is to use forest growth simulators to do long-term risk analyses of bark beetle outbreaks in Norway under different forest conditions. The modelled damage estimates presented in this paper would be a good starting point for such analyses. Such studies could provide more realistic estimates of bark beetle risks in Norway and bring new insight on how to adapt forest management in the country to face future disturbances. Additionally, establishing a systematic approach to collect data on the volume of trees damaged by beetles is essential for refining these models and enhancing their accuracy in Norway.

Conclusions

Most of the NNFI plots included in this study fall within the range of climatic and forest conditions used by Seidl et al. (2009) to fit their model. We believe the results presented here provide useful initial estimates for probability and volume of I. typographus damage throughout Norway. However, a future improved model adapted to Norwegian conditions should include data from older stands and lower temperatures. Unfortunately, the possibility to validate a model tailored to Norwegian conditions is hindered by the lack of high-resolution data on infestation rates and number of trees killed by I. typographus. This lack of data is partly a consequence of the limited damage caused by I. typographus in Norway over the last 40 years.

The most important spruce-growing regions in Norway (Trøndelag, Innlandet, and Vestviken) currently appear to have a low risk of experiencing severe bark beetle damage. However, climate change is expected to double the annual damage volume caused by I. typographus. Therefore, forest managers should incorporate projections about future bark beetle damages as described in this study and adapt their short- and long-term management methods accordingly.

Acknowledgements

All authors contributed to the study conception and design. Data analyses were performed by JMNR, JG, and CAF. The manuscript was written by JMNR, JG, PK, and TE. The manuscript has been read and approved by all authors.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

Funding for this study was provided by the Norwegian Research Council through the Climate Smart Forestry Norway project [NFR 302701], the Spruce Beetle project [NFR 325319], and the Escalera Excelencia Internacional Junta de Castilla y León (CL-EI-2021-05) relativo al (CLU-2019-01).