2.1 Data

2.2 The forest gap model ForClim

ForClim is a dynamic vegetation model that simulates the processes of growth, mortality, and regeneration (often also called “establishment”) of individual trees via species-specific size cohorts (Bugmann, 1994). ForClim is classified as a forest gap model (Shugart, 1984) and simulates forests in independent patches, each with a size of 800 m². By default, 100 patches (i.e., 8 ha) comprise a forest stand, which is used to obtain realistic averages of forest dynamics across patches. The model uses an annual time step and represents trees as cohorts with the properties number of trees (Trs), diameter at breast height (DBH), height, leaf area, and stress level. Here, we used two variants of the regeneration module within ForClim v4.0.1 (Huber et al., 2020).

2.2.1 The two regeneration models

The ForClim regeneration module initiates new cohorts of trees based on (a) site variables for climate and soil in combination with (b) species traits and (c) state variables of forest structure. The species traits of drought tolerance, temperature, and light requirements originate from the indicator values of Ellenberg (1986) (see the latest translated edition, Leuschner and Ellenberg, 2017) and the FORECE model (Kienast, 1987). They define thresholds (so-called “establishment flags”, EFs) that must be fulfilled for a species to qualify for establishment at a DBH of 1.27 cm. For example, if a species has an EF of 10 % of the available light to regenerate, the EF will be fulfilled if the available light is ≥ 10 %. If the available light is < 10 %, the species EF is not fulfilled (see the detailed description of establishment flags below). In this study, we focus on a simple variant and a more complex variant of the regeneration model. Below, a brief summary of the two models is provided, followed by an explanation of the EFs investigated here. For more details, see Fig. 1, Appendix B, and the original documentation for the simple and complex models in Bugmann (1994) and Huber et al. (2020), respectively.

The simple model simulates tree regeneration for each species independently and corresponds to the original ForClim establishment module in Bugmann (1994), which is the same as model variant 1 in Huber et al. (2020). EFs in this model indicate either “suitable” or “not suitable”; i.e., they are binary. In a first step, the annual regeneration probability (kEstP), modulated by species-specific EFs, determines whether regeneration takes place for each species. Second, if the regeneration of a species takes place, the potential maximum number of new trees for that species is calculated from (1) a regeneration intensity parameter (kEstDens, which is the maximum tree establishment density per species [m⁻² yr⁻¹]), and (2) the species-specific successional strategy (i.e., shade-tolerant species have a higher number of seeds and thus offspring compared to shade-intolerant species). Third, the actual number of new trees per species is derived by drawing a random number between 1 and the potential maximum number of trees for each species.

The complex model includes a mechanism for competition and was first introduced as variant 11 in Huber et al. (2020). EFs in this model are continuous, which allows for a more nuanced gradient from “suitable” to “not suitable”. In the first step, kEstP, as modulated by a drought index and degree-days, is used to determine if regeneration takes place for any species. Second, if regeneration does take place, the total potential number of new trees over all species is calculated from (1) a regeneration intensity parameter (kTrMax, which is the absolute maximum number of trees [ha⁻¹]) and (2) a drought index, degree-days, and the continuous EFs. The actual number of new trees over all species is then derived by drawing a random number between 1 and the potential maximum number of trees over all species. Third, the number of new trees per species is calculated by multiplying the actual number of trees over all species by the species-specific ratio of each species' EF and the sum of the EFs of all species.

Figure 1Simplified visual representation of the simple and complex regeneration model variants of ForClim. A full description of the models is provided in Sect. 2.2 below.

Download

2.2.2 Establishment flags regarding light, temperature, and soil moisture

In the present study, we focus on three of the five EFs that are used in the two models (Table 1). The definitions of these three EFs (for light, drought, and degree-days) are given below.

The available light establishment flag (ALEF) evaluates whether the sunlight available at the forest floor (gAL₀; see p. 63 in Bugmann, 1994 for details) matches a parameter for the species' light requirements to regenerate (kL_y,s). kL_y,s is derived from indicator values (ranging from 1 to 9) regarding the light requirements of young trees (L_s) (Leuschner and Ellenberg, 2017). Here,

