An evaluation of multi-species empirical tree mortality algorithms for dynamic vegetation modelling

Tree mortality is key for projecting forest dynamics, but difficult to portray in dynamic vegetation models (DVMs). Empirical mortality algorithms (MAs) are often considered promising, but little is known about DVM robustness when employing MAs of various structures and origins for multiple species. We analysed empirical MAs for a suite of European tree species within a consistent DVM framework under present and future climates in two climatically different study areas in Switzerland and evaluated their performance using empirical data from old-growth forests across Europe. DVM projections under present climate showed substantial variations when using alternative empirical MAs for the same species. Under climate change, DVM projections showed partly contrasting mortality responses for the same species. These opposing patterns were associated with MA structures (i.e. explanatory variables) and occurred independent of species ecological characteristics. When comparing simulated forest structure with data from old-growth forests, we found frequent overestimations of basal area, which can lead to flawed projections of carbon sequestration and other ecosystem services. While using empirical MAs in DVMs may appear promising, our results emphasize the importance of selecting them cautiously. We therefore synthesize our insights into a guideline for the appropriate use of empirical MAs in DVM applications.


Results
MA behaviour under present climate. Across the six species, large differences in simulated mortality patterns occurred under present climatic conditions (Fig. 1). MA structure (i.e., the explanatory variables considered to predict mortality) was an important determinant of the simulated tree cohort half-life time (expressed as MT 50% , i.e., the mean time until 50% of the initially present trees died) and response to competition (expressed as ∆MT 50% , the relative percentage change in MT 50% compared to a 'no competition' scenario). With the exception of Betula pendula, 'CI-based MAs showed higher MT 50% than 'Growth-based' and 'Size-only'-based MAs, and were highly sensitive to changes in competition (Fig. 1). This pattern occurred consistently at both sites (i.e., the moister site Bern, Fig. 1, and the more xeric site Basel, Fig. S2.1) and for all investigated tree sizes (Fig. 1). The difference between 'CI-based' , 'Growth-based' and 'Size-only-based' MAs was most pronounced for low values of competition and decreased with increasing competition intensity (Fig. S2.4). A notable exception to the difference between the MA groups was found for Picea abies, where the 'CI-based' MA by Monserud & Sterba (1999) behaved similarly as the 'Growth-based' MAs ( Fig. 1). At the species level, no clear distinction in terms of tree cohort half-life times and response to competition was evident between early-and late-successional species due to high within-species variability ( Fig. 1 and Fig. S2.1). Simulations  www.nature.com/scientificreports/ the 'warmer' and 'warmer and moister' scenario (e.g., mean ∆MT 50% of −3.4% and + 1%, respectively for Basel) compared to present climate conditions (Fig. 2, Fig. S2.2). At the more xeric site (Basel), a 'warmer and drier' climate however resulted in a considerable increase in mortality for the 'Growth-based MAs' (mean ∆MT 50% of + 11% over all 'Growth-based' MAs), while the opposite response occurred for the 'CI-based' MAs (less mortality, mean ∆MT 50% of −41%) (Fig. 2). Results for the mesic site (Bern) showed similar patterns, but of a much smaller magnitude (mean ∆MT 50% of + 1% and −7%, respectively, see Fig. S2.2). Not surprisingly, MAs based solely on tree size ('Size-only') showed only a minimal and inconsistent response to the CC scenarios ( Fig. 2). At the level of the tree species, no clear relationship emerged between ∆MT 50% and species drought tolerance, even within the group of 'Growth-based' MAs ( Fig. 2).

MA behaviour under climate change.
MA performance: empirical evaluation. Simulated stand structure (basal area and stem density) in dynamic equilibrium showed large differences between MAs and partly implausible stand structures when compared with empirical data from old-growth forests (Fig. 3). Substantial overestimations to the point of implausible basal areas for old-growth forest (e.g., > 100 m 2 per ha) occurred for MAs of the 'CI-based' as well as for the 'Growth-based' group. When comparing the percentage bias between simulated and observed forest structure across all species, 'Growth-based' MAs showed an overprediction of mean basal area by 36% and 'CI-based' MAs an overprediction of 162% across all species (Table S2.9). At the level of individual species, implausible overestimations of basal area were found for all species, except for Quercus petraea (where however only 3 MAs were available).

Figure 1.
Comparison of model behaviour based on the projected tree cohort half-life times (MT 50% , time (years) until 50% of the initially present trees died) and the response to competition (∆MT 50% , i.e. the relative change of MT 50% in percentage in response to a increase in BAL by 10 m 2 ha −1 compared to a 'no competition' scenario) for six tree species at the mesic site (Bern). Note that for response to competition, positive values indicate a relative increase in mortality and negative values a decrease in mortality. Symbol size indicates the three initial tree sizes used in the simulations, and symbol color indicates MA structure (differentiating between 'Growth-based' , 'Competition-index' (CI)-based and 'Size-only' based MAs). Numbers refer the MA sources: 1: ForClim default MA (based on theoretical assumptions 9 ), 2: Hülsmann, et al. 18

