Simulating Ips typographus L. outbreak dynamics and their influence on carbon balance estimates with ORCHIDEE r8627

New (a)biotic conditions resulting from climate change are expected to change disturbance dynamics, such as windthrow, forest fires, droughts, and insect outbreaks, and their interactions. These unprecedented natural disturbance dynamics might alter the capability of forest ecosystems to buffer atmospheric CO₂ increases, potentially leading forests to transform from sinks into sources of CO₂. This study aims to enhance the ORCHIDEE land surface model to study the impacts of climate change on the dynamics of the bark beetle, Ips typographus, and subsequent effects on forest functioning. The Ips typographus outbreak model is inspired by previous work from Temperli et al. (2013) for the LandClim landscape model. The new implementation of this model in ORCHIDEE r8627 accounts for key differences between ORCHIDEE and LandClim: (1) the coarser spatial resolution of ORCHIDEE; (2) the higher temporal resolution of ORCHIDEE; and (3) the pre-existing process representation of windthrow, drought, and forest structure in ORCHIDEE. Simulation experiments demonstrated the capability of ORCHIDEE to simulate a variety of post-disturbance forest dynamics observed in empirical studies. Through an array of simulation experiments across various climatic conditions and windthrow intensities, the model was tested for its sensitivity to climate, initial disturbance, and selected parameter values. The results of these tests indicated that with a single set of parameters, ORCHIDEE outputs spanned the range of observed dynamics. Additional tests highlighted the substantial impact of incorporating Ips typographus outbreaks on carbon dynamics. Notably, the study revealed that modeling abrupt mortality events as opposed to a continuous mortality framework provides new insights into the short-term carbon sequestration potential of forests under disturbance regimes by showing that the continuous mortality framework tends to overestimate the carbon sink capacity of forests in the 20- to 50-year range in ecosystems under high disturbance pressure compared to scenarios with abrupt mortality events. This model enhancement underscores the critical need to include disturbance dynamics in land surface models to refine predictions of forest carbon dynamics in a changing climate.

1 Introduction

Future climate will likely bring new abiotic constraints through the co-occurrence of multiple connected hazards, e.g., “hotter droughts”, which are droughts combined with not only heat waves (Allen et al., 2015; Zscheischler et al., 2018), but also new biotic conditions from interacting natural and anthropogenic disturbances, such as insect outbreaks following windthrow or forest fires (Seidl et al., 2017). Unprecedented natural disturbance dynamics might alter biogeochemical cycles – specifically the capability of forest ecosystems to buffer the CO₂ increase in the atmosphere (Hicke et al., 2012; Seidl et al., 2014) and the risk that forests are transformed from sinks into sources of CO₂ (Kurz et al., 2008a). The magnitude of such alteration, however, remains uncertain principally due to the lack of impact studies that include disturbance regime shifts on a global scale (Seidl et al., 2011).

Land surface models are used to study the relationships between climate change and the biogeochemical cycles of carbon, water, and nitrogen (Ciais et al., 2005; Cox et al., 2000; Friedlingstein et al., 2006; Luyssaert et al., 2018; Zaehle and Dalmonech, 2011). Many of these models use background mortality to obtain an equilibrium in their biomass pools. This classic approach towards forest dynamics, which assumes steady-state conditions over long periods of time, may not be suitable for assessing the impacts of disturbances on shorter timescales under a fast-changing climate. This could be considered a shortcoming in land surface models because disturbances can have significant impacts on ecosystem services, such as water regulation, carbon sequestration, and biodiversity (Quillet et al., 2010). Mechanistic approaches that account for a variety of mortality causes, such as age, size, competition, climate, and disturbances, are now being considered and tested to simulate forest dynamics more accurately (Migliavacca et al., 2021). For example, the land surface model ORCHIDEE accounts for mortality from interspecific competition for light in addition to background mortality (Naudts et al., 2015). Implementing a more mechanistic view on mortality is thought to be essential for improving our understanding of the impacts of climate change on forest dynamics and the provision of ecosystem services.

Land surface models also face the challenge of better describing mortality, particularly when it comes to ecosystem responses to “cascading disturbances”, where legacy effects from one disturbance affect the next (Buma, 2015; Zscheischler et al., 2018). Biotic disturbances, such as bark beetle outbreaks, strongly depend on previous disturbances as their infestation capabilities are higher when tree vitality is low, for example, following drought or storm events (Seidl et al., 2018). This illustrates how interactions between biotic and abiotic disturbances can have substantial effects on ecosystem dynamics and must be accounted for in land surface models to improve our understanding of the impacts of climate change on forest dynamics (Seidl et al., 2011; Temperli et al., 2013). While progress has been made towards including abrupt mortality from individual disturbance types such as wildfire (Lasslop et al., 2014; Migliavacca et al., 2013; Yue et al., 2014), windthrow (Chen et al., 2018), and drought (Yao et al., 2022), the interaction of biotic and abiotic disturbances remains both a knowledge and modeling gap (Kautz et al., 2018).

Bark beetle infestations are increasingly recognized as disturbance events of regional to global importance (Bentz et al., 2010; Kurz et al., 2008a; Seidl et al., 2018). Notably, a bark beetle outbreak ravaged over 90 % of Engelmann spruce trees across approximately 325 000 ha in the Canadian and American Rocky Mountains between 2005 and 2017 (Andrus et al., 2020). In Europe, the bark beetle Ips typographus has been involved in up to 8 % of total tree mortality due to natural disturbances from 1850 to 2000 (Hlásny et al., 2021a). A recent increase in beetle activity, particularly following mild winters (Andrus et al., 2020; Kurz et al., 2008b), windthrow (Mezei et al., 2017), and droughts (Nardi et al., 2023) has been well documented (Hlásny et al., 2021a; Pasztor et al., 2014), underscoring the need to integrate bark beetle dynamics into land surface modeling.

Past studies used a variety of approaches to model the impacts of bark beetles on forests. While some models treated bark beetle outbreaks as background mortality (Luyssaert et al., 2018; Naudts et al., 2016), others dynamically modeled these outbreaks within ecosystems (Jönsson et al., 2012; Seidl and Rammer, 2016; Temperli et al., 2013). Studies with prescribed beetle outbreaks tend to focus on the direct effects of the outbreak on forest conditions and carbon fluxes but are likely to overlook more complex feedback processes, such as interactions with other disturbances and longer-term impacts. Conversely, dynamic modeling of beetle outbreaks, provides a more comprehensive view by incorporating the life cycle of bark beetles and tree defense mechanisms and ensuing alterations in forest composition and functionality.

Simulation experiments for Ips typographus outbreaks using the LPJ-GUESS vegetation model highlighted regional variations in outbreak frequencies, pinpointing climate change as a key exacerbating factor (Jönsson et al., 2012). Simulation experiments with the iLand landscape model suggested that almost 65 % of Ips typographus outbreaks are aggravated by other environmental drivers (Seidl and Rammer, 2016). A 4 °C temperature increase could result in a 265 % increase in disturbed area and a 1800 % growth in the average patch size of the disturbance (Seidl and Rammer, 2016). Disturbance interactions were 10 times more sensitive to temperature changes, boosting the disturbance regime's climate sensitivity. The results of these studies justify the inclusion of interacting disturbances in land surface models, such as ORCHIDEE, which are used in future climate predictions and impact studies (Boucher et al., 2020).

The objectives of this study are to (1) develop and implement a spatially implicit outbreak model for Ips typographus in the land surface model ORCHIDEE inspired by the work from Temperli et al. (2013) and (2) use simulation experiments to characterize the behavior of this newly added model functionality.

2 Model description

2.1 The land surface model ORCHIDEE

ORCHIDEE is the land surface model of the IPSL (Institut Pierre-Simon Laplace) Earth system model (Boucher et al., 2020; Krinner et al., 2005). ORCHIDEE can, however, also be run uncoupled as a standalone land surface model forced by temperature, humidity, pressure, precipitation, and wind fields. Unlike the coupled setup, which needs to run on the global scale, the standalone configuration can cover any area ranging from a single grid point to the global domain. In this study, ORCHIDEE was run as a standalone land surface model.

ORCHIDEE does not enforce any particular spatial resolution. The spatial resolution is an implicit user setting that is determined by the resolution of the climate forcing (or the resolution of the atmospheric model in a coupled configuration). ORCHIDEE can run on any temporal resolution. This apparent flexibility is somewhat restricted as processes are formalized at given time steps: half-hourly (e.g., photosynthesis and energy budget), daily (i.e., net primary production), and annual (i.e., vegetation demographic processes). With the current model architecture, meaningful simulations should have a temporal resolution of 1 min to 1 h for the calculation of energy balance, water balance, and photosynthesis.

ORCHIDEE utilizes metaclasses to discretize the global diversity in vegetation. The model includes 13 metaclasses by default, including one class for bare soil, eight classes for various combinations of leaf-type and climate zones of forests, two classes for grasslands, and two classes for croplands. Each metaclass can be further subdivided into an unlimited number of plant functional types (PFTs). The current default setting of ORCHIDEE distinguishes 15 PFTs where the metaclass of C₃ grasslands have been separated into boreal, temperate, and tropical C₃ grassland PFTs.

At the beginning of a simulation, each forest PFT in ORCHIDEE contains a monospecific forest stand that is structured by a user-defined but fixed number of diameter classes (three by default). Throughout the simulation, the boundaries of the diameter classes are adjusted to accommodate changes in the stand structure, while the number of classes remains constant. Flexible class boundaries provide a computationally efficient approach to simulate different forest structures. For instance, an even-aged forest is simulated using a small diameter range between the smallest and largest trees, resulting in all trees belonging to the same stratum. Conversely, an uneven-aged forest is simulated by applying a wide range between diameter classes so that different classes represent different canopy strata.