$$\begin{array}{}\text{(1)}& {\mathrm{kL}}_{y,s}=\left\{\begin{array}{ll}\mathrm{0.025}\cdot \left({\mathrm{L}}_s-\mathrm{1}\right)& L_s<\mathrm{5}\\ \mathrm{0.1}\cdot L_s-\mathrm{0.4}& \text{else}.\end{array}\right.\end{array}$$

For each species s, the binary EF (ALEF_b,s) in the simple model is calculated with

$$\begin{array}{}\text{(2)}& {\mathrm{ALEF}}_{\mathrm{b},s}=\left\{\begin{array}{ll}\mathrm{0}& {\mathrm{gAL}}_{\mathrm{0}}<{\mathrm{kL}}_{y,s}\\ \mathrm{1}& \text{else},\end{array}\right.\end{array}$$

while the continuous EF (ALEF_c,s) in the complex model is calculated with

$$\begin{array}{}\text{(3)}& \begin{array}{rl}& {\mathrm{ALEF}}_{\mathrm{c},s}=\\ & \left\{\begin{array}{ll}\mathrm{0}& {\mathrm{gAL}}_{\mathrm{0}}<{\mathrm{kL}}_{y,s}\\ \mathrm{1}& {\mathrm{gAL}}_{\mathrm{0}}\ge {\mathrm{kL}}_{y,s}+\mathrm{0.5}\\ \frac{{\mathrm{gAL}}_{\mathrm{0}},-{\mathrm{kL}}_{y,s}}{\mathrm{0.5}}& \text{kLy}\le {\mathrm{gAL}}_{\mathrm{0}}<{\mathrm{kL}}_{y,s}+\mathrm{0.5},\end{array}\right.\end{array}\end{array}$$

where the value 0.5 refers to the highest kL_y,s (see Eq. 1) and serves as a buffer for the transition from 0 to 1.

The degree-days establishment flag (DDEF) evaluates whether the annual degree-day sum (gDD; see p. 81 in Bugmann, 1994 for details) matches the species' minimum degree-day requirement (kDDMin). The values of kDDMin originate from Kienast (1987), who derived climatic variables from Müller (1982) and Rudloff (1981) for multiple geographic locations and elevations within the species' distribution range (Ellenberg and Klötzli, 1972; Meusel et al., 1965). This approach was further improved by applying a site-specific bias correction (Bugmann, 1994). Note that this parameter has never been modified to reflect possible deviations regarding the regeneration niche. In this study, we distinguish between the original parameter (kDDMin of adults) and kDDMin_y, which applies to the regeneration. For each species s, the binary EF (DDEF_b,s) in the simple model is calculated with

$$\begin{array}{}\text{(4)}& {\mathrm{DDEF}}_{\mathrm{b},s}=\left\{\begin{array}{ll}\mathrm{0}& {\mathrm{kDDMin}}_{y,s}\ge \mathrm{gDD}\\ \mathrm{1}& \mathrm{else},\end{array}\right.\end{array}$$

while the continuous EF (DDEF_c,s) in the complex model is calculated with

$$\begin{array}{}\text{(5)}& \begin{array}{rl}& {\mathrm{DDEF}}_{\mathrm{c},\mathrm{s}}=\\ & \left\{\begin{array}{ll}\mathrm{0}& {\mathrm{kDDMin}}_{y,s}\ge \mathrm{gDD}\\ \mathrm{1}& \mathrm{gDD}\ge {\mathrm{kDDMin}}_{y,s}+\mathrm{256}\\ \frac{\mathrm{gDD}-{\mathrm{kDDMin}}_{y,s}}{\mathrm{256}}& {\mathrm{kDDMin}}_{y,s}<\mathrm{gDD}<{\mathrm{kDDMin}}_{y,s}+\mathrm{256},\end{array}\right.\end{array}\end{array}$$

where the value 256 refers to the lowest kDDMin_y,s and serves as a buffer for the transition from 0 to 1.

Lastly, the soil moisture establishment flag (SMEF) evaluates whether the drought index (gDr), defined as the ratio of actual evapotranspiration and water demand by the atmosphere (i.e., potential evapotranspiration), matches the species' threshold for this index, i.e., the drought tolerance (kDrTol). The original trait values for drought tolerance range from 1 to 5 (Leuschner and Ellenberg, 2017) and were scaled between 0.06 and 0.3 (i.e., 30 %). The evolution of the formulation of the drought index, including its integration into the regeneration model (Didion et al., 2009a), is documented in Bugmann (1994), Bugmann and Cramer (1998), and Bugmann and Solomon (2000). Similar to the DDEF, the original parameter (kDrTol of adults) and kDrTol_y are distinguished here, and the latter applies to regeneration only. For each species s, the binary EF (SMEF_b,s) in the simple model is calculated with

$$\begin{array}{}\text{(6)}& {\mathrm{SMEF}}_{\mathrm{b},s}=\left\{\begin{array}{ll}\mathrm{0}& \mathrm{gDr}>{\mathrm{kDrTol}}_{y,s}\\ \mathrm{1}& \text{else},\end{array}\right.\end{array}$$

while SMEF_c,s in the complex model is calculated with

$$\begin{array}{}\text{(7)}& \begin{array}{rl}& {\mathrm{SMEF}}_{\mathrm{c},\mathrm{s}}=\\ & \left\{\begin{array}{ll}\mathrm{0}& \mathrm{gDr}>{\mathrm{kDrTol}}_{y,s}\\ \mathrm{1}& \mathrm{gDr}\le {\mathrm{kDrTol}}_{y,s}-\mathrm{0.08}\\ \frac{{\mathrm{kDrTol}}_{y,s}-\mathrm{gDr}}{\mathrm{0.08}}& {\mathrm{kDrTol}}_{y,s}>\mathrm{gDr}>{\mathrm{kDrTol}}_{y,s}-\mathrm{0.08},\end{array}\right.\end{array}\end{array}$$

where the constant of 0.08 indicates the lowest kDrTol_y,s and serves as a buffer for the transition from 0 to 1, i.e., the EF being fulfilled or not (see Huber et al., 2020).

The two EFs in the model that are not considered here are the winter temperature establishment flag (WTEF), which depends on the minimum tolerated winter temperature and chilling requirements (Bugmann, 1994; Bugmann and Cramer, 1998; Kienast, 1987), and the browsing pressure flag, which depends on the species' susceptibility to ungulate browsing (Didion et al., 2009b). WTEF is correlated with DDEF and excluded from the calibration to avoid too many degrees of freedom. We therefore used the default parameterization for WTEF. Because no site information on browsing pressure was available, we decided against using this factor in the calibration. Instead, we kept browsing pressure constant across all sites at its default value of 20 %.

Table 1Description of ForClim model parameters that are considered for calibration.

Download Print Version | Download XLSX

2.3 Calibration approach

We used Bayesian inference to estimate the unknown parameters of the two ForClim regeneration models and their uncertainties. We estimated three species-specific (kL_y,s, kDDMin_y,s, and kDrTol_y,s) and two general (kEstDens, kTrMax) regeneration parameters of ForClim based on recruitment data from the EuFoRIa reserves. The species parameters were estimated for 11 out of 30 simulated species for which the data covered sufficient environmental variation. The species not considered for calibration were simulated with their default parameters (see Huber et al., 2020).

2.4 Calibration target

The calibration target was to obtain recruitment rates in the model that match to observations. Tree recruitment was quantified as the number of trees that pass an inventory-specific DBH threshold. Observed decadal tree recruitment rates R_i,s were calculated for each plot i and species s with

where p is the inventory period and T_i is the total number of years between the first and the last inventory at plot i. Simulated tree recruitment rates $\widehat{R}$ were calculated as

where ${\mathrm{nTrs}}_{i,s,p,k,j}$ is the number of recruited trees for one patch j, inventory period p, and repetition k. Each simulation during the calibration was conducted on > 100 patches (i.e., ca. 8 ha) to reduce the variability caused by the k stochastic realizations of the ForClim model. The trees in the initial forest inventory were randomly distributed to each of the 100 patches (each with a size of 0.08 ha) proportionally to the actual plot size until a full repetition exceeded 100 patches. If one repetition was not a multiple of the patch size of 0.08 ha, the difference in exceeded plot area determined the proportion of additional trees drawn from all trees in the initial forest inventory to populate the patches. The number of repetitions Nrep_i for each plot i emerges from the next integer above $\frac{\mathrm{8}\mathrm{ha}}{A_i}$, where A_i is the plot size in ha. This resulted in an average of eight repetitions k across all sites, although k ranged from 40 for plot sizes of 0.2 ha to 2 for two sites with plot sizes of > 4 ha. The number of patches j (${\mathrm{Npatch}}_{i,p,k}$) within one repetition k is the next integer above $\frac{\mathrm{100}\text{patches}}{{\mathrm{Nrep}}_i}$.

2.5 Model initialization

To initiate the calibration runs, we had to resolve the issue that trees below the inventory-specific DBH threshold (i.e., small trees), which may have been present in reality, are obviously not present in the data. Ignoring these trees in the initialization would create a temporal lag in tree recruitment (i.e., trees surpassing the DBH threshold), which is connected with a potentially significant underrepresentation of tree regeneration directly after model initialization. To overcome this problem, we initialized unobserved trees below the DBH threshold with the model's steady state (i.e., equilibrium of regeneration).

This steady state was determined by running the simulation with the stand structure of the initial forest inventory for 50 years and suspending all processes affecting trees above or equal to the DBH threshold. During this “spin-up” phase, trees above the DBH threshold that were included in the initial forest inventory could neither die nor grow, but they still modulated the variables of stand structure that affect regeneration (i.e., gAL₀, Trs). In contrast, trees below the DBH threshold were allowed to grow and die under the conditions observed in the initial inventory. If these newly regenerated trees grew larger than the DBH threshold during the spin-up, they were removed. This means that these trees did not die from mechanisms that simulate tree mortality in the model; they were forcefully removed from the simulation to avoid the accumulation of trees with a DBH close to the DBH threshold. Visual inspection of the simulation results showed that an equilibrium of the stand structure below the DBH threshold was reached after approximately 50 simulation years, which was the reason to fix the spin-up at this time period. After the spin-up, the simulation was continued with all trees and by running all model processes (i.e., regeneration, growth, and mortality).

2.6 Definition of goodness of fit

The goodness of fit was quantified by a (pseudo-)likelihood. We assumed that the observations R_i,s for site i and species s with the parameter vector θ were negatively binomially distributed, leading to a log-likelihood per observation of

Here, the mean ${\widehat{R}}_{i,s}$ is predicted by ForClim based on the chosen values for the model parameters kL_y, kDDMin_y, and kDrTol_y, the respective regeneration intensity parameters (kEstDens and kTrMax) of the two models (as explained above), and the dispersion parameter ϕ of the negative binomial distribution, which can be interpreted as a measure of residual variation and is a free parameter that needs to be estimated in addition to the model parameters.

We assumed that the dispersion may vary with species, DBH, and plot size according to the formula

where ϕ_s is the species-specific dispersion, ϕ_DBH is the effect of the diameter threshold DBH_i, and ϕ_A is the effect of the plot size A_i on the dispersion. We used an exponential function to only allow for positive values, as required by the negative binomial distribution. In sum, the likelihood depends on a parameter vector θ that includes the model parameters as well as the parameters that control the value of ϕ based on the DBH and plot size.

Summing the expressions in Eq. (10) over all plots i and over all species s, we arrive at the joint log-likelihood

We then additionally re-scaled the joint log-likelihood (Eq. 13) by a factor of $\mathrm{1}/\mathrm{12}$ (reflecting Nspecies) to arrive at a pseudo-likelihood that we used as the calibration target in our MCMCs.

The effects of this re-scaling are that the uncertainties get larger and the absolute magnitude of (stochastic) likelihood differences decreases. For both reasons, MCMC samplers are less prone to get stuck in local optima, whereas the shape of the posterior surface as well as the posterior optimum remain unchanged. However, the re-scaling (which essentially corresponds to a down-weighting of the observational evidence by a factor of 12) can be interpreted as accounting for possible non-independence of the data and structural model error (see Oberpriller et al., 2021). The value $\mathrm{1}/\mathrm{12}$ was chosen ad hoc, but, given that some structural error and non-independence are likely present in the observations, we believe it is of the right order of magnitude. Given the approximate nature of this correction, our scaled likelihood should be viewed as a pseudo-likelihood or an informal likelihood (Smith et al., 2008), and the same labels should be applied to the posterior.

In total, our pseudo-likelihood depends on a vector θ consisting of 48 parameters, as we estimated three ecological threshold parameters (kL_y, kDDMin_y, and kDrTol_y; see Fig. 1) and a separate dispersion parameter ϕ_s for each of the 11 species, one extra dispersion parameter for all other species, two parameters for the effects of DBH and plot area (ϕ_DBH and ϕ_A), and one parameter for the regeneration intensity (kEstDens or kTrMax for the simple or the complex model, respectively). We defined wide uniform priors for each parameter that comprise the full range from the species' lowest value to its highest value in the default parameterization in ForClim (see Table 1).

2.7 Posterior estimation

We calibrated the model using the differential evolution sampler (DEzs; ter Braak and Vrugt, 2008) implemented in the R package BayesianTools (Hartig et al., 2019). We sampled with two independent sets of three chains (i.e., a total of six chains) for all 353 training plots. The z matrix was re-initialized at the beginning of the sampling procedure (at 5000 to 6800 and 2000 to 3700 iterations for the simple model and the complex model, respectively) to improve the mixing of the chains. This was necessary because of the very wide prior range for the dispersion parameters, which led to a degenerate z matrix. The same procedure was applied to improve chain mixing after 120 000 to 139 300 (simple model) and 120 000 to 145 500 iterations (complex model). For the simple model, the upper prior range for the parameter kEstDens had to be adjusted from 0.02 to 0.2 after 139 300 iterations. Ultimately, after 191 600 (simple model) and 200 900 iterations (complex model), one set of three independent chains converged, as judged from a visual inspection of the chains and Gelman and Rubin's MCMC convergence diagnostic (see Table A4). The other set of independent chains did not fully converge, mostly because one chain was stuck for the species-specific dispersion parameters. Computational constraints did not allow us to run the sampler for even longer. A single simulation took 3 s per plot; i.e., 6 chains × 353 calibration plots resulted in a total computation time of 1.765 h per iteration without parallelization. Fortunately, the Euler High Performance Computing Cluster of ETH Zürich enabled us to use 1000 cores (500 per model variant) and sufficient memory. The effective computing time, including the overhead when utilizing all resources, was 10–15 s per iteration and ca. 25 d in total for each model variant.

The posterior distribution from the calibration consisted of 1000 samples drawn from the last 32 300 (simple model) and 45 400 (complex model) iterations. The simulations from the posterior parameter distribution provided posterior estimates of decadal tree recruitment rates (${\widehat{R}}_{i,s}$) for all 343 test plots i and species s. The mean posterior estimate (MPE) and the 80 % credible intervals (CIs) of ${\widehat{R}}_{i,s}$ from the 1000 posterior simulations were used to assess residuals and evaluate model performance (see Table A4). The MPE and CIs for the parameter estimate of the posterior distribution were used to compare the trait-based model and the calibrated model.

2.8 Performance comparison using the RMSE and marginal pseudo-likelihood (ML)

Model performance was assessed with the RMSE and the marginal pseudo-likelihood (ML) on both training and hold-out data. The difference between the two metrics is that the RMSE is a general metric of fit, while the marginal pseudo-likelihood is a Bayesian metric that relates to the Bayes factor and posterior model weights and thus allows the support for two alternative models to be compared based on a specific likelihood.

For the simulations from the calibrated and trait-based models, the RMSE was calculated for different DBH thresholds: for variable thresholds between plots and harmonized DBH thresholds of 7 and 10 cm. This harmonization was done by artificially increasing the DBH threshold in the observed and simulated data to mimic a consistent inventory design with a common DBH threshold. We calculated the RMSE from the training and test data based on a comparison of observed and simulated recruitment in terms of the species-unspecific, the species-specific, and the average over the species-unspecific and species-specific RMSEs (Table 2).

The ML was calculated only for the simulations of the calibrated models because the pseudo-likelihood relies on the dispersion parameters, which were not estimated for the trait-based model. The ML is the average pseudo-likelihood of the model given the training or the test data and is averaged over the posterior parameter uncertainty (see Delpierre et al., 2019). We evaluated the ML in both cases on the validation data based on the posterior distribution inferred from the training data. This approach, which corresponds to the fractional Bayes factor (O'Hagan, 1995), avoids inconsistencies when comparing models with weak or uninformative priors. The Bayes factor is then obtained by taking the ratio between two marginal likelihoods with e^M1−M2. This provides the relative posterior support of M1 relative to M2 by the data (Kass and Raftery, 1995).

When interpreting the results, it is important to remember that both the RMSE and the ML, as evaluated here, will typically be higher for more complex models of the training data, so the comparison of models with these metrics of the training data is of limited use. However, models can be sensibly compared by comparing their performances for the hold out, and it is also informative to look at the reduction in performance between training and the hold out, which gives an indication of overfitting.

2.9 Posterior sensitivity

Global sensitivity analyses allow for the assessment of model behavior across large parameter spaces. However, large parameter spaces may also cover unrealistic parameter configurations, and computational requirements are high. Therefore, a strategy for constraining the parameter space to a relevant location is required (see Huber et al., 2018). We combined the benefits of a global and a local sensitivity analysis by constraining the parameter space by deriving a posterior distribution from the observations. This allowed us to evaluate model sensitivity with respect to the uncertainty derived from observed recruitment patterns in European forest reserves.

To analyze the sensitivity of regeneration to changes in the model parameters within the posterior distribution, we analyzed the effect of increased tolerance of trait values (i.e., a lower kL_y, higher kDrTol_y, and lower kDDMin_y) on simulated recruitment within the posterior parameter range. This was done by modeling $\widehat{R}$ with a GLM and a negative binomial distribution with z-scaled values of negative kL_y, negative kDDMin_y, and positive kDrTol_y as predictors. This model was implemented using glmmTMB (Brooks et al., 2017).

3.1 Species traits and regeneration intensity

The trait-based regeneration niche differs from the regeneration niche that emerges from the model which was calibrated with the observations made in unmanaged forest reserves. Variation in these differences is evident between trait types, species, and model variants (Fig. 2 and Table A2).

Light requirements stand out from the other traits because they were most sensitive during model calibration, as indicated by the narrow posterior distributions compared to the prior parameter range (Fig. 2, left; Fig. A2a). Most posterior estimates were not only narrow but also supported by the trait values (Fig. 3, left). Estimates of the complex model were generally closer to the trait values, with a good rank correlation (i.e., Spearman's rho = 0.57). Light requirements defined by traits for the simple model were systematically lower compared to the values emerging from the calibration. Nevertheless, a Spearman's rho of 0.94 indicates that the calibration put the species in a plausible order (Fig. 3, top left). The estimates of the species-specific light requirements were more similar for both approaches for shade-tolerant tree species in general. Values from the calibration for Quercus spp., Pinus sylvestris, and Betula spp. were much lower compared to the trait values in the complex model but were more similar in the simple model. For Tilia cordata, the calibration led to much higher light requirements compared to the trait-based values with either model (see Fig. 2, left). Drought tolerance values from both approaches matched moderately well when using the complex model (Spearman's rho = 0.46) with a relatively large CI, and the expectation was within the CI for eight species (Fig. 2, center). The MPE from the calibration was close to the definitions of the trait values for five species. However, drought tolerance estimates from the calibration did not match the trait values well for the simple model (Spearman's rho = −0.18). Only for six species did the wide CIs include the trait-based values, and the MPE matched the trait values in the simple model for two species only. For Quercus spp. (complex model), Tilia cordata (simple model), and Pinus sylvestris (both models), the calibration led to much lower values compared to the trait values. Conversely, for Alnus glutinosa and Fraxinus excelsior, the calibration resulted in higher values compared to the trait-based values (see Fig. 2, center).

Estimates for the minimum degree-days had a wide CI for both models and a low rank correlation between the approaches (Spearman's rho = 0.19 for the complex model and −0.3 for the simple model; Fig. 2, right). Although the posterior CI of the calibration included the trait-based values for many species, the intervals were wide and the MPE values were close to the trait-based values for only two and three species in the complex and simple model, respectively. This indicates that the calibration values neither fully disagree nor perfectly match the trait-based values. For Tilia cordata (both models) and Carpinus betulus (simple model), the calibrated values were much lower compared to the trait-based values. Conversely, the calibrated degree-day values for Pinus sylvestris (simple model) and Abies alba (both models) were considerably lower compared to the trait-based values (see Fig. 2, right).

Figure 2Mean posterior estimate, including the 80 % CI, of the species-specific parameters for the complex model (red) and the simple model (blue). Point type indicates the approach (inverse calibration approach (ICA) or trait-based approach (TBA)). The panels show values for three species traits: light requirements (kL_y), drought tolerance (kDrTol_y), and minimum degree-days (kDDMinkL_y). The uniform prior parameter range (min, max) of each species trait (x axis) is indicated by the gray rectangle in the background.

Download

Figure 3Comparison of expected values for species-specific establishment thresholds based on ecological knowledge (TBA, trait-based approach; see Leuschner and Ellenberg, 2017) and MPE (ICA, inverse calibration approach; this study). The range displays the 80 % CI for the complex model (top panels) and the simple model (bottom panels). The 1:1 relationship is indicated by the black line. Spearman’s rank correlation of the MPE and the trait-based approach (TBA) values are shown in each panel.

Download

The estimates of the calibration for general regeneration intensity were narrow compared to the prior (Fig. 4) and therefore considerably sensitive. While the species-unspecific parameter kTrMax from the calibration (complex model) was significantly lower than its default value for the trait-based model, the species-specific parameter kEstDens (simple) was higher than its default value for the trait-based model (Table A1). Considering the interaction between kEstP (which was reduced by a factor of $\mathrm{1}/\mathrm{5}$) and kEstDens or kTrMax (see the model description in Appendix B), the overall amount of regeneration is generally lower for the calibrated model compared to the trait-based model. Specifically, for the complex model, the MPE (kTrMax = 8762) suggests 25 times less maximum regeneration compared to the default values of the trait-based models (kTrMax = 50 000 × 5 = 250 000). For the simple model, the MPA (kEstDens = 0.022) suggests a reduction of ca. 24 % per species compared to the default trait-based model (kEstDens = 0.006 × 5 = 0.03). However, it is noteworthy that these parameters modulate regeneration differently and that the magnitude of the deviation between the parameters in the calibrated model, the default trait-based model, and the model variants is not directly translated into the simulated regeneration amount within the model (see Fig. 1 and the full set of equations in Appendix B).

The coefficient for the effect of the DBH threshold in the species-specific dispersion was significantly different from zero, with an MPE of −0.39 and −0.41 for the complex model and simple model, respectively. This indicates that dispersion increases for higher DBH thresholds (Fig. 4c). The coefficient for the effect of plot size is slightly negative but not significantly different from zero (Fig. 4d), which indicates that there is no significant effect of plot size on dispersion. However, a very weak positive effect of larger plots on dispersion is visible. At the species level, dispersion effects differed significantly, with the lowest dispersion seen for Fagus sylvatica and the highest for Quercus spp. (Fig. A1 and Table A3). These findings emerged from both the simple and the complex model (see Fig. 4c and d).

Figure 4Posterior distributions of the non-species-specific parameters determining the amount of regeneration: (a) kTrMax (complex model, red) and (b) kEstDens (simple model, blue). Effects on the dispersion parameter ϕ: (c) DBH threshold and (d) plot size. The prior parameter in (a) and (b) is given by the extent of the graph. The prior range for the dispersion parameters was −5 to 5 and is not shown. Note that lower values of the dispersion parameter indicate higher dispersion. Consequently, negative estimates for dispersion are positive effects on the actual dispersion. Species' dispersion parameters are presented in Table A3.

Download

3.2 Model performance

The calibration led to somewhat better performance compared to the trait-based approach (Fig. 5 and Table 2). Both model variants performed better when calibrated and revealed the uncertainty of the posterior simulations. However, performance differed strongly between species. The most gains in performance coupled with a high degree of uncertainty were evident for Abies alba and Tilia cordata with both models (Fig. 5). No increase in performance was evident for Quercus spp. and Alnus glutinosa, although high uncertainty of the simulations was evident for Alnus glutinosa with the complex model. Slight but distinct gains in performance were found for Fagus sylvatica and Picea abies with the complex model. Note that not only the intercept for the comparison of observations and simulations but also the slope changed, which indicates that this was not only due to the regeneration intensity parameters. Overall performance increased distinctly for the complex model and considerably for the simple model (see the RMSE in Table 2). In summary, the calibration clearly improved model performance for almost all species in the complex model and to a limited extent in the simple model. In addition, the uncertainty based on the posterior parameter distribution was clearly visible in the simulations. The RMSE decreased with increasing DBH threshold (Fig. A3) for both models.

The ML confirmed the higher performance of the complex model (ML = −391.61) compared to the simple model (ML = −401.86), which was also evident from the RMSE (Table 2 and Fig. A2b). According to Kass and Raftery (1995), the fact that the Bayes factor from the training data is above 200 suggests strong support for the complex model (Bayes factor = 28 282.54).

Table 2RMSEs of the simulated and predicted R and the pseudo-likelihoods for the different models and approaches (the inverse calibration approach (ICA) and the trait-based approach (TBA)). The marginal pseudo-likelihood (ML) was derived from the posterior parameter distribution of the training data (N=353) and the test data (N=343). Species-specific RMSE values are presented in Fig. A4.

Download Print Version | Download XLSX

Figure 5Simulated vs. observed recruitment rates at 7 cm DBH of the 11 species for (a) the simple model and (b) the complex model. The red lines show linear regressions from 1000 simulated recruitment rates using the posterior parameter distribution at the 343 test plots. The data points used for the regression are indicated by the blue color of the hexagons, where light blue indicates fewer points and dark blue indicates more points. The yellow lines characterize the recruitment rates for the same test plots based on the default parameter setting (TBA, trait-based approach).

Download

3.3 Posterior sensitivity

The variation of the parameter values characterizing the species' light requirements had the strongest effect on conspecific regeneration for each model variant (Fig. 6a and b, respectively). The variation of drought tolerance within the posterior, which was rather high, had a much weaker effect on conspecific recruitment, and the minimum degree-days had a very low effect on recruitment within the posterior parameter range. Interestingly, positive effects for heterospecific recruitment were evident in the case of higher shade tolerance of Quercus spp. and Acer pseudoplatanus, but only in the complex model (Fig. 6b). Picea abies (simple model) and Fagus sylvatica (complex model) showed the strongest negative effects on the emergence of heterospecific recruitment (Fig. 6).

Figure 6Estimates of the effects of increased tolerances for kL_y, kDrToly_y, and kDDMiny_y on the emergence of conspecific and heterospecific recruitment for the simple model (a) and the complex model (b). The effect sizes correspond to the coefficients from a GLM that predicts recruitment with scaled tolerance parameters so that they have a mean of 0 and a standard deviation of 2; an increase in the value always indicates an improvement for the species.

Download

4.1 Species traits

The species trait values varied considerably between the calibrated models, the trait based-models, and the two model variants. We aim to explain these differences by reflecting on (a) the theoretical expectations of the two approaches and (b) the structural differences between the model variants.

Differences between the trait-based and calibrated models can be expected based on ecological considerations regarding the regeneration niche (Grubb, 1977) and ontogenetic niche shifts (Werner and Gilliam, 1984) as well as methodological aspects such as the importance of context for modeling trait–demography relationships (Yang et al., 2018). Specifically, we expected high sensitivity and a good match for the shade tolerance of the species because light availability is a key determinant of tree regeneration on small spatial scales (Collins and Good, 1987) and its context is modeled explicitly and in rather high detail in ForClim (note the direct link between stand structure and light availability in Eqs. 1 to 3). In contrast, the context of trait values related to climate (drought tolerance and temperature requirements) is only vaguely defined by species distribution limits (Meusel et al., 1965) along macroclimatic gradients (Rudloff, 1981). In addition, the traits for light requirements are differentiated between juveniles and adults (Leuschner and Ellenberg, 2017), while those related to climate are not. Therefore, we expected less agreement between the calibrated and trait-based models for the latter drivers. Our results supported these expectations, as shown by the mostly good agreement between the light requirement parameters compared to climatic parameters in the two approaches. Notably, this pattern was found in both the simple and the complex model, and it thus appears to be a robust feature, irrespective of the structure of the regeneration model. Furthermore, the drought-related parameters matched better between the calibrated and the trait-based models in the complex variant, which indicates that trait values embedded in a model that connects drought effects and competition during regeneration have more support from the data. This interlink of competition and drought has also been demonstrated in grassland communities (Grant et al., 2014; Levine et al., 2022) and tree species mixtures (Jucker et al., 2014; Grossiord, 2020; Young et al., 2017; Clark et al., 2016; Ruiz-Benito et al., 2013; Haberstroh and Werner, 2022).

4.2 Model performance

Generally, the calibration resulted in much better performance for the complex model and a moderate improvement for the simple model compared to the trait-based approach. This is consistent with previous studies using Bayesian calibration of dynamic models (Augustynczik et al., 2017; Cailleret et al., 2019; Trotsiuk et al., 2020; Van Oijen et al., 2005). The main reason for this improvement in both model variants is the overall lower regeneration amount, which results from the combination of establishment probability and regeneration intensity. Thus, our results suggest that calibration can help to sharpen the estimates of regeneration parameters that are not well constrained by standard empirical data.

The somewhat better performance of the complex model is best explained by the way the species-specific amount of regeneration is determined. While the simple model does this uniformly, the complex model distributes the regeneration to the species according to their environmental suitability. This is also reflected in the more realistic estimates of the regeneration niche along the drought gradient. Overall, our findings corroborate the considerations of Huber et al. (2020), who suggested the simultaneous use of different model variants. In our study, the regeneration patterns across the very heterogeneous forest types in our dataset were captured much better by the complex model, which implicitly allows for differentiating processes (captured via the EFs) in the regeneration layer. In contrast, the ideas underlying the simple regeneration model, which was originally developed for multi-species forests with high evenness (Botkin et al., 1972), turned out to be less suitable for reproducing the observed regeneration patterns.

In addition to established performance measures such as the RMSE or Bayes factor, the comparison of the default trait values and inversely calibrated trait values allowed us to evaluate whether the calibrated parameters are just “degrees of freedom” that are used to make the model fit better to the data or whether their estimates are plausible from an ecological perspective (see Hellegers et al., 2020). Based on the discussion of species trait values above, this leads to the conclusion that the simple model is more realistic for the factor of light, whereas the complex model captures processes related to drought better.

4.3 Methodological considerations