Understanding uncertainties and sensitivities of projected ecosystem dynamics under environmental change is of immense value for research and climate change policy. Here, we analyze sensitivities (change in model outputs per unit change in inputs) and uncertainties (changes in model outputs scaled to uncertainty in inputs) of vegetation dynamics under climate change, projected by a state-of-the-art dynamic vegetation model (LPJ-GUESS v4.0) across European forests (the species Picea abies, Fagus sylvatica and Pinus sylvestris), considering uncertainties of both model parameters and environmental drivers. We find that projected forest carbon fluxes are most sensitive to photosynthesis-, water-, and mortality-related parameters, while predictive uncertainties are dominantly induced by environmental drivers and parameters related to water and mortality. The importance of environmental drivers for predictive uncertainty increases with increasing temperature. Moreover, most of the interactions of model inputs (environmental drivers and parameters) are between environmental drivers themselves or between parameters and environmental drivers. In conclusion, our study highlights the importance of environmental drivers not only as contributors to predictive uncertainty in their own right but also as modifiers of sensitivities and thus uncertainties in other ecosystem processes. Reducing uncertainty in mortality-related processes and accounting for environmental influence on processes should therefore be a focus in further model development.
Terrestrial ecosystem models have emerged in the last 3 decades as a central tool for decision making and basic research on vegetation ecosystems (Cramer et al., 2001; Fisher et al., 2018; IPCC, 2014; Smith et al., 2001; Snell et al., 2014). Projections from different vegetation models, however, often disagree on important details, for example regarding the observable past (Bastos et al., 2020) or the future carbon uptake of forest ecosystems (Huntzinger et al., 2017; Krause et al., 2019). Among the possible reasons for such differences is the uncertainty in climate scenarios (Saraiva et al., 2019), model structural uncertainty (Bugmann et al., 2019; Oberpriller et al., 2021a; Prestele et al., 2016), initial condition uncertainty (Dietze, 2017b), and uncertainty about the model parametrization (Grimm, 2005), which in turn make models' projections themselves uncertain (Dietze, 2017a). It is widely appreciated that understanding which exact factors drive these uncertainties is of immense value for directing research (Tomlin, 2013) and also for interpreting and understanding projections (Dietze et al., 2018). For example, the IPCC started in its Fifth Assessment Report to systematically analyze uncertainties and attribute them to model inputs (IPCC, 2014) similar to other predictive sciences (e.g., nuclear reactor safety (Chauliac et al., 2011), energy assessment for buildings (Tian et al., 2018), or policy analysis (Maxim and van der Sluijs, 2011).
The two main tools to understand how uncertainties in model inputs (drivers, parameters, and model structure) affect model outputs are sensitivity analysis (SA) and uncertainty analysis (UA) (Cariboni et al., 2007; Caswell, 2019; Saltelli, 2002; Saltelli et al., 2008). The key difference between these two methods is that in a UA the central starting point is the quantification of uncertainty in the model inputs (e.g., parameters, typically determined via expert elicitations and previous studies Matott et al., 2009). This uncertainty is then propagated to the model outputs and back-attributed to the different inputs. An SA, on the other hand, calculates how the model output changes per unit or percentage of change in the respective input (Jørgensen and Bendoricchio, 2001). This calculation is primarily independent of the inputs' uncertainties, although local SAs can be affected by the reference point, and global SAs can be affected by the range over which the sensitivity is calculated. Overall, however, both methods share the goal of identifying inputs with a high influence on model outputs, with the underlying idea that better constraining these will increase robustness and reliability of model projections (Balaman, 2019).
Although the benefits for understanding model behavior and predictive uncertainties are obvious, relatively few SAs and UAs have been applied to complex ecosystem models and especially the widely used dynamic global vegetation models (DGVMs) that project terrestrial ecosystem responses to climate change or land management (see, e.g., Courbaud et al., 2015; Cui et al., 2019; Huber et al., 2018; Reyer et al., 2016; Tian et al., 2014; Wang et al., 2013). A reason for this is arguably the complex structure of most DGVMs (Fer et al., 2018), which makes SAs and UAs computationally demanding and difficult to interpret, especially when performing state-of-the-art global SAs and UAs that compute sensitivities and uncertainties across the entire parameter space (Saltelli et al., 2008) rather than just locally around a reference parameter set (see, e.g., Hamby, 1994). Moreover, several studies highlight that sensitivities and uncertainties of DGVMs also exist with respect to environmental drivers (Barman et al., 2014; Wu et al., 2017, 2018), especially solar radiation (Barman et al., 2014; Wu et al., 2018), temperature (Barman et al., 2014), and precipitation (Wu et al., 2017), and it is reasonable to expect that there can be interactions between parameter and environmental sensitivities, meaning that certain parameters are more sensitive in some environments than in others. It therefore seems important to investigate parametric sensitivities in conjunction with their environmental sensitivities in one combined analysis.
In this study, we concentrate on a well-established and widely applied DGVM, the Lund–Potsdam–Jena General Ecosystem Simulator (LPJ-GUESS) (Gerten et al., 2004; Sitch et al., 2003; Smith et al., 2001). Three previous SAs or UAs for the LPJ family identified the intrinsic quantum efficiency of CO
Since the publication of these studies, however, the structure of the LPJ-GUESS model changed substantially. The most important changes are the inclusion of the nitrogen cycle (Smith et al., 2014) and new management modules (Lindeskog et al., 2021). Since these changes, no study has systematically examined how model sensitivities and uncertainties were affected by the new model structure. Moreover, previous SAs and UAs ignored management parameters, which are expected to have large impacts on carbon pools and fluxes (Lindeskog et al., 2021).
A further limitation of most previous studies for LPJ-GUESS and other models (e.g., Mäkelä et al., 2020) is that they either analyzed sensitivities and uncertainties to parameter changes or changes in the environmental drivers (but not both). As discussed earlier, however, there are good reasons to expect that the sensitivity of parameters will change if environmental drivers change. Given that previous sensitivity analyses used different choices for these boundary conditions (different sensitivities for the climate scenarios and sites in Jiang et al., 2012; for different elevations in Pappas et al., 2013; different sites in Wramneby et al., 2008), this not only limits the comparability between studies but also questions the generality of the results for all climatic conditions. Only Jiang et al. (2012) combined parameter and driver sensitivities but used for the latter only a number of fixed climate scenarios instead of a range of possible values, which prohibits a systematic joint analysis. Moreover, it would be interesting to compare the relative importance of drivers and parameters for the predictive uncertainty of model simulations and how these change between environmental zones (here we use the classification of Metzger et al., 2005) and thus on an environmental gradient. When sensitivities or uncertainties of parameters belonging to a specific process increase on an environmental gradient, this indicates that the process itself becomes more important on the gradient (Saltelli, 2002). By comparing such changes to existing ecological hypotheses, we can test if model sensitivities and thus process descriptions are in line with ecological expectations.
To answer these questions, we analyzed sensitivities and uncertainties in LPJ-GUESS for 200 randomly distributed sites across Europe (see Appendix A1.1). We address the issue of interactions between environmental and parametric sensitivities by simultaneously investigating uncertainty in environmental drivers (precipitation, temperature, solar radiation, CO
LPJ-GUESS is a process-based ecosystem model that simulates vegetation growth, vegetation dynamics, biogeography, and biogeochemical (e.g., nitrogen and carbon) and water cycles (Lindeskog et al., 2013; Olin et al., 2015; Smith et al., 2014). Ecosystem dynamic processes in the model include establishment; growth; mortality; and competition for light, space, and soil resources. To simulate these processes, the model combines time steps on different scales from daily (e.g., phenological and photosynthesis processes) to yearly (e.g., allocation of net primary production to tree carbon components). LPJ-GUESS includes forest gap dynamics succession of cohorts (each represented by an average individual) of different plant functional types (PFTs) or species. Each PFT and species has a unique parameter set.
In this study, we use a model version that was slightly modified from Lindeskog et al. (2021), which is based on the LPJ-GUESS 4.0 version, with a re-parameterization for spruce (Picea abies), pine (Pinus sylvestris), and beech (Fagus sylvatica) (see Appendix A1.2 for Pin. syl. and Pic. abi.). To account for the stochastic components of establishment, mortality, and patch-destroying disturbances, LPJ-GUESS simulates several replicate patches (25 for the simulation with the reference parametrization and 1 for each simulation in the SA and UA) representing "snapshots" of the grid cell. In this model version, fire is based on the BLAZE model (Rabin et al., 2017). Thus, annually burned area is generated based on fire weather and fuel continuity and distributed to monthly intervals based on climatology (Giglio et al., 2010). Tree mortality is then estimated by computing fire lines based on weather and is converted into height-dependent survival probabilities (see Haverd et al., 2014) depending on empirical biome-specific parameters.
A first set of key parameters from our expert elicitation (see below) for establishment are the bioclimatic limits (i.e., minimum growing degree days (gdd5min_est), minimum 20-year coldest month (tcmin_est), maximum 20-year coldest month (tcmax_est), and minimum forest photoactive radiation at forest floor (parff_min)), which build the environmental envelope for establishment. Given that the bioclimatic limits are fulfilled, at regular intervals new PFTs are established (here 1 year) given enough space, light, soil, water, and photoactive radiation are all available for establishment at the forest floor (Smith et al., 2001). Moreover, each of our three investigated species has a maximum establishment rate (est_max) (Smith et al., 2001).
The structure of trees in the model is mainly linked to the simulated growth of trees, which is triggered by allocating all net primary production (NPP) besides a reproduction debt of 10 % (reprfrac) to tree components, thereby satisfying the mechanical balance, e.g., the allometric equation for the relationship between height and diameter with allometric parameters (k_allom2, k_allom3; e.g., Huang et al., 1992), the relationship between tree leaf to sapwood area (k_latosa; e.g., Robichaud and Methven, 1992), the relationship between crown area and height (k_rp or packing constraint; Zeide, 1993), the maximum crown area (crownarea_max), and leaf longevity (leaflong), and functional balance, as well as demographic constraints (Sitch et al., 2003). Each living tissue is assigned a turnover rate transferring sapwood into heartwood (turnover_sap) and leaves (turnover_leaf) and fine roots (turnover_root) to litter. Investment into aboveground and belowground growth is influenced by the resource stress as individuals are competing for light, space, nitrogen, and water. Competition for light is determined by the photosynthetic response and light extinction in the canopy. Competition for space (self-thinning) is represented in the model via allometric equations between crown area and stem diameter (Sitch et al., 2003). Competition for nitrogen and water is determined by individual tree demand for nitrogen and water, soil availability of nitrogen and water, and the PFT-specific root profile. Competition between species will favor certain life history strategies in particular situations, for example shade-tolerant (e.g., Fagus sylvatica and Picea abies) or intermediate-shade-tolerant (e.g., Pinus sylvestris) growth responses and dynamically changing root-to-shoot ratios.
Tree mortality (natural or via harvest) in the model responds to growth efficiency (ratio of annual NPP to leaf area) being too low over a 5-year period, e.g., due to light competition, maximum longevity of a PFT, or changes in environmental conditions (e.g., tolerance to drought (drought_tolerance) changes water uptake) exceeding the species suitable range. Light competition is modeled using the foliage projective cover (FPC), defined as the area of ground with foliage directly above it using Beer's Law (Smith et al., 2011). The resulting shading mortality is distributed proportional to species' FPC growth in the respective year due to their biomass increase. Mortality is modeled inversely proportional to the growth efficiency (with a given species-specific threshold (greff_min), e.g., Waring, 1983). Moreover, negative NPP of a species kills all individuals of the respective cohort. Background mortality probability increases with tree age, reaching 1 at the maximum longevity (longevity). Mortality has also a stochastic component. Natural disturbances are implemented in the model as process-based wildfires (with a given fire resistance for each species (fireresist)) and as patch-destroying disturbances (e.g., windthrow and landslides) with the same yearly occurrence probability for all patches (inverse of distinterval). Additional mortality arises from forest management activities, determined by thinning intensity (percentage of all trees cut, thinning_intensity) and cutting intervals (cut_interval), which can be set for each species individually. For a more detailed description of the management module and the additional management parameters, see Lindeskog et al. (2021).
Nitrogen input is implemented in the model through nitrogen deposition (prescribed) and biological nitrogen fixation. The latter is simulated empirically as a linear function with intercept (nfix_a) and slope (nfix_b) of the 5-year-averaged actual evapotranspiration (Cleveland et al., 1999). The resulting amount of nitrogen accumulates in the ecosystem equally over the year and directly adds to the available mineral soil nitrogen pool. When nitrogen is in living tissue, a fraction (nrelocfrac) is re-translocated before leaf and root shedding.
Photosynthesis is modeled as a function of absorbed photosynthetically active radiation, temperature (optimum temperature range for photosynthesis determined by pstemp_low and pstemp_high, Larcher, 1983), intercellular CO
Water availability for plants is based on precipitation and snowmelt in the two-layer soil hydrology submodule (for details, see Hickler et al., 2004; Smith et al., 2001). Vegetation transpiration and evaporation (with a maximum evapotranspiration rate (emax)) from bare ground and leaves both reduce water availability and runoff from saturated soil (Sitch et al., 2003). Water vapor exchange by the vegetation canopy is calculated on a daily basis within the photosynthesis scheme (e.g., minimum canopy conductance not associated with photosynthesis (gmin)). The water supply and transpiration demand are calculated on a daily basis and converted into a drought stress coefficient. Given this coefficient, the investment in roots at the costs of leaves is calculated.
We selected 200 study sites (see Appendix A1.1) spatially and environmentally stratified over Europe by applying random stratified sampling (using the R package splitstackshape, Mahto, 2019) with longitudinal and latitudinal coordinates and mean precipitation, solar radiation, and temperature as categories based on IPSL-CM5 Earth System Model CMIP5 (Dufresne et al., 2013) climate data. We chose 200 sites as a compromise between the high computational demand of running LPJ-GUESS multiple times for all sites and a good spatial and environmental coverage of Europe. For these sites, we performed simulations for each of the three most common species in Europe (Fagus sylvatica, Pinus sylvestris, and Picea abies) as monospecific stands and additionally all three species together as mixed stands.
The simulation period was from 1861 to 2199. To start the simulations with equilibrium carbon pools and fluxes, we spun up LPJ-GUESS vegetation and soil carbon and nitrogen pools to pre-industrial equilibrium by recycling the 1861 to 1900 climate, the 1861 CO
The a priori selection of the most influential parameters that can be specified in the parameter file and their ranges was based on our expert knowledge (following the SHELF expert elicitation protocol; see Gosling, 2018) and a literature review. The resulting 11 (33 %) parameters common for all species and 22 (20 %) species-specific parameters (see Table 1) were grouped to the specific processes they contribute most to (Table 1, Grouping).
Table 1 The model inputs investigated in the sensitivity analysis can be grouped into (a) common parameters, (b) species-specific parameters, and (c) drivers. The ranges for the parameters have been determined from experts and literature, and default parameter values that changed from Hickler et al. (2012) due to the re-parameterization are explained in Appendix A1.2.
Grouping Parameter Explanation Unit Default value Min. value Max. value – – – – – – Literature sources Mortality andmanagement distinterval Average return time for genericpatch-destroying disturbances year 920 200 1000 – – – – – – – Nitrogen nfix_a First term in N fixation equation – 0.102 0.102 0.367 – – – – – – – Nitrogen nfix_b Second term in N fixation equation – 0.524 0.524 – – – – – – Cleveland et al. (1999) Nitrogen nrelocfrac Fraction of N re-translocated prior toleaf and root shedding – 0.5 0.1 0.8 – – – – – – – Photosynthesis andlight lambda_max Non-water-stressed ratio ofintercellular to ambient CO pp – 0.8 0.6 0.8 – – – – – – Pappas et al. (2013) Structure andphenology reprfrac Fraction of NPP allocated toreproduction – 0.1 0.05 0.3 – – – – – – – Structure andphenology turnover_root Rate of fine root turnover 1 yr 0.7 0.65 0.75 – – – – – – – Structure andphenology crownarea_max Maximum crown area m 40 20 60 – – – – – – – Structure andphenology k_allom2 height kallom diameter – 60 30 80 – – – – – – – Structure andphenology k_rp crown area kallom height – 1.6 1.3 1.6 – – – – – – – Water emax Maximum evapotranspiration rate mm d 5 2 6 – – – – – – –
Table 1 Continued.
Group Parameter Explanation Unit Default Min. Max. Default Min. Max. Default Min. Max. Literature value value value value value value value value value sources – – – – – Establishment parff_min Min. forest floor PAR for grassgrowth and tree establishment J m d 2 500 000 1 500 000 3 500 000 1 000 000 750 000 1 600 000 1 000 000 750 000 1 600 000 – Establishment gdd5min_est Min. GDDs on 5 C base forestablishment C d 500 250 700 350 300 700 1300 1050 1450 – Establishment tcmin_est Min. 20-year coldest-month mean temp for establishment C – Establishment tcmax_est Max. 20-year coldest-month mean temp for establishment C 5.5 6 3 6 7 5 8 Schilbaslki et al. (2017) Establishment est_max Max. sapling establishment rate 1 m yr 0.2 0.1 0.25 0.1 0.05 0.2 0.2 0.05 0.25 – Establishment alphar Shape parameter for recruitment andgrowth rate relationship – 10 4 15 4 2 5 2 0.8 5 – Mortality andmanagement longevity Expected longevity under lifetime non-stressed conditions (yr) year 500 300 900 300 200 1000 400 250 650 – Mortality andmanagement fireresist Fire resistance – 0.4 0.05 0.7 0.1 0.05 0.8 0.1 0.05 0.8 – Mortality andmanagement cutinterval Time until trees are cut year 90 40 140 90 60 120 105 80 140 – Mortality andmanagement greff_min Threshold for growth suppression mortality kgC m yr 0.21 0.07 0.26 0.135 0.03 0.19 0.02 0.001 0.13 Pappas et al. (2013) Mortality andmanagement drought_tolerance Implements drought-limited establishment plus water uptake, from 0 (total) to 1 (not at all) – 0.25 0.1 0.4 0.48 0.2 0.65 0.39 0.2 0.49 – Mortality andmanagement thinning_intensity Percentage of threshold crown coverage that is kept after thinning – 0.9 0.45 1 0.9 0.5 1 0.9 0.55 1 – Photosynthesis andlight respocoeff Respiration coefficient – 1 0.8 2.2 1 0.8 2.2 1 0.5 1.5 – Photosynthesis andlight pstemp_low Approx. lower range of temp optimum for photosynthesis C 10 6.75 15 10 6.75 14 15 8 20 Renwick et al. (2019);Pallandt-Vermeulen et al.(2015) Photosynthesis andlight pstemp_high Approx. higher range of temperature optimum for photosynthesis (C) C 25 16 30 25 16 30 25 20 30 Zhang et al. (2014) Structure andphenology cton_leaf_min Minimum leaf ratio – 31.90 27.32 38.37 38.37 31.9 43.16 24.06 22.7 27.19 – Structure andphenology sla Specific leaf area m kgC 8.56 7.812 9.3 11.52 8.7 15.1 43.08 28.33 48.23 Mencuccini and Bonosi(2001); Pallandt-Vermeulen et al. (2015);Xiao et al. (2006) Structure andphenology turnover_sap Rate of sapwood turnover fraction peryear 0.085 0.05 0.1 0.065 0.04 0.09 0.085 0.05 0.1 – Structure andphenology k_latosa Tree leaf to sapwood area ratio – 3000 1800 5200 4000 2500 7000 5000 2500 8000 – Structure andphenology k_allom2 height kallom diameter – 30 15 60 Values from common parameters Values from common parameters Zhang et al. (2014) Structure andphenology k_chillb Coefficient in equation for budburstchilling time requirement – 100 80 800 100 80 800 600 250 800 – Water gmin Minimum canopy conductance not associated with photosynthesis mm s 0.3 0.22 0.38 0.3 0.22 0.38 0.5 0.42 0.58 Pappas et al. (2013)
Table 1 Continued.
Environmental drivers insol Mean deviation of solar radiation fromstandard scenario RCP 4.5 per site W m RCP 4.5 65.2 – – – – – – – Environmental drivers temp Mean deviation of temperature fromstandard scenario RCP 4.5 per site C RCP 4.5 5.82 – – – – – – – Environmental drivers prec Mean deviation of precipitation fromstandard scenario RCP 4.5 per site mm per month RCP 4.5 6.27 – – – – – – – Environmental drivers CO Mean deviation of CO from standard scenario RCP 4.5 per site ppm RCP 4.5 237 – – – – – – – Environmental drivers ndep Mean deviation of nitrogen depositionfrom standard scenario RCP 4.5 per site g mm yr RCP 4.5 – – – – – – –
From the environmental drivers of the model, we selected incoming solar radiation, temperature, precipitation, atmospheric CO
LPJ-GUESS predicts a substantial number of output variables, which could all be examined regarding their sensitivities and uncertainties. Here, we concentrate on carbon outputs (gross primary production, GPP; total standing biomass, TSB; and net biome productivity, NBP) because of forests' role for carbon cycling (Bonan, 2008), their large contribution to the land carbon sink (Pugh et al., 2019), and the economic importance of tree growth for forest owners (Pearce, 2001).
Sensitivities and uncertainties were calculated by Monte Carlo sampling from the assumed multivariate parameter and climate uncertainty. For the monospecific and mixed simulations, we randomly drew 10 000 and 50 000 parameter and climate combinations, respectively, randomly from the prespecified uncertainty ranges and ran the model based on these combinations for each of the 200 sites. For the mixed simulations we individually drew parameter combinations for each species for each simulation, meaning that the same parameter could be different for different species. In total, this means that 16 million (
We quantified sensitivity and uncertainty indices by running multiple linear regressions with the model output averaged over time as a response and parameters and drivers, as well as their second-order interactions, as predictors. With 200 sites, each having three monospecific and one mixed stand setup, we ran
It is important to note that uncertainties and sensitivities have different interpretations, and which of these two is more relevant strongly depends on the purpose. The calculated percentual sensitivities can be interpreted as percentage change in the corresponding output when changing a parameter value 1 % in the prespecified range. The calculated uncertainties per parameter or driver can be interpreted as relative proportion of the overall uncertainty budget coming from environmental drivers and parameters. For scenario analysis, e.g., comparing different cut intervals of forests, sensitivities provide a direct estimate of the model response, e.g., how much biomass changes when the cut interval is changed. For a comparison of different model forecasts, uncertainties are usually more relevant. If a reduction of uncertainty via a model–data comparison is the purpose, both measures are important, as parameters with high sensitivities can contribute more or less predictive uncertainty depending on their input uncertainty.
Regardless of the output variable, LPJ-GUESS was most sensitive to photosynthesis-related parameters (respcoeff, lambda_max), parameters controlling the wood turnover (turnover_sap) and tree allometry (k_rp), water-related parameters (emax), mortality-related parameters (greffmin), and environmental drivers (temperature, CO
Graph: Figure 1 Relative sensitivities (percent output change per percent parameter change) of the individual parameters and environmental drivers regarding (a) total standing biomass, (b) net biome productivity, and (c) gross primary production. Sensitivities were not substantially different between Fag. syl. (green squares), Pic. abi. (blue circles), and Pin. syl. (red triangles), but parameter sensitivities were stronger for monospecific stands than mixed stands (purple asterisks). The height of the bar reflects the mean over monospecific and mixed stands. Positive values for points and bars indicate a positive and negative values a negative relationship with the corresponding output.
Mixed stands were less sensitive to changes in parameters than monospecific stands (Fig. 1). For monospecific simulations, species sometimes showed different magnitudes and even directions of sensitivities. In comparison to the other species, Fag. syl. was more strongly affected by bioclimatic limits, while Pin. syl. showed higher sensitivity to environmental drivers (temperature and solar radiation). Moreover, TSB and GPP are negatively sensitive to temperature, except for Fag. syl. For NBP, the direction of sensitivities changes between species for the non-water-stressed ratio of intercellular to ambient CO
Looking at uncertainties, we found that environmental drivers contributed most of all processes and drivers to the predictive uncertainty (Fig. 2), regardless of the considered model output. For TSB projections, CO
Graph: Figure 2 Uncertainty contributions in percent of the individual parameters and environmental drivers regarding (a) total standing biomass, (b) net biome productivity, and (c) gross primary production showed no strong differences between Fag. syl. (green squares), Pic. abi. (blue circles), and Pin. syl. (red triangles) and were stronger for monospecific stands than mixed stands (purple asterisks). The height of the bars reflects the mean over monospecific and mixed stands. Positive values for points and bars indicate a positive and negative values a negative relationship with the corresponding output.
By analyzing uncertainty contributions on a species level, a more diverse picture emerged. Fag. syl. was more affected by temperature and less affected by solar radiation than the other species. Additionally, we found that uncertainty contributions of environmental drivers were substantially higher for mixed stands than for monospecific stands.
To project the uncertainties of TSB (for GPP and NBP, see Appendix 1.4) into the European environmental space, we filtered stands according to environmental zones and then calculated mean uncertainties per environmental zone and aggregated these per process.
The broad pattern of TSB uncertainty contributions for all three monospecific and mixed stands remains similar in all environmental zones. On average across all environmental zones, stands and species about 45 % of the uncertainty was due to environmental drivers, 15 % was due to mortality-related parameters, 14 % was due to photosynthesis-related parameters, 12 % was due to structure-related parameters, 7 % was due to water-related parameters, and 7 % was due to nitrogen-related parameters (Fig. 3).
Graph: Figure 3 The aggregated relative uncertainties of total standing biomass per environmental zone (with more than five sites) show a higher importance of drivers in the south than in the north. The environmental zones are from Metzger et al. (2005) and are defined as follows: ALN, Alpine north; ALS, Alpine south; ANA, Anatolian; ATC, Atlantic central; ATN, Atlantic north; BOR, Boreal; CON, Continental; LUS, Lusitanian; MDM, Mediterranean mountains; MDN, Mediterranean north; MDS, Mediterranean south; NEM, Nemoral; and PAN, Pannonian. In the radar plots of each environmental zone, the color and percentage value of the process label indicates which simulation setup (monospecific with corresponding species or mixed stands) has contributed most uncertainty and how much.
For the individual environmental zones, however, there were subtle differences. In the Mediterranean mountain (MDN) and Pannonian (PAN) zones, environmental-driver-induced uncertainty was higher than on average, especially for monospecific stands (Fig. 3). In the Boreal (BOR), Atlantic central (ATC), and Atlantic north (ATN) zones, tree-structure-related uncertainty increased compared to the average pattern (Fig. 3). In the Atlantic central (ATC) and Atlantic north (ATN) zones, nitrogen-related uncertainty increased for all species and stands (Fig. 3).
To examine this spatial pattern further, we investigated the change of uncertainties across a temperature gradient. To this end, we aggregated the uncertainties per site and process or driver and then fitted a linear regression with the process or driver as a predictor and the aggregated uncertainties as dependent variables.
For TSB, we found that increasing mean annual temperature increased the uncertainty contributions of environmental drivers and water and establishment parameters, while the uncertainty due to nitrogen-related and tree-structure-related parameters decreased (Fig. 4a). Thereby, the uncertainty contributions of environmental drivers (
Graph: Figure 4 The uncertainty contributions to total standing biomass projections of parameters and environmental drivers change across a mean annual temperature gradient across Europe from north to south (with p values and R2 for the processes and drivers). With increasing temperature, the importance of drivers and establishment became higher for total standing biomass, while the uncertainty contributions from nitrogen and structure declined (a). The uncertainty contributions due to temperature increased on the temperature gradient, and the contributions from solar radiation decreased (c).
Looking in more detail at the environmental drivers, uncertainty induced by temperature (
The above geographical and correlative observations of changing uncertainties across Europe receive further support when looking at the interactions between uncertainties of different drivers and parameters (Fig. 5). Interaction indices were calculated by averaging the interactions found in the linear regression over all sites and species (Fig. 5b). Moreover, to investigate the overall influence on other parameters or drivers we summed the absolute individual interaction indices of each parameter together (Fig. 5a).
Graph: Figure 5 The induced uncertainty of environmental drivers and mortality- and photosynthesis-related parameters changed the most depending on other parameters (a). Strong individual interactions between parameters and environmental drivers in monospecific projections of total standing biomass were rare (b). If strong interactions occurred, these were mainly between two environmental drivers or environmental drivers and parameters and only rarely between two parameters (b).
We found that environmental drivers (temperature, solar radiation, CO
In this study, we analyzed sensitivities and uncertainties of the LPJ-GUESS vegetation model due to environmental driver and parameter variations across European forests. We found that the model is most sensitive to relative (percentage) changes in photosynthesis-related parameters, structure-related parameters controlling the wood turnover and tree allometry, water-related parameters, mortality-related parameters, and environmental drivers (Fig. 1), irrespective of the considered output variable. When considering the different uncertainties (i.e., the entire plausible range) in these parameters and the environmental inputs, we found that environmental drivers and parameters controlling evapotranspiration, background mortality, and nitrogen cycling contribute most to predictive uncertainty (Fig. 2). When correlated against a temperature gradient and thus geographically from north to south, uncertainty contributions to TSB increased for environmental drivers and decreased for tree-structure- and nitrogen-related parameters (Figs. 3, 4). Interactions between the uncertainty contributions were mainly between different drivers or between model parameters and drivers, whereas only a few parameter–parameter interactions were present (Fig. 5).
Our finding that average sensitivities of carbon-related projections across European forests were highest for photosynthesis-related parameters amplifies the evidence from earlier studies (Pappas et al., 2013; Zaehle et al., 2005), although we have used different parameter ranges. In addition, the finding about high sensitivity of LPJ-GUESS to parameters controlling tree structure and especially carbon turnover (turnover_sap) (Fig. 1) is in line with results reported for a previous version of LPJ-GUESS (Pappas et al., 2013) and its important role for carbon allocation in trees found in empirical studies (e.g., Herrero de Aza et al., 2011). The finding that carbon-related projections are very sensitive to mortality-related parameters (greffmin) is also supported by previous studies on the sensitivity of vegetation models and underlines the importance of improving mortality submodules for generating precise projections of vegetation dynamics (Bugmann et al., 2019; Hardiman et al., 2011). Moreover, sensitivities in mixed stands were lower than in monospecific stands for NBP and GPP (Fig. 1) (in line with Wramneby et al., 2008). The reason for that imbalance may be that other species can dampen and even benefit from non-optimal life history strategies of an individual species (Loehle, 2000). Another reason might be that we sampled parameters for each species individually for mixed simulations, which reduces the influence of each parameter on stand-level carbon projections.
We found that uncertainty contributions of environmental drivers were comparable to the uncertainty contributions of all parameters together (Figs. 2–5; see also Snell et al., 2018, for the FLMs model; in addition, see Petter et al., 2020, who found that most uncertainty is induced by the choice of the forest model). Particularly high uncertainty contributions arose from temperature (negative effect for TSB, GPP positive for NBP), CO
Several of our results suggest that environmental context influences the sensitivity of LPJ-GUESS model parameters. First, we found changing uncertainties across different vegetation zones (Fig. 3) and on an environmental gradient (Fig. 4) and that most interactions occurred with environmental drivers (Fig. 5). Moreover, uncertainty contributions analyzed by a random forest were similar to the linear regression results but assign higher importance to environmental drivers (see Appendix A1.3). All of these findings indicate that environmental context can change the importance of different processes in the model, which is in line with the biological expectation that the environment affects the physiology of organisms directly, and thus indirectly affects fitness and biotic interactions (e.g., Seebacher and Franklin, 2012; Tylianakis et al., 2008), and the expectation that environmental responses can be particularly nonlinear (e.g., Burkett et al., 2005) or show higher-order interactions.
Interestingly, our results of decreased uncertainty contributions of structure-related parameters and increased contributions of environmental drivers on the temperature gradient (Fig. 4) also seem in line with the stress-gradient hypothesis (Maestre et al., 2009), an empirically observed pattern that states that in stressful environments, positive interactions should occur more often than in benign environments (e.g., Callaway, 2007). For the ecosystem that we consider, we interpret increasing temperature as increasing stress (e.g., Ruiz-Pérez and Vico, 2020) and structure as the best indicator for competitive interactions as the structure dictates resource allocation (e.g., bigger crown but identical stem diameter leads to more photosynthesis, or more sapwood to heartwood turnover requires less NPP). With this interpretation, one would conclude that under increasing stress the importance of competition-related parameters decreases in the model, as expected from the stress-gradient hypothesis. We acknowledge that a fair amount of interpretation is needed to arrive at this conclusion, and we do not claim that this result lends evidence to the empirical discussion about the generality of the stress-gradient hypothesis, but we find it noteworthy that such a large-scale pattern emerges in the model from lower-level processes without having been imposed (see also Levin, 1992).
The management and the nitrogen cycling module are the most recent improvements of the LPJ-GUESS model (Smith et al., 2014; Lindeskog et al., 2021). Compared to previous sensitivity and uncertainty analysis, the high contributions of the nitrogen fixation to the predictive uncertainty of TSB and GPP (Fig. 2a, c) are novel (albeit not surprising), as nitrogen is an important factor for the productivity of most temperate and boreal ecosystems (Vitousek and Howarth, 1991). The main reason why few earlier studies report those uncertainties is that vegetation models have only recently begun to integrate nitrogen cycling and limitation (e.g., Smith et al., 2014). The management module showed only small uncertainties, which could be due to the narrow parameter ranges for the cut interval and thinning intensity reflecting typical forest owners' choices. As forest owners usually try to maximize their profits (Johansson, 1986; Brazee and Amacher, 2000) and thus biomass production, the low sensitivity of the management module is not surprising. A more suitable and important test case and application of the management module would be a historical reconstruction of foliage projective cover data or similar outputs of the LPJ-GUESS model.
Our study helps to guide the model application, discussion of uncertainties, and model development of LPJ-GUESS and other DGVMs. First, future model applications and model comparisons should focus on mortality as this process contributes high uncertainties for carbon-related projections (see Figs. 1–3; see also Fisher et al., 2018). Thereby, it should be investigated if these uncertainties stem from the intra-specific variability of the parameters itself (Bolnick et al., 2011) or the parameters not being identifiable (see Marsili-Libelli et al., 2014) or if a model–data comparison could reduce uncertainties in the parameters (e.g., Hartig et al., 2011; Dietze, 2017b). Using time series inventory data might help as they are informative for constraining mortality modules (Cailleret et al., 2020). Second, lower sensitivities of establishment-related parameters are surprising as we know that not all three investigated species can effortlessly establish across all of Europe, e.g., Fag. syl. can only establish in locations with no extreme drought and heat and no extreme winter frosts (Bolte et al., 2007). Thus, either we missed important parameters of this module or the parametrization of the model needs to be updated. Third, when introducing new processes or coupling with other models (e.g., Forrest et al., 2020) calculating interactions helps to get a first impression as to how these new processes influence other model processes and to potentially detect missing links. Moreover, future model applications can interpret their results with regard to the sensitivities in different factors (Saltelli et al., 2019) and discuss uncertainties and the causing factors when used in policy advice (Laberge, 2013).
We caution that our results regarding the importance of different factors for predictive uncertainties (but not sensitivities) depend on the a priori-defined uncertainty range of the contributing factors (see Wallach and Genard, 1998), as well as on several other technical choices in our study. For determining uncertainty ranges of the drivers, we used RCP scenarios; however, these were not created as probabilistic minimum and maximum ranges. For the model parameters, we relied on expert guesses, reducing subjectivity as far as possible by following the SHELF expert elicitation protocol (Gosling, 2018). Future studies could include more experts and their opinion on parameter distributions to reduce variability in this protocol. As the model is sensitive to parameters and environmental drivers and these elements influence each other, we treated them in a combined sensitivity and uncertainty analysis (Saltelli et al., 2019); however, when interpreting it should be kept in mind that the one group relates to uncertainties in the model, while the other is external, meaning that the two are conceptually very different (see also Dietze, 2017b). A certain ambiguity also arises from the definition of the indicators: here, we calculated sensitivities and uncertainties by capturing only linear components and second-order interactions, and we may therefore miss highly nonlinear (and in particular hump-shaped) responses in LPJ-GUESS (Roux et al., 2021). However, our comparison to uncertainties calculated with random forest variable importance, a method that would also capture nonlinearities, did not reveal any qualitative differences in the ranking of parameter importance (Appendix A1.3). Overall, while we acknowledge that a certain amount of subjectivity exists in the choice of input uncertainty and calculation of indices, we believe that our results are quantitatively robust to those choices.
Moreover, we acknowledge that LPJ-GUESS is known to be sensitive to the scaling parameters alpha_a and alpha_C3 (Pappas et al., 2013; Zaehle et al., 2005), which we have omitted from our analysis. However, these parameters are not accessible in the parameter input file. Instead, they are hard coded in the model's source code, and therefore a normal user would not change them. We argue that these parameters should thus be counted towards the more general and here neglected contribution of structural uncertainty (i.e., the uncertainty regarding the functional form of processes or even to entire modules) to the joint model uncertainty. Several previous studies suggest that the sensitivity of vegetation models to structural changes can be large and often larger than that to parameters (e.g., Bugmann et al., 2019), and it would certainly be useful (although very complicated) to explore these uncertainties together with the factors considered here in a joint analysis. In the present study, however, we considered only the parameters that would be accessible to normal LPJ-GUESS users and neglect structural uncertainty that could be explored by changing the source code.
Our findings highlight the relative importance of parametric uncertainties in different processes and their interactions with uncertainties in environmental drivers for carbon projections with LPJ-GUESS. Our results demonstrate that environmental context changes uncertainty contributions of other processes across the European environmental gradient. The pattern of decreasing importance of competition towards the warmer areas is in line with the stress-gradient hypothesis, which posits that the importance of competition decreases with increasing environmental stress. Our findings improve our understanding of forest ecosystem models, enable pathways for future ecosystem model development, and thus build a basis for more realistic projections. In the future, parametric uncertainties could be reduced by model–data fusion (e.g., Trotsiuk et al., 2020) of LPJ-GUESS, concentrating on the parameters contributing most uncertainty in each geographic region (Fig. 3). Reducing uncertainties in the drivers is more difficult. To some extent, environmental drivers are themselves influenced by the vegetation (Strengers et al., 2010), and thus model–data fusion on a fully coupled model including feedback loops between vegetation and climate and a general improvement of climate models could reduce driver uncertainty to some degree. However, much of the uncertainty in this section effectively arises from potential greenhouse gas emission trajectories, for which a probabilistic assignment is difficult due to their dependency on human decision-making.
We sampled 200 geographically and environmentally stratified sites over Europe and thereby avoided sites near the sea. The corresponding sites with their average temperatures are shown in Fig. A1 and cover most European climates and vegetation zones.
Table A1 Differences in parametrization of Hickler et al. (2012) and our study for the investigated species (Fag. syl., Pic. abi. and Pin. syl).
Parameters Hickler et al. (2012) Our study Hickler et al. (2012) Our study Hickler et al. (2012) Our study drought_tolerance 0.3 0.3 0.43 0.48 0.25 0.25 fireresist 0.1 0.1 0.1 0.1 0.2 0.4 leaflong 0.5 0.5 4 7 2 4 turnover_leaf 1 1 0.33 0.1429 0.5 0.25 turnover_sap 0.085 0.085 0.05 0.065 0.065 0.085 est_max 0.05 0.1 0.05 0.1 0.2 0.2 alphar 3 10 2 4 6 10 parff_min 1 250 000 1 000 000 1 250 000 1 000 000 2 500 000 2 500 000 tcmin_surv (min. 20-year coldest-month mean temperature for survival) 3.5 7.5 30 30 30 30 tcmin_est (min. 20-year coldest-month mean temperature for establishment) 3.5 6.5 29 29 30 29 tcmax_est (max. 20-year coldest-month temperature for establishment) 6 7 1.5 3 1 5.5 twmin_est (min. warmest-month mean temperature for establishment) 5 1000 5 1000 5 8 k_chillb 600 600 100 100 100 100 sla 43? 43.08 11 11.52 8 8.56 k_allom2 40 60 40 60 40 60 wooddens 200 293 200 185 200 211 longevity 500 400 500 300 500 500 ga (aerodynamic conductance) 0.04 0.04 0.14 0.14 0.14 0.14 gdd5min_est 1500 1300 600 350 500 500
Graph: Figure A1 Our 200 sampled sites, which were geographically and environmentally stratified over Europe to cover the most important countries and climate and temperature zones.
There are several technical and methodological reasons requiring a re-parametrization of LPJ-GUESS for our study. First, most European forests are managed and species are planted far outside of their natural distribution. Second, the introduction of the nitrogen cycle (Smith et al., 2014) changed the model structure, and thus the parameters require an adjustment. Third, the productivity of trees in managed forests did not fit to the reported inventory data (Fig. A2). To account for all these issues, we adjusted the parametrization of (Hickler et al., 2012) to allow species growing according to their actual (i.e., caused by forest management) distribution instead of their natural distribution.
Picea abies and Pinus sylvestris are planted especially far outside their natural distribution (Fig. A2). Therefore, we specifically adjusted bioclimatic limits, drought tolerances, longevity, leaf turnover, disturbance intervals, and allometry for these species.
Graph: Figure A2 Simulated (black points), observed (blue), and natural distributions (green) of the adjusted parametrization (b, d) compared to applying the parametrization from Hickler et al. (2012) (a, c) for Picea abies and Pinus sylvestris. EUFO indicates data from EUFROGEN (2008 and 2013), Mauri indicates data from Mauri et al. (2017), and Caudullo indicates data from Caudullo (2017). The simulations were run from 1600 to 2010 without management and without competition between species. The plotted biomasses were averages over the last 20 years.
To check the consistency of the results obtained via linear regressions, we compare them to the variable importance of random forest. The variable importance measures additionally nonlinear effects, and thus it should be able to deal with nonlinear models like DGVMs. We calculated the variable importance the same way as we did for the linear regression by fitting a random forest with all parameters against the sum of differences between model outputs with default values and model outputs with sampled parameters. As our parameters were sampled from a uniform distribution with no correlation between the individual parameters, random forest variable importance can be compared to linear regression results.
The ranking is very similar to the ranking of the parameters and environmental drivers obtained via linear regression (Fig. A3). There is, however, a difference in the magnitude of the uncertainty induced by drivers, which is higher compared to linear regression (Fig. A3). The higher uncertainty due to drivers is thus a nonlinear effect and stresses our conclusion that environmental conditions change the uncertainty contributions of other parameters.
Graph: Figure A3 Results of the random forest uncertainty contributions. The uncertainties due to environmental drivers are higher than the uncertainties due to parameters compared to linear regression results, but the ranking of parameters is similar to linear regression results.
Interactions of gross primary production (Fig. A4a, b) and net biome production (Fig. A4c, d) are similar to the interactions of total standing biomass. These interactions are mostly between environmental drivers and environmental drivers or between environmental drivers and parameters (Fig. A4). Some strong interactions are between parameters and parameters, however, in such interactions there are always parameters included that have strong interactions with environmental drivers (Fig. A4).
Graph: Figure A4 Interactions of uncertainty contributions of GPP and total standing biomass are similar to net biome productivity, with most interactions arising from environmental drivers.
High sums of strong interactions arise from temperature, precipitation, solar radiation, greffmin, emax, and respcoeff (Fig. A4a, b).
We acknowledge funding from the Bavarian Ministry of Science and the Arts in the context of Bavarian Climate Research Network (bayklif). We thank the LPJ-GUESS developers for developing and maintaining the LPJ-GUESS model. We also thank two anonymous reviewers for their valuable comments and feedback on an earlier version of the manuscript.
By Johannes Oberpriller; Christine Herschlein; Peter Anthoni; Almut Arneth; Andreas Krause; Anja Rammig; Mats Lindeskog; Stefan Olin and Florian Hartig
Reported by Author; Author; Author; Author; Author; Author; Author; Author; Author