The model uses allometric relationships to link tree height and crown diameter to stem diameter. Individual tree canopies are not explicitly represented; instead, a canopy structure model based on simple geometric forms (Haverd et al., 2012) has been included in ORCHIDEE (Naudts et al., 2015). Diameter classes represent trees with different mean diameter and height, which informs the user about the social position of trees within the canopy. Intra-stand competition is based on the basal area of individual trees, which accounts for the fact that trees with a higher basal area occupy dominant positions in the canopy and are therefore more likely to intercept light and thus contribute more to stand-level photosynthesis and biomass growth compared to suppressed trees (Deleuze et al., 2004). If recruitment occurs, diameter classes evolve into cohorts. However, in the absence of recruitment, all diameter classes contain trees of the same age.

Individual tree mortality from self-thinning, wind storms, and forest management is explicitly simulated. Other sources of mortality are implicitly accounted for through a so-called constant background mortality rate. Furthermore, age classes (four by default) can be used after land cover change, forest management, and disturbance events to explicitly simulate the regrowth of the forest. Following a land cover change, biomass and soil carbon pools (but not soil water columns) are either merged or split to represent the various outcomes of a land cover change. The ability of ORCHIDEE to simulate dynamic canopy structures (Chen et al., 2016; Naudts et al., 2015; Ryder et al., 2016), a feature essential to simulate both the biogeochemical and the biophysical effects of natural and anthropogenic disturbances, is exploited in other parts of the model, i.e., precipitation interception, transpiration, energy budget calculations, radiation scheme, and calculation of the absorbed light for photosynthesis.

Since revision 7791, mortality from Ips typographus outbreaks is explicitly accounted for and thus conceptually excluded from the so-called environmental background mortality. Subsequently, changes in canopy structure resulting from growth, forest management, land cover changes, wind storms, and Ips typographus outbreaks are accounted for in the calculations of the carbon, water, and energy exchanges between the land surface and the atmosphere. For details on the functionality of the ORCHIDEE model that is not of direct relevance for this study, e.g., the energy budget calculations, soil hydrology, snow phenology, albedo, roughness, photosynthesis, respiration, phenology, carbon and nitrogen allocation, land cover changes, product use, and nitrogen cycle, readers are referred to Krinner et al. (2005), Zaehle and Friend (2010), Naudts et al. (2015), Vuichard et al. (2019).

2.2 Origin of the bark beetle (Ips typographus) model – the LandClim legacy

Although mortality from windthrow (Chen et al., 2018) and forest management (Luyssaert et al., 2018; Naudts et al., 2016) was already accounted for in ORCHIDEE prior to r8627, insect outbreaks and their interaction with other disturbances were not. The LandClim model (Schumacher, 2004) and, more specifically, the Ips typographus model developed by Temperli et al. (2013) have been used as a basis to develop the Ips typographus model in ORCHIDEE r8627.

LandClim is a spatially explicit stochastic landscape model in which forest dynamics are simulated at a yearly time step for 10–100 km² landscapes consisting of 25 m × 25 m patches. Within a patch recruitment, growth, mortality, and competition among age cohorts of different tree species are simulated with a gap model (Bugmann, 1996) in response to monthly mean temperature, climatic drought, and light availability. LandClim, for which a detailed description can be found in Schumacher (2004), Temperli et al. (2013), includes the functionality to simulate the decadal dynamics and consequences of Ips typographus outbreaks at the landscape scale (Temperli et al., 2013). In the LandClim approach, the extent, occurrence and severity of beetle-induced tree mortality are driven by the landscape susceptibility, beetle pressure, and infested tree biomass. While the LandClim beetle model was designed and structured to be generally applicable for Northern Hemisphere climate-sensitive bark beetle–host systems, it was parameterized to represent disturbances by the Ips typographus in Norway spruce (Picea abies Karst.; Temperli et al., 2013).

As LandClim and ORCHIDEE are developed for different purposes, their temporal and spatial scales differ. These differences in model resolution justify developing a new model while following the principles embedded in the LandClim approach. LandClim assesses bark beetle damage at 25 m × 25 m patches, and, to do so, it uses information from other nearby patches as well as landscape characteristics, such as slope, aspect, and altitude. The susceptibility of a landscape to bark beetle infestations is calculated using multiple factors, such as drought-induced tree resistance, age of the oldest spruce cohort, proportion of spruce in the patch's basal area, and spruce biomass damaged by windthrow. These drivers are presented as sigmoidal relationships ranging from 0 to 1 (denoting none to maximum susceptibility, respectively) that are combined in a susceptibility index for each Norway spruce cohort in a patch. Bark beetle pressure is quantified as the potential number of beetles that can infest a patch, and its calculation considers, among others, previous beetle activity, maximum possible spruce biomass that beetles could kill, and temperature-dependent bark beetle phenology. Finally, the susceptibility index and beetle pressure are used to estimate the total infested tree biomass and total biomass killed by bark beetles for each cohort within a patch.

In ORCHIDEE, however, the simulation unit is about 6 orders of magnitude larger, i.e., 25 km × 25 km. Hence, a single grid cell in ORCHIDEE exceeds the size of an entire landscape in LandClim. Where landscape characteristics in LandClim can be represented by statistical distributions, the same characteristics in ORCHIDEE are represented by single values. These differences between LandClim and ORCHIDEE imply that the original bark beetle model cannot be implemented in ORCHIDEE without substantial adjustments; the model at the ORCHIDEE grid cell should work without requiring spatial information and statistical distributions of landscape characteristics because those are not available in ORCHIDEE.

2.3 Bark beetle outbreak development stages

Bark beetle outbreak development stages (Fig. 1) are useful to understand the dynamics of an outbreak (Edburg et al., 2012; Hlásny et al., 2021a; Wermelinger, 2004). Nonetheless, the outbreak model in ORCHIDEE r8627 simulates the dynamic of the Ips typographus outbreak as a continuous process. Hence, endemic, epidemic, buildup and post-epidemic stages are not explicitly simulated. In this study, outbreak development stages were only introduced to structure the model description. If needed, these stages could be distinguished while post-processing the simulation results if (arbitrary) thresholds are set for specific variables such as DR_beetles.

Table 1 List of symbols.

Figure 1 Dynamic interplay of the different host and beetle characteristics during a bark beetle (Ips typographus) outbreak. The time window spans four outbreak development stages: buildup, epidemic, post-epidemic, and endemic. The curves represent key characteristics, showing the growth in beetle population and subsequent decline in host population. I_hosts dead characterizes the presence of defenseless uprooted or cut spruce trees, i_{hosts  alive},characterizes living spruce trees that could become hosts for the bark beetles, i_{hosts susceptibility} susceptibility of spruce trees to bark beetle attacks, i_{beetles mass attack} quantifies the capability of the bark beetles to mass attack, and i_{beetles  survival} characterizes the survival of bark beetles. Host and bark beetle characteristics are detailed in the subsequent text. When the density of the host trees is declining due to an increased host mortality from the bark beetle outbreak itself, the competition between trees for light and nutrients declines as well. As a consequence, the host susceptibility decreases, which in ORCHIDEE is the main pathway for an outbreak to move back to the endemic phase. After 1 year, the wood from a storm is not fresh enough for bark beetles to breed in, so i_{host  dead} goes to zero. In ORCHIDEE, the bark beetle population needs to be capable of mass attacking living trees within a year to make the transition from the buildup to the epidemic phase.

2.4 Bark beetle (Ips typographus) damage in ORCHIDEE

The biomass of trees killed by bark beetles in 1 year and one grid cell (B_{beetles kill}) is calculated as the product of the biomass of trees attacked by bark beetles (B_{beetles attacked}) by the probability of a successful attack (P_success, a) averaged over the number of spruce age classes and weighted by their actual fraction ($F_{\text{spruce,}a}/F_{\text{spruce}}$). The approach assumes that a successful beetle colonization always results in the death of the attacked tree, which is a simplification of reality (Schlyter et al., 1989).

$$\begin{array}{}\text{(1)}& B_{\text{beetles kill}}={\sum }_{a=\mathrm{1}}^{\mathrm{3}}{P}_{\text{success,}a}\times B_{\text{beetles attacked}}\times {\displaystyle \frac{F_{\text{spruce,}a}}{F_{\text{spruce}}}}\end{array}$$

During the endemic stage, B_{beetles attacked} and B_{beetles kill} are at their lowest values, and the damage from bark beetles has little impact on the structure and function of the forest. Losses from B_{beetles kill} can be considered contributors to the background mortality.

The biomass of trees attacked by bark beetles (B_{beetles attacked}) is the outcome of bark beetles that successfully overcame the tree defenses and succeeded in boring holes in the bark in order to reach the sapwood. B_{beetles attacked} is calculated at the grid cell by multiplying the actual stand biomass of spruce (B_total) and the probability that bark beetles attack spruce trees in the grid cell (P_attacked).

$$\begin{array}{}\text{(2)}& B_{\text{beetles attacked}}=B_{\text{total}}\times P_{\text{attacked}}\end{array}$$

P_attacked represents the ability of the bark beetles to spread and to locate new suitable spruce trees as hosts for breeding. P_attacked is calculated by the product of two indexes (all indexes in this study are denoted i and are analogue to the susceptibility indexes from Temperli et al., 2013): (1) the beetle pressure index (i_{beetles  pressure}), which is a proxy for the bark beetle population, and (2) the stand attractiveness index (i_{hosts attractiveness}) is related to its health and reflects the ability of the forest to resist an external stressor such as bark beetle attacks.