Discussion
In contrast to our expectations, the results showed (1) high variability and inconsistent mortality patterns between alternative MAs for multiple species under current climatic conditions; (2) opposite responses to climate change whose directions were associated with MA structure but independent of species characteristics; and (3) high variability of simulated forest structure for all species, which in most cases deviated substantially from empirically supported ranges. We therefore recommend that empirical MAs for multiple species should be applied in DVMs only if their behaviour has been tested for different scenarios in advance, particularly once they are applied in a climate change context. Below, we discuss each aspect in more detail and conclude with specific recommendations for selecting appropriate empirical MAs for DVM applications (Table 1). Under present climatic conditions, several studies 12,37,46 have shown that empirical MAs based on large inventory datasets are usually well-suited to capture the relationship of mortality with tree size, competition and site conditions for a wide range of species, thus making them promising candidates for improving the empirical foundation of DVMs 46 . A previous study by Hülsmann, et al. 37 analysed empirical MAs from different sources outside a DVM framework and found consistent species-specific patterns of mortality despite large differences in the underlying methodology (i.e., model structure, parameterization), datasets and geographic origin. However, when applying MAs from the same sources as Hülsmann, et al. 37 within a DVM framework, we found highly variable and largely inconsistent results at the level of individual species, indicating complex feedbacks between MAs and DVMs. Two major reasons complicate the straightforward application of empirical MAs in DVMs.
First and most evidently, the difference in MA structure ('CI-based' , 'Growth-based' , 'Size-only based') led to considerably different mortality responses in the DVM. CI-based MAs are typically well-suited to capture mortality primarily in response to the effects of changes in competition 19,44 , a major driver of mortality in dense young www.nature.com/scientificreports/ stands undergoing self-thinning, but less so in old-growth forests 47 . In contrast, 'Growth-based' MAs consider competition and other environmental factors via decreased changes in radial growth rates 46 , reflecting the fact that reduced growth is often linked to a higher mortality probability 16,48 . 'Growth-based' MAs thus appear to be more flexible in capturing mortality processes, e.g. in older life stages (for instance drought impacts in larger trees 31 ) and thus may be structurally better suited for DVM applications over longer simulation intervals 22 . A second reason for the large differences of MA performance lies in the datasets used for their calibration. In particular, some of the MAs in our study covered different diameter ranges 37 , which can lead to substantially different mortality predictions (e.g., due to different mortality causes for different life stages 20 ). For instance, the calibration dataset for the MA of Picea abies by Monserud and Sterba 39 included observations from large trees Figure 3. Comparison of forest structure (basal area and stem density, dbh > 7 cm) simulated by ForClim in dynamic equilibrium at the European study sites (see Appendix S1 and Fig. S1.1 for details) with measured data from old-growth forests (black dots indicating empirical measurements and grey areas indicating the envelope of measured data). Sources for empirical data for the different species are provided in Appendix S1. www.nature.com/scientificreports/ that featured a high mortality probability (i.e., a U-shaped function of DBH 12 ), thus resulting in a more realistic basal area in the dynamic equilibrium simulations (Fig. 3). Similarly, the 'Size-only' based MAs included in our study were based on nation-wide datasets including also large trees (e.g., MAs for Abies and Quercus 39 ), thus explaining their relatively good performance (i.e., basal area close to empirical observations) for the dynamicequilibrium simulations.
In summary, we found that large differences in the datasets as well as in the structure and parameterization of the MAs overruled species-specific mortality responses in the DVM application under current climate. Although the MAs show species-specific patterns when applied within the range of their calibration data (e.g., the observed diameter range and stand structure) outside a DVM 37 , our findings emphasize that substantial extrapolation effects 23 can occur when applying empirical MAs in a DVM framework 9,49 .
Under future warmer and drier climatic conditions, our expectation was to find mortality responses in accordance with the species' ecological characteristics, in particular their drought tolerance 4,31 . In contrast, the observed response was driven mainly by MA structure. For 'Growth-based' MAs, the drought-induced reduction in growth under the 'drier' scenario resulted in an increase in mortality, which is in agreement with observations from irrigation experiments 50 , long-term monitoring 51 and studies along environmental gradients 1,4,51 . For the 'CI-based' MAs, however, a growth reduction led predominantly to a decrease in competition intensity (relative to the present climate), which resulted in a decrease in mortality, since other environmental variables (e.g., site index) were insufficient to capture the effect of increasing drought 8 . This pattern has previously been found for Pinus sylvestris 22 , but its generality in the context of climate change responses has not been demonstrated in a DVM context for a wider suite of species.
The absence of a relationship between a tree species' drought tolerance and its mortality response to climate change for the 'Growth-based' MA group is contrary to empirical evidence, which suggests drastically different mortality rates of different tree species under a warmer and drier climate 1 , such as a high drought-sensitivity of Picea abies 52 or a lower sensitivity of the drought-tolerant species Quercus petraea and Pinus sylvestris 53 . The absence of this pattern in our analysis is likely related to the large differences in the datasets underlying the calibration of the MAs 37 , which can contain complex site-specific patterns 32 or effects of biotic mortality agents which can confound the patterns related to single traits 6 . This problem can partly be overcome if MAs were derived using longer-term datasets that span larger regions, thus better capturing climate-change related signals such as drought 24,31 , and by using inverse calibration approaches to better integrate empirical data in a DVM framework 14 . Our findings thus underline the need for further research in developing empirical mortality models for multiple species with a suitable foundation for climate change applications in DVMs. Besides developing empirical MAs, the improvement of mechanistic, eco-physiological MAs that are applicable across species and regions will also be central, particularly for DVM applications under 'no-analogue' climatic conditions (e.g., unprecedented warming or drought conditions) 5,23 .
The evaluation of simulated forest dynamics with empirical data from old-growth forests revealed large variations in simulated forest structure, which in most cases substantially exceeded the observed data ranges. The partly large overestimation of basal area is of particular concern for DVM applications, since it implies similarly overestimated timber stocks in mature forests, an increasingly important topic in the context of the bio-economy 54 , and implausible values of carbon sequestration in living biomass 11 . Overestimation of basal area and biomass is furthermore problematic because DVMs are progressively used as tools for projecting biodiversity and ecosystem service provisioning of forests 55,56 via indicators derived from simulated forest structure 57 . Underestimations of mortality for mature forests, as observed for some MAs in our study, may therefore lead to flawed conclusions for biodiversity (e.g., amount of deadwood 58 ) and ecosystem service provisioning. Our findings thus highlight the importance of research on tree mortality in mature and old-growth stages, where besides competition, direct ecophysiological impacts that can be mediated strongly by biotic agents 6 play a central role 47,59,60 .
Although inventory-based variables such as tree size and stand density can be more robust predictors of mortality rates than physiological indicators 15 , our study implies that the use of empirical MAs for multiple species based on large-scale datasets is no guarantee by itself for improving the robustness of DVM projections. We thus advise caution towards their use and emphasize that empirical MAs should be selected carefully according to the specific DVM application (cf. Table 1). Based on our findings, we consider three aspects to be critical for selecting an appropriate MA: (1) the timespan of the simulation (since the importance of the MA typically increases with DVM simulation time 61 ), (2) the climatic conditions the DVM is applied to, and the degree to which they differ from the conditions the MA was developed for 22 , and (3) the species set used in the DVM simulation (see Table 1). We furthermore emphasize the importance of the spatial scale of the DVM application (e.g., from local to global) and region of application, which need to be considered when selecting empirical MAs (see also Adams, et al. 23 and Hülsmann,et al. 37 ). These topics were however beyond the scope of the present study and require further research. Given that the availability of empirical MAs is still limited, further MAs should be developed, preferably covering additional species while providing structurally similar models. Such developments of empirical MAs should aim at more flexible mortality models based on extensive, long-term datasets along broad environmental gradients 62 that may become increasingly available via novel forest monitoring and data homogenization activities.
We conclude that even small differences in species-specific mortality probabilities can lead to drastically different DVM projections, which can result in biased and erroneous projections of future forest structure 9 , successional trajectories, species composition 34 and ecosystem functioning, e.g., carbon uptake 11 . A prudent selection and further development of empirical MAs is therefore an important step towards improving the robustness of the projections of DVMs, from the stand to the global scale. Considering the recent rapid increase in global datasets about tree demography 24,63 , empirical MAs have a large potential in bridging the knowledge gaps in mechanistic understanding and considerably improve future DVM projections alongside process-based physiological mortality models 23  which was originally developed for the simulation of forest dynamics in Central Europe and has been continuously extended for a wide range of applications (short-and long-term simulations of managed and unmanaged forests) across Europe and beyond 61,64 . ForClim is based on the so-called 'gap model' approach 65 and shares large conceptual similarity with many other DVMs 66 . In this approach, growth, mortality and regeneration of tree cohorts (i.e., groups of trees of the same age and species) are simulated in response to environmental conditions within small patches, based on the theory of patch dynamics 67 . Tree growth is based on the concept of 'constrained optimum growth' 68 , which means that growth under optimum conditions is reduced by environmental constraints (i.e., availability of light, water, temperature and nutrients) based on the species' environmental tolerances (e.g., shade-or drought-tolerance). Tree mortality is by default expressed via a 'theoretical' MA, representing a combination of stress-related and size-dependent mortality 35 , see also Appendix S1.
Empirical mortality algorithms (MAs). We investigated a total of 22 inventory-based empirical MAs for six species (Pinus sylvestris, Picea abies, Abies alba, Betula pendula, Quercus petraea, Fagus sylvatica) from various sources 18,20,[38][39][40][41][42][43][44][45] . All empirical MAs were based on logistic regression models 69 where mortality probability (p) of a tree i is expressed as with X i denoting the design matrix of the linear predictors and β the respective parameter vector. Almost all MAs use tree size (typically diameter at breast height, DBH) as a predictor variable and treat competition either directly using a competition index (typically basal area of larger trees, BAL; 'CI-based' MAs,) or indirectly using tree radial growth as a predictor ('Growth-based' MAs) 19 . Similarly, the effect of environmental conditions is included indirectly in the case of the Growth-based MAs and often directly via a site index in the case of the CIbased MAs. Site index is a typical measure of site quality in forestry and is expressed as the expected height of trees at a certain reference age 19 . An overview of all empirical MAs, their origin, diameter ranges and respective predictors is provided in Appendix S1,  (1) and (2), a full factorial design with the factors competition, tree size and climate was used in order to test the MAs for a wide range of conditions under clearly defined boundary conditions. For this purpose, a cohort of initial trees was planted at the beginning of the simulation (250 trees per ha, see Thrippleton, et al. 22 Table S1.1 and Hülsmann, et al. 37 . To avoid confounding effects from tree regeneration, no regeneration was simulated in steps 1 and 2. Simulations under present climate and climate change were carried out at two sites, representing mesic conditions with near-optimal growing conditions (Bern) and more xeric conditions (Basel; see Appendix S1 for further details). At both sites, four climate scenarios were investigated (with the levels: 'present' , 'warmer' , 'warmer, moister' and 'warmer drier' climate), assuming a change of annual temperature and precipitation of + 4 °C for the 'warmer' scenario (no precipitation change), + 4 °C and + 20% precipitation for the 'warmer, moister' scenario as well as + 4 °C and −20% precipitation for the 'warmer, drier' climate scenario. For the transition to these future climates, a linear change of temperature and precipitation for 100 years (corresponding to the time period of year 2000 to 2100) was assumed until a new constant climate was reached for the rest of the simulation period (see Bircher,et al. 9 ). The assumption of a temperature increase by 4 °C and a precipitation decrease by 20% was based on the 'high impact' RCP 8.5 projections until the end of the twenty-first century for Central Europe 70 . The assumption of a precipitation increase by 20% for the 'warmer, moister' scenario was based on the relatively large uncertainty regarding precipitation changes until the end of the twenty-first century, which also include the possibility of a precipitation increase at this order of magnitude 70 . For details about the present climate see Appendix Table S1.2, for details about the climate change simulations, cf. Appendix S1.
For step (3), the emergent forest structure in dynamic equilibrium was evaluated for each species, assuming monospecific stands. Simulations were started from bare ground (in contrast to the planted initial tree cohorts in steps 1 and 2) for a time span of 1000 years, under present climate, unmanaged conditions and assuming natural regeneration, as done in other long-term DVM applications 71 . Simulated forest structure at the end of the simulation time was compared to observed ranges (basal area and stem density) from largely undisturbed, old-growth forest reserves in Europe (see Appendix S1 for the respective references and Fig. S1.1 for their locations). For each species, the most representative study site was selected among the forest reserves available for each species to conduct the long-term simulations (see Table S1.2). The selection was based on site characteristics reflecting old-growth forest conditions and a clear dominance by the respective species of interest. Details about the study sites and the empirical data can be found in Appendix S1. www.nature.com/scientificreports/ Evaluation of MA behaviour and accuracy. MA behaviour for simulation steps 1 and 2 was characterized by two indicators as in Thrippleton, et al. 22 : (1) the tree cohort half-life time (MT 50% ), defined as the mean time until 50% of the initially present trees died, and (2) the response of the MA to changing competition or climate conditions (∆MT 50% ), measured as the relative percentage change in MT 50% to a baseline scenario; positive ∆MT 50% values indicate higher mortality in comparison to a scenario without competition or climate change, respectively. For the evaluation of the long-term DVM simulations (step 3), simulated forest structure was compared to observed ranges in empirical measurements (basal area and stem density) from largely undisturbed, old-growth forest reserves (see above and Appendix S1 for further details). All analyses and visualizations were carried out using R version 3.6.1 72 .