$$\begin{array}{}\text{(3)}& P_{\text{attacked}}=i_{\text{hosts attractiveness}}\times i_{\text{beetles\hspace{0.17em} pressure}}\end{array}$$

2.5 Host attractiveness

The stand attractiveness index (i_{hosts attractiveness}) represents how interesting a stand is for a new bark beetle colony. When i_{hosts attractiveness} tends to 0, the stand is constituted mainly by healthy trees, which are less attractive for beetles, whereas an i_{hosts attractiveness} value approaching 1 represents a highly stressed spruce stand suitable for colonization by bark beetles. Factors that contribute to the stress of a forest are nitrogen limitation, limited carbohydrate reserves, and monospecific spruce forest. Trees experiencing extended periods of environmental stress are expected to have lower carbon and nitrogen reserves available for defense compounds, making them vulnerable to bark beetle attacks even at relatively low beetle population densities (Raffa et al., 2008). Nonetheless, reserve pools in ORCHIDEE r8627 have not yet been evaluated, so, instead, proxies were used, such as long-term drought (PWS_max) and relative density index (i_{rd,  spruce}), which were already simulated in ORCHIDEE r8627.

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{i}_{\text{hosts attractiveness}}& =max(i_{\text{hosts\hspace{0.17em} competition}},i_{\text{hosts defense}})\\ & \times i_{\text{hosts share}},\end{array}\end{array}$$

where i_{hosts  competition} and i_{hosts defense} both represent proxies for the reduction in the nitrogen and carbohydrate reserve due to strong competition for light and soil resources and consecutive years that are drier than average. For this study, the maximum drought intensity during the last 3 years (PWS_max) is considered as a proxy for spruce stand healthiness:

$$\begin{array}{}\text{(5a)}& {\displaystyle }{i}_{\text{hosts defense}}=\mathrm{1}/(\mathrm{1}+e^{S_{\text{drought}}\cdot (\mathrm{1}-{\text{PWS}}_{\text{max}}-{\text{PWS}}_{\text{limit}})})\text{(5b)}& {\displaystyle }\begin{array}{rl}{\text{PWS}}_{\text{max}}& ={\sum }_{a=\mathrm{1}}^{\mathrm{3}}max({\text{PWS}}_y,{\text{PWS}}_{y-\mathrm{1}},{\text{PWS}}_{y-\mathrm{2}}{)}_{\text{spruce}},\\ & \times {\displaystyle \frac{F_{\text{spruce,}a}}{F_{\text{spruce}}}},\end{array}\end{array}$$

where PWS_y, spruce is the average daily plant water stress index over the growing season for the spruce stand and is equal to 0 when plants are highly stressed and equal to 1 when plants are not stressed. PWS_limit is the plant water stress below which the healthiness of the stand will strongly be affected. The number of age class (a) within the stand is equal to 3 in this study. In addition to drought, an overstocked forest may also decrease the overall healthiness of a spruce stand (i_{hosts  competition}).

$$\begin{array}{}\text{(6a)}& i_{\text{hosts\hspace{0.17em} competition}}=\mathrm{1}/(\mathrm{1}+e^{S_{\text{competition}}\cdot (i_{\text{rd\hspace{0.17em} spruce}}-i_{\text{rd\hspace{0.17em} limit}})})\end{array}$$

In ORCHIDEE, the relative density index (i_rd) is used to quantify the competition between trees at the stand level. At an i_rd of 1, the forest is expected to be at its maximum density given the carrying capacity of the site, implying the highest level of competition between trees. i_{rd  limit} represents the limit at which the bark beetle outbreak starts to decline because of a lack of suitable host trees. The severity of bark-beetle-caused tree mortality decreases when we increase the spatial resolution from the stand to the landscape scale. At the landscape scale, which can cover areas up to 2500 km², the duration of mortality may be longer and the severity lower because beetles disperse across the landscape and cause mortality at different times. This distinction is important for interpreting model results, particularly when considering parameters like i_{rd  limit} in the ORCHIDEE model. The parameter i_{rd  limit} describes the proportion of trees surviving after an outbreak and should therefore be adjusted for the spatial scale of a grid cell in ORCHIDEE. In a model setup, where a grid cell represents a single stand (∼1 ha), i_{rd  limit} should be close to 0, indicating that nearly all trees may be killed. However, in a simulation with grid cells representing 2500 km², not all trees will be killed, which is reflected in setting i_{rd  limit} to 0.4.

The value of i_{rd  spruce} is computed as follows:

$$\begin{array}{}\text{(6b)}& i_{\text{rd\hspace{0.17em} spruce}}={\sum }_{a=\mathrm{1}}^{\mathrm{3}}{\displaystyle \frac{D_{\text{spruce,}a}}{D_{\text{max}}}}\times {\displaystyle \frac{F_{\text{spruce,}a}}{F_{\text{spruce}}}},\end{array}$$

where D_spruce, a is the current tree density of an age class and F_spruce, a is the fraction of spruce in the grid cell that resides in this age class. D_max represents the maximum stand density of a stand given its diameter. In ORCHIDEE, D_max is calculated based on the quadratic mean diameter (cm) and two species specific parameters, α and β:

$$\begin{array}{}\text{(6c)}& D_{\text{max}}=({\text{dia}}_{\text{quadratic}}/\mathit{\alpha }{)}^{(\mathrm{1}/\mathit{\beta })}.\end{array}$$

The index i_hosts share (used in Eq. 4) takes into account that in a mixed tree species landscape, even a few non-host trees may chemically hinder bark beetles in finding their host trees (Zhang and Schlyter, 2004), explaining why insect pests, including Ips typographus outbreaks, often cause more damage in pure compared to mixed stands (Nardi et al., 2023). ORCHIDEE r8627 does not simulate multi-species stands but does account for landscape-level heterogeneity of forests with different plant functional types. The bark beetle model in ORCHIDEE assumes that within a grid cell, the fraction of spruce over other tree species is a proxy for the degree of mixture:

$$\begin{array}{}\text{(7a)}& i_{\text{hosts share}}=\mathrm{1}/(\mathrm{1}+e^{S_{\text{share}}\cdot ({\text{SH}}_{\text{spruce}}-{\text{SH}}_{\text{limit}})}),\end{array}$$

where

$$\begin{array}{}\text{(7b)}& {\text{Sh}}_{\text{spruce}}=F_{\text{non-spruce}}/F_{\text{spruce}}.\end{array}$$

2.6 Implicit representation of bark beetle populations

The bark beetle pressure index (i_{beetles  pressure}) is now formulated based on two components: (1) the bark beetle breeding index of the current year (i_{beetles  generation}) and (2) an index of the loss of tree biomass in the previous year due to a bark beetle infestation (i_{beetles  activity}). The index i_{beetles  activity} is thus a proxy for the previous year's bark beetle activity. The expression accounts for the legacy effect of bark beetle activities by averaging activities over the current and previous years. In this approach, the susceptibility index (i_{beetles  survival}) serves as an indicator for increased bark beetle survival, which could result from favorable conditions for beetle demography (see the next section).

$$\begin{array}{}\text{(8)}& \begin{array}{rl}{i}_{\text{beetles\hspace{0.17em} pressure}}& =i_{\text{beetles\hspace{0.17em} survival}}\\ & \times {\displaystyle \frac{(i_{\text{beetles\hspace{0.17em} generation}}+i_{\text{beetles\hspace{0.17em} activity}})}{\mathrm{2}}}\end{array}\end{array}$$

The model calculates i_{beetles  generation} from a logistic function, which depends on the number of generations a bark beetle population can sustain within a single year:

$$\begin{array}{}\text{(9)}& i_{\text{beetles\hspace{0.17em} generation}}=\mathrm{1}/\left(\mathrm{1}+e^{-S_{\text{generation}}\cdot (\frac{{\text{DD}}_{\text{eff}}}{{\text{DD}}_{\text{ref}}}-G_{\text{limit}})}\right),\end{array}$$

where S_generation and G_limit are tuning parameters for the logistic function, DD_eff represents the sum of effective temperature for bark beetle reproduction in °C d⁻¹, and DD_ref denotes the thermal sum of degree-days for one bark beetle generation in °C d⁻¹. Saturation of i_{beetles  generation} represents the lack of an available breeding substrate when many generations develop over a short period.

DD_eff is calculated from 1 January until the diapause of the first generation. In ORCHIDEE, a diapause is triggered when the day length exceeds 14.5 h (e.g., 27 April for France). Each day before the diapause with a daily average temperature around the bark above 8.3 °C (T_min) and below 38.4 °C (T_max) is accounted for in the summation of DD_eff (Eq. 10). This approach simulates the phenology of bark beetles, which tend to breed earlier when winter and spring are warmer, thus allowing for multiple generations in the same year (Hlásny et al., 2021a).

$$\begin{array}{}\text{(10)}& \begin{array}{rl}{\text{DD}}_{\text{eff}}& ={\sum }_{n_{\text{diapause}}}^{i=\mathrm{1}}(T_{\text{opt}}-T_{\text{min}})\cdot e^{(\mathrm{0.0288}\cdot T_{\text{bark},i})}\\ & -e^{(\mathrm{0.0288}\cdot b_{\text{eff}}-(\mathrm{40.99}-T_{\text{bark},i})/\mathrm{3.59})}-\mathrm{1.25},\end{array}\end{array}$$

where i is a day and n_diapause is the number of days between the 1 January and the day of the diapause. T_opt (30.3 °C) is the optimal bark temperature for beetle development, and T_min (8.3 °C) is the temperature below which the beetles development stops. T_bark, i is the average daily bark temperature. T_bark, i is calculated as the daily average air temperature with 2 °C subtracted. All parameter values are taken from Temperli et al. (2013).

The bark beetle activity of the previous year (i_{beetles  activity}) is calculated as follows:

$$\begin{array}{}\text{(11)}& i_{\text{beetles\hspace{0.17em} activity}}=\mathrm{1}/\left(\mathrm{1}+e^{-S_{\text{activity}}(\frac{B_{\text{kill},y-\mathrm{1}}}{B_{\text{total}}}-{\text{act}}_{\text{limit}})}\right),\end{array}$$

where i_{beetles  activity} denotes the biomass of the stand damaged by bark beetles in the previous year, B_total is the total biomass of the stand, and S_activity and act_limit are parameters that drive the intensity of this negative feedback.

During the buildup stage, the population of bark beetles can either return to its endemic stage if tree defense mechanisms are preventing bark beetles from successfully attacking healthy trees or evolve into an epidemic stage (Fig. 1) if the tree defense mechanisms fail. During the post-epidemic stage, the forest is still subject to higher mortality than usual, but signs of recovery appear (Hlásny et al., 2021a). Recovery may help the forest ecosystem to return to its original state or switch to a new state (different species, change in the forest structure) depending on the intensity and the frequency of the disturbance (Van Meerbeek et al., 2021).

2.7 Bark beetle survival

The capability of the bark beetles to survive the winter in between two breeding seasons is critical in simulating epidemic outbreaks. During regular winters, winter mortality for bark beetles is around 40 % for the adults and 100 % for the juveniles (Jönsson et al., 2012). In our scheme, this mortality rate is implicitly accounted for in the calculation of the bark beetle survival index (i_{beetles  survival}). A lack of data linking bark beetle survival to anomalous winter temperatures justifies the implicit approach and prevented including this information as a modulator of i_{beetles  survival}. The latter explains why winter temperatures do not appear in Eq. (11). Instead, the model simulates the survival as a function of the abundance of suitable tree hosts, which decreases the competition for shelter and food:

$$\begin{array}{}\text{(12)}& i_{\text{beetles\hspace{0.17em} survival}}=max(i_{\text{hosts dead}},i_{\text{hosts\hspace{0.17em} alive}}).\end{array}$$

The availability of wood necromass from trees that died recently, particularly following windstorms, plays a critical role in bark beetle survival and proliferation. In the year following a windstorm, uprooted and broken trees may offer an ideal breeding substrate for bark beetles, facilitating their population growth.

In Temperli et al. (2013), an empirical correlation between windthrow events and bark beetle susceptibility was established. ORCHIDEE enhances realism by considering the actual suitable hosts (living or recently dead trees) as the primary driver of bark beetle survival. To avoid overestimating bark beetle population growth, max${}_{N_{\text{wood}}}$ has been introduced. Any addition of dead trees beyond max${}_{N_{\text{wood}}}$ is considered ineffective in affecting the bark beetle population. This ensures that an excess of breeding substrate does not artificially inflate beetle numbers.

This relationship is quantitatively represented in ORCHIDEE through the dead host index, i_hosts dead, which is driven by the availability of recent dead trees. The formulation of i_hosts dead is as follows:

$$\begin{array}{}\text{(13)}& i_{\text{hosts dead}}=min\left({\displaystyle \frac{N_{\text{wood}}}{B_{\text{wood}}}}/\underset{N_{\text{wood}}}{max},\mathrm{1}\right).\end{array}$$

Here, N_wood represents the quantity of woody necromass from the current year, B_wood is the total living woody biomass in the stand, and max${}_{N_{\text{wood}}}$ is the threshold of the ratio $N_{\text{wood}}/B_{\text{wood}}$ signifying the maximum level. This index captures the immediate increase in dead trees suitable for bark beetle breeding following a windthrow event. However, it takes about a year for dead wood to lose its freshness and suitability for bark beetle breeding. This is accounted for by excluding woody necromass that is older than 1 year from the i_hosts dead calculation.

max${}_{N_{\text{wood}}}$ can also be considered a parameter that depends on the spatial scale of the simulation. The mortality rate of trees (DR_windthrow) that will trigger an outbreak is very different across spatial scales. Where a relatively high share of dead wood is needed to trigger an outbreak on the patch scale, a much lower share of dead wood suffices on the landscape scale to trigger a widespread bark beetle outbreak. Therefore, these parameters must be set according to the spatial resolution of the simulation experiment.

The index i_{hosts  alive} denotes the survival of bark beetles, which is facilitated by the abundance of suitable trees that reduces the competition among bark beetles for breeding substrates and therefore increases their survival.

$$\begin{array}{}\text{(14)}& i_{\text{hosts\hspace{0.17em} alive}}=i_{\text{beetles mass attack}}\times i_{\text{hosts susceptibility}}\end{array}$$

The number of suitable tree hosts, i_{hosts  alive}, is driven by two factors: (1) the abundance of weak trees, which can be more easily infected by bark beetles. ORCHIDEE does not explicitly represent weak trees, but tree health is thought to decrease with an increasing density given the stand diameter. The index for host suitability is thus calculated by making use of the relative density index (i_{rd  spruce}).

$$\begin{array}{}\text{(6a')}& i_{\text{hosts susceptibility}}=\mathrm{1}/\left(\mathrm{1}+e^{S_{\text{susceptibility}}\cdot (i_{\text{rd\hspace{0.17em} spruce}}-i_{\text{rd\hspace{0.17em} susceptibility}})}\right)\end{array}$$

Equation (6a^′) is close to Eq. (6a), but the parameter S_{susceptibility} has been reduced by a factor of 2 in order to reflect that i_{hosts susceptibility} is more sensitive to i_{rd  spruce} than i_{hosts  competition}. The second factor is (2) i_{hosts  mass  attack}, which represents the ability of bark beetles to attack healthy trees when the number of bark beetles is large enough. This index only depends on the size of the bark beetle population (i_{beetles  pressure}; see Eq. 8).

$$\begin{array}{}\text{(15)}& i_{\text{hosts\hspace{0.17em} mass\hspace{0.17em} attack}}=\mathrm{1}/\left(\mathrm{1}+e^{S_{\text{mass\hspace{0.17em} attack}}\cdot (i_{\text{beetles\hspace{0.17em} pressure}}-B{P}_{\text{limit}})}\right),\end{array}$$

where S_{hosts  mass  attack} and BP_limit are parameters. S_{mass  attack} controls the steepness of the relationship, while BP_limit is the bark beetle pressure index at which the population is moving from endemic to epidemic stage where mass attacks are possible.

The epidemic stage corresponds to the capability of bark beetles to mass attack healthy trees and overrule tree defenses (Biedermann et al., 2019). At this point in the outbreak, all trees are potential targets irrespective of their health. Three causes have been suggested to explain the end of the epidemic phase: (1) the most likely cause is a high interspecific competition among beetles for the tree host when the density is decreasing (decreasing i_{hosts  alive}) (Komonen et al., 2011; Pineau et al., 2017), (2) a series of very cold years will decrease their ability to reproduce (decreasing i_{beetles  generation}), and (3) a rarely demonstrated increasing population of beetle predators (Berryman, 2002). In ORCHIDEE r8627, the first two causes are represented, but the final one, i.e., the predators, are not.

2.8 Tree mortality from bark beetle infestation

When bark beetles attack a tree, the success of their attack will likely depend on the capability of the tree to defend itself from the attack. Trees defend themselves against beetle attacks by producing secondary metabolites (Huang et al., 2020). The high carbon and nitrogen costs of these compounds limit their production to periods with environmental conditions favorable for growth (Lieutier, 2002). The probability of a successful bark beetle attack is driven by the size of the bark beetle population (i_{beetles  pressure}) and the health of each tree. ORCHIDEE, however, is not simulating individual trees but rather diameter classes within an age class. An index of tree health for each age class (i_{hosts  health, a}) was calculated as follows:

$$\begin{array}{}\text{(16)}& P_{\text{success,}a}=i_{\text{hosts\hspace{0.17em} health,}a}\times i_{\text{beetles\hspace{0.17em} pressure}}.\end{array}$$

A tree rarely dies solely from bark beetle damage (except during mass attacks) as female beetles often carry blue stain fungi, which colonize the phloem and sapwood, blocking the water-conducting vessels of the tree (Ballard et al., 1982). This results in tree death from carbon starvation or desiccation. As ORCHIDEE r8627 does not simulate the effects of changes in sapwood conductivity on photosynthesis and the resultant probability of tree mortality, the index of weakened trees (i_{hosts  health, a}) makes use of two proxies similarly to Eqs. (5) and (6) but simplified to be calculated only for one age class at a time:

$$\begin{array}{}\text{(17)}& {\displaystyle }{i}_{\text{hosts\hspace{0.17em} health,}a}={\displaystyle \frac{(i_{\text{hosts\hspace{0.17em} competition,}a}+i_{\text{hosts defense},a})}{\mathrm{2}}},\text{(5a')}& {\displaystyle }{i}_{\text{hosts defense},a}=\mathrm{1}/(\mathrm{1}+e^{S_{\text{drought}}\cdot (\mathrm{1}-{\text{PWS}}_{y,\text{spruce},a}-{\text{PWS}}_{\text{limit}})}).\end{array}$$

Contrary to Eq. (5a), PWS${}_{y,\text{spruce}{,}_a}$ is the plant water stress from the current year.

$$\begin{array}{}\text{(6a'')}& {\displaystyle }{i}_{\text{hosts\hspace{0.17em} competition,}a}=\mathrm{1}/(\mathrm{1}+e^{S_{\text{competition}}\cdot (i_{\text{rd\hspace{0.17em} spruce},a}-i_{\text{rd\hspace{0.17em} limit}})}),\text{(6b'')}& {\displaystyle }{i}_{\text{rd\hspace{0.17em} spruce},a}={\displaystyle \frac{D_{\text{spruce,}a}}{D_{\text{max}}}}.\end{array}$$

To assess the bark beetle damage rate (DR_beetles), B_{beetles kill} has to be divided by B_total.

2.9 Flow of the calculations

The equations presented above contain feedback loops which are visualized in Fig. 2. In ORCHIDEE, these feedback loops are accounted for in subsequent time steps rather than the same time step.

Figure 2 Order of the calculations and feedback in the Ips typographus outbreak model of ORCHIDEE. The numbers correspond to the equation numbers in this study. The dotted line boxes represent the five main concepts of the outbreak model described in Sect. 2.4, 2.5, 2.6, 2.7, and 2.8. The variable name definitions are listed in Table 1.

3 Methods and material

3.1 Model configuration

Given the large-scale nature of the ORCHIDEE, a sensitivity experiment of the bark beetle outbreak functionality was carried out rather than focusing the model evaluation on matching observed damage volumes at specific case studies. Focusing on model sensitivity for a range of environmental conditions is thought to reduce the risk of overfitting the model to specific site conditions (Abramowitz et al., 2008).

ORCHIDEE r8627, including the bark beetle model, was run at the location of eight FLUXNET sites selected to simulate a credible temperature and precipitation gradient for spruce (see below). For each location, the half-hourly meteorological data from the flux tower were gap-filled and reformatted so that they could be used as climate forcing by ORCHIDEE. Boundary conditions for ORCHIDEE, such as soil texture, pH, and soil color were retrieved from the US Department of Agriculture map for the corresponding grid cell. The observed land cover and land use for the grid cell were ignored and set to pure spruce because this study did not investigate the effect of species mixture in the simulation experiments. The resolution of the grid cell chosen for this analysis is 2500 km². Although this corresponds to a high resolution for large-scale simulations with ORCHIDEE, it is a coarse resolution for studying bark beetle outbreaks.

The climate forcings were looped over as much as needed to bring the carbon, nitrogen, and water pools to equilibrium during a 340-year-long spin-up. Following the spin-up, a 100-year simulation was run starting with a windthrow event on the first day of the first year. The results presented in this study come from the 100-year-long simulations. Given the focus on even-aged monospecific spruce forests in regions where spruce growth is not constrained by precipitation, variables such as i_hosts share and i_{hosts defense} were omitted from this study. Note that ORCHIDEE does not account for possible acclimation e.g., temporal changes in bark beetle behavior or bark beetle resistance to external stressors such as winter temperature.

3.2 Selection of locations

Bark beetle populations are known to be sensitive to temperature as they are more likely to survive a mild winter (Lombardero et al., 2000) and tend to breed earlier when winter and spring are warmer than usual, allowing for multiple generations in the same year (Hlásny et al., 2021a; see also Eq. 10 in Sect. 2.6). In order to assess the temperature effect of the bark beetle outbreak model in ORCHIDEE, eight locations in Europe were selected (Table 2) which represent the range of climatic conditions within the distribution area of Norway spruce (Picea abies H. Karst L.), the main host plant for Ips typographus, the bark beetle species under investigation.

Table 2 Climate characteristics of the eight locations used in the simulation experiments. The acronyms refer to the site names used in the FLUXNET database (Pastorello et al., 2020).

For these eight locations, half-hourly weather data from the FLUXNET database (Pastorello et al., 2020) were used to drive ORCHIDEE. Some of these locations (FON, SOR, HES, COL, WET) are in reality not covered by spruce, but all sites are, however, located within the distribution of Norway spruce. In this study, locations were selected to use the observed weather data to simulate a credible temperature and rainfall gradient for spruce. HES location is no longer part of the FLUXNET network, but the previous data are still available and relevant for this analysis.

3.3 Sensitivity to model parameters

The sensitivity assessment evaluates the responsiveness of four key variables (i_{hosts susceptibility}, i_{beetles mass attack}, i_{beetles  generation}, and i_{beetles  activity}) of the Ips typographus outbreak model implemented in ORCHIDEE. The assessment aims to demonstrate the ability of ORCHIDEE to simulate diverse dynamics of bark beetle infestations. The selection of i_{hosts susceptibility}, i_{beetles  activity}, i_{beetles mass attack}, and i_{beetles  generation} was based on two criteria: (1) their substantial influence on the dynamics of the Ips typographus outbreak noted during model development and (2) their independence from direct measurable data, rendering them less suitable for evaluation through literature review.

For each of the four variables, three distinct values were assigned to two parameters labeled “shape” and “limit”. The shape parameter determines the shape of the logistic relationship, with three values tested: (a) shape = −1.0, yielding a linear relationship; (b) $-\mathrm{5.0}<\text{Shape}<-\mathrm{30.0}$, resulting in a logistic curve; and (c) shape = −500.0, turning the logistic relationship into a step function. For the logistic curve, the exact shape value between −30.0 and −5.0 is chosen according to each index under study: (1) $S_{\text{susceptibility}}=-\mathrm{5.0}$, (2) $S_{\text{activity}}=-\mathrm{20.0}$, (3) $S_{\text{mass\hspace{0.17em} attack}}=-\mathrm{30.0}$, and (4) S_generation=5.0. For S_{mass  attack} and S_activity, higher values were chosen because the slope of the logistic curve has significant impact in order to trigger an outbreak.

The second parameter called limit determines the threshold, derived from expert insights, at which the logistic relationship will reach its midpoint value of 0.5 (i_{rd  susceptibility}, BP_limit, Act_limit, or G_limit). For instance, i_{rd  susceptibility} is set to 0.55, indicating i_{hosts susceptibility} midpoint sensitivity (Eq. 6^′). Setting BP_limit to 0.12 results in an i_{beetles mass attack} midpoint when i_{beetles  pressure} is 0.12, selected for its proximity to scenarios where i_hosts dead equals 1.0 (Eq. 14). The variable act_limit was positioned at 0.06, signifying the i_{beetles  activity} midpoint at DR_beetles=6 % from the preceding year, exceeding endemic levels yet not reaching epidemic outbreaks (Eq. 10). Lastly, G_limit is fixed at 1.0, denoting the midpoint for i_{beetles  generation} upon completing one generation annually, underpinning the rarity of bark beetle outbreaks with fewer than one generation per year (Eq. 9). Starting from these reference values, a restrictive simulation was run in which the limit parameter values were reduced by 50 %. Likewise, a permissive simulation was run to test 50 % higher values for limit.

The sensitivity analysis of the model parameters explores 36 (product of three shapes, three limits, and four equations) combinations of parameter values named a set, but the full design of the experiment is 8³=512 sets (eight parameters, with three values for each). This deliberate choice was made because of the computation time cost of a single run. In order to reduce the number of runs from 512 to 36, we had to make simplifications: (1) one equation at the time is studied, reducing the number of sets necessary to realize the sensitivity analysis to nine, and (2) every other parameters from the remaining equation is set to default value; e.g., limit values are set to their reference values, and shape values are set to their a priori assumption (Table 4). The major drawback of this approach is that interaction effects between equations cannot be investigated in the study. Nonetheless, this sensitivity analysis aims to document model behavior, rather than seek precise parameter values which can be achieved with the main effect of each equation only (see Sect. 3.4).

The simulations were run for the THA site, where they were repeated for two prescribed windthrow events with a different intensity, i.e., a DR_windthrow of 0.1 % and 10 %. The effect of the parameters with a negligible windthrow event, i.e., killing only 0.1 % of the trees, was tested to confirm that the selected parameters did make ORCHIDEE simulate a bark beetle outbreak in the absence of windthrow (score 5 in Sect. 3.4).

3.4 Parameter tuning and credibility score

The results of the sensitivity experiment were used to select key model parameters. Selecting the values for the Shape and Limit parameters (see Sect. 3.3) used in the calculation of the variables i_{hosts susceptibility}, i_{beetles mass attack}, i_{beetles  generation}, and i_{beetles  activity} has been carried out in order to reproduce the observed dynamics of bark beetle outbreaks. Observed dynamics were compiled through a literature search for peer-reviewed papers that reported quantitative characteristics of bark beetle outbreaks (Table 3). Five characteristics could be documented and used to calculate the score:

• the delay between the windthrow event and the start of the bark beetle outbreak (score 1);

• the length of the bark beetle outbreak, defined by the number of years required for a bark beetle population to go back to its endemic level (score 2);

• the cumulative number of trees per unit area, killed by the bark beetles at the end of an outbreak (score 3);

• the average tree mortality rate (DR_beetles) during an endemic stage (score 4);

• no outbreak triggered when DR_windthrow=0.1 %, which is a mandatory characteristic in order to avoid outbreak for no specific event (score 5).

Based on Table S1 in the Supplement and the reference range in Table 3, scores are calculated for each parameter set, and values are either 0 or 1. The credibility score (CS) is the sum of four scores multiplied by a fifth score. The CS is computed as follows: CS = (score 1 + score 2 + score 3 + score 4) × score 5. Only parameter sets that achieve a CS of 4 are selected. If multiple parameter values are possible for a given equation, the most frequently selected value is preferred.

Table 3 Literature-based summary of characteristics of large-scale bark beetle outbreaks. Due to data scarcity, the characteristics combine outbreak dynamics of different bark beetle species, different host species, and different locations. The reference range is used to calculate the credibility score (CS) of each set of parameters (see Table S1).

3.5 Sensitivity to climate and windthrow

In this simulation experiment, the influx of fresh dead tree hosts (N_wood) used for bark beetle breeding was controlled by modifying the maximum damage rate of a windthrow event (DR_windthrow) in ORCHIDEE. Seven DR_windthrow values were simulated (i.e., 0.1 %, 5 %, 7.5 %, 10 %, 15 %, 20 %, and 35 %). Given the monotonic nature of the relationships between DR_windthrow and i_hosts dead (Eq. 12), each event triggers a proportional increase in the dead host availability (i_hosts dead) scaling between 0 and 1 (Fig. 3). Through its equations, ORCHIDEE assumes that for damage rates above 20 %, the variable i_hosts dead (N_wood) will always be equal to 1.0. The index i_{rd  spruce}, however, may further decrease with increasing windthrow damage, which makes the 35 % damage rate still interesting to investigate. Although the simulations were run for all DR_windthrow, only four windthrow damage rates, including a windstorm resulting in a 35 % damage rate (Fig. 3), were presented to enhance the readability of the Results section.

Figure 3 Relationship between windthrow damage rate (DR_windthrow) and the dead host index (i_hosts dead). For each site, a DR_windthrow=0.1 % was used as the reference simulation because an endemic bark beetle population is expected following such a low-intensity windthrow event. The four DR_windthrow values shown in blue were selected for the subsequent presentation of the results because they cover the entire range for i_hosts dead.

The main driver of the number of generations a bark beetle population can achieve in 1 year is the number of days with a temperature higher than 8.3 °C during winter time (Temperli et al., 2013), which is the reason why temperature is so important for bark beetle reproduction. By taking REN, THA, WET, and HES, the number of bark beetle generations ranged from 0.8 to 3.5 (Fig. 4), which is similar to the number of generations observed across Europe (Faccoli and Stergulc, 2006; Jönsson et al., 2009, 2011). Limiting the analysis to only four sites simplifies the presentation without affecting the range under investigation.

Figure 4 Average number of bark beetle generations during the 5 years following the wind storm for eight locations along a climate gradient. The HYY location in Finland was selected as the reference for the REN, THA, WET, and HES locations. Only results from the reference and four selected locations (shown in blue) are shown to enhance readability.

For the climate gradient, the simulation for HYY served as a reference since the number of generations is lower than 1, for which no outbreak should happen under any circumstances. Under present climate conditions, an outbreak in HYY should be considered an undesirable model result. Likewise, a DR_windthrow=0.1 % is considered too low to trigger an outbreak and was therefore used as the reference for the wind damage rate tests.

The experiment consisted of 40 simulations, i.e., right sites (including the reference) multiplied by five wind damage rates (including the reference). Although the simulations were also run for SOR, COL, and FON, their results were found to be too similar to the results of selected sites to present them as well. Hence, the Results section presents only 25 out of the 40 simulations. Three output variables were assessed: bark beetle damage rate (DR_beetles), total biomass (B_total), and net primary production (NPP). Total biomass was investigated over 100 years, whereas DR_beetles and NPP were assessed for the first 20 years following a windthrow.

3.6 Continuous vs. abrupt mortality

Where most land surface models use a fixed turnover time to simulate continuous mortality (Pugh et al., 2017; Thurner et al., 2017), ecological reality is better described by abrupt mortality events. An idealized simulation experiment was used to qualify the impact of abrupt mortality on net biome productivity by changing from a framework in which mortality is approximated by a constant background mortality to a framework in which mortality occurs in abrupt, discrete events. The impact of a change in the mortality framework was assessed with an idealized simulation experiment that compared three configurations of ORCHIDEE: (1) a configuration that simulates mortality as a continuous process, labeled “the continuous configuration”, which corresponds to previous versions of ORCHIDEE; (2) a configuration capable of simulating abrupt mortality from windthrow and subsequent bark beetle outbreaks, labeled “the abrupt configuration”; and (3) a configuration in which windthrow is activated but a bark beetle outbreak is implicitly accounted for in the background mortality. This third configuration enabled attributing the impact to windthrow. The effect of simulating abrupt mortality was evaluated over 20-, 50-, and 100-year time horizons.

The impact of changing the mortality framework from continuous to abrupt was quantified on the basis of 120 simulations (eight locations multiplied by seven windthrow damage rates and two configurations, with added eight sites multiplied by one configuration) of 100 years each.

The simulations with abrupt mortality were run first. Subsequently, the number of trees killed was quantified and used as a reference value for the continuous mortality setup. This approach resulted in the same quantities of dead trees at the end of the simulation for both frameworks, which then differed only in the timing of the simulated mortality.  This precaution is necessary to avoid comparing two different mortality regimes where the result would mainly be explained by the intensity of the mortality rather than by its underlying mechanisms.

Changes in forest functioning were evaluated through the temporal evolution of accumulated net biome productivity (NBP) over a 100-year time frame. NBP is defined as the regional net carbon accumulation after considering losses of carbon from fire, harvest, and other episodic disturbances. In ORCHIDEE, NBP is calculated following the definition by Chapin et al. (2006) as the carbon remaining in the biomass, litter, and soil after accounting for photosynthesis and respiration because fire, harvest, leaching, and volatile emissions were not accounted for in this simulation experiment.

4 Results

4.1 Sensitivity to model parameter sets

The impact of spruce stand competition (i_{hosts susceptibility}) on outbreak dynamics was examined by adjusting the parameters S_{susceptibility} and i_{rd  susceptibility} in Eq. (6a^′). When S_{susceptibility} resulted in a linear relationship ($S_{\text{susceptibility}}=-\mathrm{1.0}$), no peak in bark beetle damage occurred for the three tested values of i_{rd  susceptibility} (permissive, reference, and restrictive) at a 10 % windthrow damage rate (Fig. 5h). However, employing a step function $(S_{\text{susceptibility}}=-\mathrm{500.0}$) led to either sporadic peaks of bark beetle damage with a permissive i_{rd  susceptibility} or a 2-year outbreak with a maximum damage rate of 60 % and a restrictive i_{rd  susceptibility} (Fig. 5h), neither of which aligns with the observations summarized in Table 3.

The outcome closest to observations from Table 3 was obtained with a logistic relationship ($S_{\text{susceptibility}}=-\mathrm{5.0}$), where i_{rd  susceptibility} determined the duration of the outbreak: 11, 16, and 25 years for restrictive, reference, and permissive parameter values, respectively (Fig. 5h). Either the restrictive or the reference parameter value could be utilized since a range of 11–16 years aligns with the observations (Table 3). To examine the occurrence of improbable outbreaks, sensitivity tests were repeated for a 0.1 % windthrow damage rate. None of the nine parameter combinations triggered an outbreak (Fig. 5g), suggesting that improbable outbreaks due to the calculation of i_{hosts susceptibility} are unlikely.

From the calculation of the credibility score, only one set obtains a score of 4 ($S_{\text{susceptibility}}=-\mathrm{5.0}$, i_{rd  susceptibility}=0.55; Table S1). The concerning parameter value was selected and is reported in Table 4.

The effect of the capability of bark beetle to mass attack (i_{beetles mass attack}) when the population exceeds a threshold was evaluated by varying S_{mass  attack} and BP_limit (Eq. 14). Linear relationships ($S_{\text{mass\hspace{0.17em} attack}}=-\mathrm{1.0}$) resulted in similar outbreak dynamic for all BP_limit values, with the model settling on a constant endemic damage following an outbreak, though higher than observed (Table 3; Fig. 5f). Introducing a logistic or step function slightly altered outbreak dynamics except when assuming a step function for the restrictive value, which prevented an outbreak. Repeating sensitivity tests for a 0.1 % windthrow damage rate showed that assuming linear or logistic relationships could trigger an outbreak (Fig. 5e), indicating that improbable outbreaks may arise from the calculation of i_{hosts  mass  attack}.

From the calculation of the credibility score, three sets obtain a score of 4, but only set 4.6 was chosen because of its intermediate position compared to sets 4.9 and 4.5 (Table S1). The concerning parameter values ($S_{\text{mass\hspace{0.17em} attack}}=-\mathrm{30.0}$, BP_limit=0.06) were selected and are reported in Table 4.

The impact of bark beetle activities from the previous year (i_{beetles  activity}) on outbreak dynamics was investigated by varying S_activity and act_limit (Eq. 10). Linear or logistic relationships resulted in excessively long outbreaks (>30 years) compared to observations (Table 3; Fig. 5b), whereas assuming a step function relationship simulated a decline in the outbreak after 14 years. Sensitivity tests repeated for a 0.1 % windthrow damage rate showed that assuming a linear relationship could trigger an improbable outbreak (Fig. 5a) through the calculation of i_{beetles  activity}.

From the calculation of the credibility score, only one set obtains a score of 4 ($S_{\text{activity}}=-\mathrm{500.0}$, act_limit=0.12; Table S1). The concerning parameter value was selected and are reported in Table 4.

To explore the effect of the numbers of generation (i_{beetles  generation}) on the outbreak dynamics, S_generation and G_limit from Eq. (9) were varied. The bark beetle damage rate was more sensitive to G_limit than S_generation, but only a linear relationship with the reference G_limit=1.0 yielded an intermediate outbreak intensity consistent with the continental climate at the test location (i.e., THA; Fig. 5d). Other combinations resulted in either too strong of a peak or no peak during the outbreak. Repeating sensitivity tests for a 0.1 % windthrow damage rate showed that none of the nine parameter combinations triggered an outbreak (Fig. 5c), indicating that improbable outbreaks from the calculation of i_{beetles  generation} are unlikely.

From the calculation of the credibility score, three sets obtain a score of 4, but only set 1.4 was chosen because of its intermediate position compared to sets 1.1 and 1.5 (Table S1). The concerning parameter values (S_generation=1.0, G_limit=1.0) were selected and are reported in Table 4.

Figure 5 Simulation results from the sensitivity experiment at the THA site. Eight parameters from four equations were evaluated. Each equation represents an index from the bark beetle outbreak model (i_{hosts susceptibility}, i_{hosts  mass  attack}, i_{beetles  activity}, and i_{beetles  generation}). Each index is represented by a logistic function defined by a shape parameter (Shape) and a limit parameter (Limit). Three values were chosen for each parameter, resulting in nine pairs of parameters for each index. Colored lines represent the shape parameter varying from linear (Shape = −1.0, red; logistic $-\mathrm{5.0}<\text{Shape}<-\mathrm{30.0}$, green) to step function, where Shape = −500.0 (blue). The line type represents three different values for the Limit parameter, where references (dashed line) are values of i_{rd  susceptibility}, BP_limit, act_limit, and G_limit (given in Table 4), whereas permissive (full line) and restrictive (dashed–dotted) represent a 50 % decrease or increase, respectively.

4.2 Model tuning

By comparing the outcomes of the sensitivity tests (Sect. 4.1) to a compilation of observations (Table 3), a first estimate for several parameters was proposed (Table 4). Details about the other values used in the sensitivity test can be found in Table S1.

Table 4 Parameter values from the bark beetle model based on the score obtained in the sensitivity analysis.

^∗ Parameter values deliberately fixed and excluded from the sensitivity analysis (Sect. 3.3 for justification).

4.3 Impact of climate and windthrow on bark beetle damage

In ORCHIDEE, the warmest sites, HES and WET, experienced significant bark beetle outbreaks across a wide spectrum of windthrow mortality rates, whereas colder sites like REN and THA saw outbreaks only in response to the most severe windthrow events (Fig. 6b, c). A greater average number of bark beetle generations in the years following windthrow events led to higher bark beetle damage rates at the peak of outbreaks. For instance, at a 35 % windthrow mortality rate, HES reached a maximum bark beetle damage rate of 50 %, whereas REN's maximum was 22 % (Fig. 6a, b).

Interestingly, high tree mortality rates from windthrow could also lead to delays and a lower maximum DR_beetles (Fig. 6). For instance, at the HES site, 10 %, 20 %, and 35 % windthrow damage rates triggered maximum DR_beetles values of 50 %, 43 %, and 37 %, respectively (Fig. 6a). Conversely, low DR_windthrow, such as that of 5 % at WET, delayed the peak of bark beetle outbreaks by 9 years (Fig. 6d). Additionally, the model simulated a post-epidemic stage during which the outbreak damage rate remained relatively low (<10 %) and lasted between 3 and 10 years (Fig. 6). Overall, the simulated outbreaks lasted between 11 and 20 years, which is consistent with field observations (Table 3).

Figure 6 Simulation results of 16 simulations (four locations with four windthrow damage rate values, DR_windthrow). Lines represent the annual bark beetle damage rate as a fraction of the total biomass (DR_beetles). N_bgen is the average number of bark beetle generations during 5 years after the windthrow event. DR_windthrow represents the percentage of biomass loss by a windthrow event at the start of the simulation.

At the coldest site, HYY, ORCHIDEE simulated only a small number of bark beetle generations, preventing outbreaks from occurring. This observation validates the initial parameter tuning (Table 4), indicating that it is robust enough to prevent improbable outbreaks, such as the model triggering outbreaks in sites where bark beetles cannot reproduce.

4.4 Impact of climate and windthrow on stand biomass and net primary production

All locations experienced a 10- to 20-year decrease in total biomass until at most 9 kg C m⁻², at which time the outbreak ended (Fig. 7a, b, c, d). The model can simulate significant epidemic events even if the initial trigger, such as the windthrow event in our study, is not particularly intense. Once the bark beetles can mass attack living trees, the bark beetle population (i_{beetles  pressure}) will increase and kill more and more trees until so many trees are killed that the stand density of the remaining living trees drops below the threshold of $i_{\text{rd\hspace{0.17em} spruce}}=i_{\text{rd\hspace{0.17em} limit}}=\mathrm{0.4}$. In ORCHIDEE, an i_{rd  limit}=0.4 for spruce forest corresponds to a biomass of around 9 kg C m⁻², which in ORCHIDEE is too low to maintain an epidemic population of bark beetles at the 2500 km² grid cell. Interestingly, for the climate observed at REN, where the number of generations is approximately one, the bark beetle population can only become an epidemic following an intense windthrow event with a 35 % damage rate (Fig. 7).

Throughout the outbreak period, there was a notable decrease in net primary production (NPP) (Fig. 7). This decrease is primarily attributed to a sharp decline in the leaf area index (not shown). Following the epidemic phase, the leaf area recovers. Following the outbreak, the reduction in stand tree density due to bark beetle damage decreases autotrophic respiration (not shown), and the sparser canopy allows for more light to reach the forest floor, where it fosters recruitment (not shown), resulting in a higher NPP or forest growth (Fig. 7). Consequently, carbon use efficiency tends to be higher in sparsely populated stands compared to densely populated ones.

Figure 7 Simulation results of 16 simulations (four sites with four windthrow mortality rate values). Lines represent the annual average net primary production (NPP) in g C m⁻² yr⁻¹ or total stand biomass (B_total) in kg C m⁻². N_bgen is the average number of bark beetle generations during the 5 years after the windthrow event. DR_windthrow represents the percentage of biomass loss by a windthrow event at the start of the simulation. Grey areas represent the epidemic phase.

Figure 8 Difference in cumulative net biome production at three discrete time horizons (i.e., 20, 50, and 100 years) between a fixed continuous mortality rate (blue, n=8), abrupt tree mortality from a windstorm and the subsequent bark beetle outbreak (red, n=56), and abrupt mortality from a windstorm not followed by a bark beetles outbreak (green, n=56). Note that in the continuous mortality configuration, the mortality rate was adjusted to obtain a similar number of trees killed after 100 years as in the abrupt mortality configuration. The variation in each box plot arises due to different locations and prescribed storm intensities. Each box plot displays the median value (thick horizontal line), the quartile range (box border), and the 95 % confidence interval (vertical line). A Wilcoxon test between the three configurations at each time horizon showed significant differences (p value <0.0001) denoted by the four stars.

4.5 Continuous vs. abrupt mortality

The total accumulated net biome production (NBP) was evaluated using the ORCHIDEE model across three different time frames: 20, 50, and 100 years. At the 20-year mark, the average accumulated NBP notably differed between mortality configurations that are continuous, abrupt, and abrupt without a bark beetle outbreak (“no-beetle”): $-\mathrm{7.12}\pm \mathrm{0.97}$, $-\mathrm{1.37}\pm \mathrm{0.28}$ and 3.39±0.74 kg C m⁻² yr⁻¹ for the abrupt, no-beetle, and continuous mortality configurations, respectively. These differences were statistically significant (Wilcoxon, p value <0.0001), indicating a substantial initial reduction in NBP with the abrupt configurations as ecosystems behaved as carbon sources, whereas under the continuous configuration, they acted as carbon sinks (Fig. 8). The variability in NBP demonstrated the broad temperature gradient in Europe and indicated that despite many locations potentially acting as sources under the abrupt configuration, some may transition to carbon sinks within the first 20 years following a disturbance.

Moving to the 50-year horizon, the difference between the three frameworks decreased, with net biome productions of $-\mathrm{0.81}\pm \mathrm{0.60}$, 4.43±0.15, and 5.61±0.18 kg C m⁻² yr⁻¹ for the abrupt, no-beetle, and continuous mortality configurations, respectively. The difference in sink strength remained statistically significant (Wilcoxon; p value <0.001), with the NBP in the abrupt configuration approaching carbon neutrality, while without the consecutive bark beetles outbreak, the ecosystems already became a carbon sink. The climate conditions had a lasting effect on the responses, with the abrupt configuration showing a greater range in responses compared to the continuous one.

At the 100-year mark, the average cumulative NBP values for the abrupt and continuous configurations approached each other with values of 4.85±0.26, 7.09±0.17, and 7.73±0.40 kg C m⁻² yr⁻¹, respectively (Fig. 8) but were still significantly different (Wilcoxon, p value <0.001). ORCHIDEE simulated a return to a carbon sink (indicated by positive cumulative NBP values) suggesting a long-term recovery and potential return to pre-disturbance productivity levels within a century following the windthrow and beetle outbreak event. The continuous configuration displayed a consistently higher median value, suggesting weaker impact of tree mortality dynamics on the long-term carbon cycle.

5 Discussion

5.1 Simulating the dynamics of bark beetle outbreaks and their interaction with windthrow

Our Ips typographus outbreak model has demonstrated its capability to simulate a broad range of disturbance dynamics. The variation in the outbreak dynamics and the response of the outbreak to its main drivers (Figs. 5 and 6) give confidence in the ability of ORCHIDEE to simulate various outbreak scenarios observed across the temperate and boreal zones under changing climate conditions.

Windthrow events have a significant ecological impact because such disturbances offer fresh breeding substrates, which in turn increase bark beetle populations (Lausch et al., 2011). Our model results align with these findings, indicating that windthrows causing damage of 5 % or more may trigger beetle outbreaks (Fig. 6). Additionally, a strong increase in bark beetle populations has been observed following a windthrow event (Wermelinger, 2004), a pattern reflected in the ORCHIDEE simulations. The model simulates a buildup stage spanning 1 to 9 years, where bark beetle numbers increase prior to peaking, with the duration influenced by the severity of the windthrow and the prevailing climate (Fig. 6).

Temperature is another critical factor affecting bark beetle life cycles. Intra- and interannual variation in temperature impact bark beetles, with warmer conditions fostering multiple generations per year, whereas cooler, damp climates slow breeding and survival rates (Benz et al., 2010). In line with these findings, the temperature dependence of the ORCHIDEE simulations show that cold winters at locations such as SOR and REN reduced bark beetle activity compared to warmer locations like THA and WET (Fig. 6). Lieutier et al. (2002) documented that if the population is large enough, bark beetles can mass attack healthy trees. Our model incorporates this dynamic, which is illustrated by epidemic stages where living trees become viable hosts, which then exacerbates the growth of the beetle population (Fig. 1).

The aftermath of a windthrow and subsequent bark beetle outbreak also affects the forest carbon and nitrogen cycles. This impact is observed in the form of snags, which are standing dead trees that undergo decomposition. Snags can temporarily disrupt the link between soil and ecosystem carbon and nitrogen dynamics (Rhoades, 2019; Custer et al., 2020). While in ORCHIDEE the decay of fallen logs does not account for snags yet, the model suggests a recovery period ranging from 5 to 15 years, contingent upon the intensity of the bark beetle outbreak (Fig. 7). As snags create gaps in the canopy, conditions favorable to natural forest regeneration emerge (Jonášová and Prach, 2004).

5.2 Emerging properties from interacting disturbances

While this study did not precisely quantify the impact of simulating abrupt mortality rather than approaching mortality as a continuous process, it demonstrated that the impact of abrupt mortality varies across location and time; i.e., ecosystem functions, such as carbon storage, are affected by natural disturbances like Ips typographus outbreaks, having significant impacts on short- to mid-term carbon balance estimates (Fig. 8). The simulation experiments also highlighted that the legacy effects of disturbances can endure for decades even for a simplified representation of forest ecosystems such as ORCHIDEE, where the recovery might be too fast due to the absence of snags (Senf et al., 2017).

The ability to simulate resistance (i.e., staying essentially unchanged despite the presence of disturbances; Grimm and Wissel, 1997) as an emerging property is evident from Figs. 6 and 7 for REN, where no bark beetle outbreaks were observed following a medium windthrow event (5 %–20 %). However, in all simulated locations that could not resist a bark beetle outbreak, the forest was resilient (i.e., returning to the reference state or dynamic after a temporary disturbance; Grimm and Wissel, 1997), and ecosystem functions were restored to the level from before the windthrow. The elasticity (the speed of return to the reference state or dynamic after a temporary disturbance; Grimm and Wissel, 1997) of the carbon sink capacity ranged from 7 to 14 years. This elasticity is in line with the little observational evidence of ecosystem shifts due to natural disturbances in forests (Millar and Stephenson, 2015). Finally, after the disturbance and the recovery of vegetation structure, the ecosystems simulated by ORCHIDEE showed persistence (i.e., continue along their initial developmental path; Grimm and Wissel, 1997).

5.3 Are cascading disturbances important for carbon balance estimates?

The enhanced complexity introduced in the ORCHIDEE model by incorporating abrupt mortality events as opposed to a continuous mortality prompts the question of whether this model refinement yields new insights into carbon balance estimates. Our century-long analysis demonstrated that the net biome production, a the metric for carbon sequestration, ultimately converges between the continuous and abrupt mortality frameworks (Fig. 8). This suggests that irrespective of the nature of the mortality events, the forest ecosystem goes through a recovery phase, marked by increased growth that compensates for the growth deficits during the disturbance.

However, our experiment has not taken into account the frequency of disturbances. Given the profound influence of disturbance legacies on carbon dynamics, a recurrence interval shorter than the recovery time of the forest might result in a tipping point. Such a scenario could diminish the carbon sequestration potential of the forest beyond the 100-year time frame and, in extreme cases, may even lead to ecosystem collapse, an outcome not explored in the current simulations nor documented in the recent literature (Millar and Stephenson, 2015).

In the mid-term, spanning 20 to 50 years, the widely used continuous mortality model appears to inflate the carbon sink capabilities of forests when juxtaposed with abrupt mortality scenarios. Since policy frameworks, including the Green Deal for Europe (2023) and the Paris Agreement (UNFCCC, https://unfccc.int/process/the-kyoto-protocol, last access: March 2023), were imposed upon these medium-term predictions, they would benefit from adopting model simulations that integrate abrupt mortality events to avoid an overestimation of carbon sink capacities of forest. Furthermore, the accuracy of carbon balance estimates strongly depends on the initial state of the forest in the model. Forest conditions markedly affect carbon uptake rates. Thus, incorporating an abrupt mortality framework into the ORCHIDEE model could substantially refine and strengthen the predictive power of our carbon balance assessments across short-, medium-, and long-term scales.

However, it is important to note that the global-scale NBP might not be significantly impacted by abrupt disturbances compared to continuous ones. This is because the scale is so vast that abrupt mortality events are constantly occurring somewhere on the planet, which smooths out the effects illustrated in Fig. 8.

5.4 Shortcomings of the bark beetle outbreak model

The bark beetle outbreak model developed in this study builds upon the strengths of the previously established LandClim model, though it also inherited some of its limitations. One notable shortcoming is the model for bark beetle phenology, which is an empirical model that makes use of accumulated degree-days. Since the conception of the phenology model a decade ago, Europe's climate has undergone substantial changes, primarily manifested in warmer winters and springs (Copernicus, 2024). Because of these changes, chances have increased that two or even more bark beetle generations occur within a calendar year (Hlásny et al., 2021a). These changes call for an update of the beetle's phenology model to align with these more recent observations (Ogris et al., 2019).

A second limitation is that our study, ORCHIDEE, has been parameterized to simulate only Ips typographus in Europe. In order to change the beetles and tree host interactions e.g., pine bark beetle in North America (Dendroctonus monticolae Hopkins), the sensitivity of indexes must be revised; for example, the pine beetle does not breed on the dead wood falling from windthrow but is very sensitive to drought events (Preisler et al., 2012). The indexes i_{hosts defense} and i_hosts dead as well as the phenology model will need to be revised.

Another issue is the model's consideration of drought. As outlined in the “Methods and material” section, drought is treated as an exacerbating factor rather than a primary trigger, as is the case for windthrow. This understanding was accurate for Ips typographus a decade ago (Temperli et al., 2013); however, emerging evidence increasingly suggests that drought events may indeed trigger bark beetle outbreaks across Europe (Nardi et al., 2023; Netherer et al., 2015). Consequently, this extreme drought as a trigger should be incorporated in a future revision of ORCHIDEE's Ips typographus outbreak model.

6 Outlook

This study simulated the one-way interaction between windthrow and Ips typographus outbreaks in unmanaged forests. Future research will incorporate additional interactions, such as the interplay between droughts, storms, and bark beetles; storms, bark beetles, and fires; and forest management, storms, and bark beetles.

The bark beetle outbreak model could also be enhanced by simulating: (a) standing dead trees (or snags), which would help account for differences in wood decomposition between snags and logs (Angers et al., 2012; Storaunet et al., 2005); (b) migration of bark beetles to neighboring locations, which becomes significant to account for in a model that operates at spatial resolutions below approximately 10 km; and (c) an up-to-date beetle phenology model which accounts for the recent change in their behavior induced by climate change.

This research provided an initial qualitative assessment of a new model feature. However, the application of the model necessitates an evaluation of the simulations against observations of cascading disturbances on the regional scale, which is the topic of an ongoing study.

7 Conclusion

Our approach enables improving the realism of the Ips typographus model in ORCHIDEE without reducing its generality (Levins, 1966). The integration of a bark beetle outbreak model in interaction with other natural disturbances such as windthrow into the ORCHIDEE land surface model has resulted in a broader range of disturbance dynamics and has demonstrated the importance of simulating various disturbance interaction scenarios under different climatic conditions. Incorporating abrupt mortality events instead of a fixed continuous mortality calculation provided new insights into carbon balance estimates. The study showed that the continuous mortality framework, which is commonly used in the land surface modeling community, tends to overestimate the carbon sink capacity of forests in the 20- to 50-year range in the region/area under the same disturbance pressure compared to scenarios with abrupt mortality events.

Apart from these advances, the study revealed possible shortcomings in the bark beetle outbreak model, including the need to update the beetle's phenology model to reflect recent climate changes and the need to consider extreme drought as a trigger for bark beetle outbreaks in line with emerging evidence. Looking ahead, future work will further develop the capability of ORCHIDEE to simulate interacting disturbances such as the interplay between extreme droughts, storms, and bark beetles and between storms, bark beetles, and fires.

The final step will be a quantitative evaluation based on observed data (Marini et al., 2017) in order to assess the capability of ORCHIDEE to simulate complex interaction between multiple sources of tree mortality affecting the carbon balance on a large scale.