Impacts of trait variation through observed trait – climate relationships on performance of an Earth system model : a conceptual analysis

In many current dynamic global vegetation models (DGVMs), including those incorporated into Earth system models (ESMs), terrestrial vegetation is represented by a small number of plant functional types (PFTs), each with fixed properties irrespective of their predicted occurrence. This contrasts with natural vegetation, in which many plant traits vary systematically along geographic and environmental gradients. In the JSBACH DGVM, which is part of the MPI-ESM, we allowed three traits (specific leaf area (SLA), maximum carboxylation rate at 25 C (Vcmax25) and maximum electron transport rate at 25 C (Jmax25)) to vary within PFTs via trait–climate relationships based on a large trait database. TheR2 adjusted of these relationships were up to 0.83 and 0.71 for Vcmax 25 and Jmax 25, respectively. For SLA, more variance remained unexplained, with a maximum R2 adjustedof 0.40. Compared to the default simulation, allowing trait variation within PFTs resulted in gross primary productivity differences of up to 50 % in the tropics, in > 35 % different dominant vegetation cover, and a closer match with a natural vegetation map. The discrepancy between default trait values and natural trait variation, combined with the substantial changes in simulated vegetation properties, together emphasize that incorporating climate-driven trait variation, calibrated on observational data and based on ecological concepts, allows more variation in vegetation responses in DGVMs and as such is likely to enable more reliable projections in unknown climates.


Introduction
Terrestrial vegetation plays a pivotal role in land-atmosphere interactions, modifying carbon, water and heat fluxes via biochemical processes such as photosynthesis and respiration or via biophysical vegetation properties such as stomatal conductance and albedo.Therefore, a correct representation of terrestrial vegetation and its dynamics in Earth system models (ESMs) is essential, especially for future climate projections.In these models, vegetation can be simulated by dynamic global vegetation models (DGVMs), which integrate vegetation dynamics with land surface models.They are developed to predict vegetation distribution and fluxes of energy, water and carbon under past, current and future climate regimes.As such, they allow analysis of transient vegetation responses and feedbacks to climate (Foley et al., 1998;Prentice et al. 2007).
Compared to earlier models, vegetation dynamics and interactions between the biosphere and atmosphere have been Published by Copernicus Publications on behalf of the European Geosciences Union.

L. M. Verheijen et al.: Trait variation impacts on global vegetation model performance
much improved over the past decade, although many issues still need to be resolved (Quillet et al., 2010).One of these issues is the way in which plant functional types (PFTs) are used to represent vegetation.PFTs are classes of plant species with presumably similar roles in an ecosystem, responding in a comparable manner to environmental conditions like water and nutrient availability (Lavorel et al., 1997).They are defined by a combination of attributes such as plant growth form (e.g., trees, shrubs, grasses, herbs), phenology (evergreen, raingreen, summergreen) and bioclimatic tolerances (e.g., minimum temperature requirements).Notably, most current PFT classifications use constant parameter values for some key plant traits, i.e., plant properties that reflect the way plants cope with their environment (Violle et al., 2007).Using constant plant traits in PFTs has serious limitations (Ordoñez et al., 2009;Van Bodegom et al., 2012), as such constant values contrast with observed trait variation (Ackerly and Cornwell, 2007;Freschet et al., 2011;Westoby et al., 2002), therefore not accounting for local environmental constraints.Furthermore, global trait database analyses have shown that the variation in plant traits is large within PFTs and often even greater than the difference in means among PFTs (Laughlin et al., 2010;Wright et al., 2005a).Even though PFTs may capture a significant part of the global plant trait variation, a large part (for some traits up to 75 %) may still be unexplained (Kattge et al., 2011) and is thus not represented in current DGVMs.
Given that plants can adjust to the environment via changes in traits, and that such changes influence ecosystem functioning (Diaz et al., 2004;Lavorel and Garnier, 2002), implementing these trait-driven interactions within PFTs of DGVMs is highly relevant.By simulating variation in plant trait responses, it is better possible to quantify plant adaptation to climate and to account for plant-atmosphere feedbacks in DGVMs.Indeed, the limitations of inflexible PFTs is increasingly acknowledged by modelers and several attempts to allow for more (trait) variation within PFTs can be found in, for example, JeDi-DGVM (Pavlick et al., 2012) or aDGVM (Scheiter et al., 2013).However, none of the approaches so far tried to maximally include trait variation based on observational trait data and capture multiple sources of this variation by relating trait data to environmental variables.This inclusion may be possible by making use of the relationship between traits and multiple environmental drivers as occurring at both regional and global scales, for instance in relation to temperature, water and nutrient availability or disturbances (Ordoñez et al., 2009;Van Ommen Kloeke et al., 2012;Wright et al., 2005b).Such relationships between environmental conditions and traits can potentially be understood via ecological assembly theory, which describes the processes that determine species assemblages (Cornwell and Ackerly, 2009;Cornwell et al., 2006;Götzenberger et al., 2012).An important abiotic assembly process is habitat filtering (Keddy, 1992), which describes how local environmental drivers (e.g., soil fertility or pre-cipitation) constrain the range of potential species and related trait values in a given habitat.For many traits, such as specific leaf area, leaf nitrogen and wood density, this contributes to trait convergence within communities, prevailing over trait divergence (Freschet et al., 2011;Swenson and Enquist, 2007) and resulting in global relationships between community trait means and climatic drivers.Community trait means match closely to the model resolution of DGVMs, given that the majority of the models provide a mean value of a trait for a PFT (in a given grid cell).By identifying the environmental drivers of variation in community trait means, multiple causes of variation are determined and the observed natural trait variation can be modeled with a high level of accuracy.
The aim of this study is thus to model climate-driven trait variation within PFTs, determined as comprehensively as possible from observations, and to identify the impacts through vegetation responses and vegetation-atmosphere feedbacks on DGVM model behavior.So far, observational trait data have been used to derive mean parameter values and uncertainties for PFTs in the context of different vegetation models (Kattge et al., 2009;Ziehn et al., 2011), but this study provides a proof of concept to predict mean parameter values per PFT for each individual grid cell, based on a systematic analysis of observed trait-climate relationships using a large observational trait database.
The trait-climate relationships were implemented in the JSBACH DGVM, which is part of the Max Planck Institute Earth System Model (MPI-ESM) (Brovkin et al., 2009;Raddatz et al., 2007;Roeckner et al., 2003) and is used in the context of ESMs and carbon cycle model intercomparisons (Friedlingstein et al., 2006).We simulate variation in three originally constant and PFT-specific key leaf traits in JS-BACH, for which sufficient georeferenced observational data could be selected to calibrate the traits.These traits were specific leaf area (SLA), maximum carboxylation rate at a reference temperature of 25 • C (Vcmax 25 ) and maximum electron transport rate at 25 • C (Jmax 25 ).To determine relationships between these three traits and climatic drivers for each PFT, community means of trait values from the TRY global plant trait database (Kattge et al., 2011) were related to global climatic data.Based on these relationships between traits and climate, SLA, Vcmax 25 and Jmax 25 were reparameterized for each grid cell on a yearly basis, depending on the local climatic conditions in a grid cell.Compared to the default model, this enabled enhanced feedbacks between plants and environment, as traits within natural PFTs could vary dynamically in space and time.A simulation with variable traits is compared for trait distribution, productivity and vegetation distribution to a default simulation with the original trait values of the model, and to an additional simulation with constant, but observation-based, trait values.
As DGVMs are already parameterized to produce approximate realistic results, our variable traits simulation will not necessarily approach reality better than the default simulation.Therefore, in this paper the focus lies on simulation intercomparisons to evaluate the impact of incorporating climate-driven trait variation.The simulations are not meant to represent current vegetation or climate states, but were run into equilibrium states with their own internal climate and vegetation composition in order to study the impacts of trait variation within PFTs.However, we include comparisons with global biomass and gross primary productivity (GPP) estimates and potential vegetation distribution to evaluate the realism of the various simulation results.

Model description
Simulations were performed with the ESM of the Max Planck Institute (MPI-ESM).The model setup consisted of the JSBACH DGVM, a land surface model (Raddatz et al., 2007) with a vegetation dynamics module (Brovkin et al., 2009), coupled to the atmosphere model ECHAM5 (Roeckner et al., 2003) to allow for full inclusion of vegetationatmosphere feedbacks.The model version on spatial resolution of T63 (approx.1.875 • , dividing the world into 18 432 grid cells) with 47 atmospheric layers was used, with atmospheric CO 2 concentration kept constant at 353.9 ppm.Seasonal sea surface temperatures and sea ice were prescribed from a default simulation with the MPI ocean model (MPI-OM) (Marsland et al., 2003).Vegetation dynamics were interactive to allow for vegetation shifts and vegetationatmosphere feedbacks through trait dynamics.Terrestrial grid cells contained multiple PFTs, each occupying a certain fraction, depending on their competitive ability (based on net primary productivity -NPP).This study does not include crops and pastures, but focuses on responses of natural vegetation, represented by eight PFTs.These were tropical broadleaved evergreen trees (TrET), tropical broadleaved deciduous trees (TrDT), extratropical (both temperate and boreal) evergreen trees (ExTrET), extratropical deciduous trees (ExTrDT), raingreen shrubs (RgSh), cold/deciduous shrubs (DSh), C3 grasses (C3G) and C4 grasses (C4G).No anthropogenic land use or land cover change was simulated.

Selected trait and climate data
All PFT-dependent trait parameters in the model for which sufficient georeferenced observational trait data could be collected were selected.This resulted in selection of 3 traits; SLA (m 2 kg −1 carbon), Vcmax 25 (µmol m −2 s −1 ) and Jmax 25 (µmol m −2 s −1 ).SLA is used to determine the amount of carbon that can be stored in the green and reserve carbon pools.In contrast to most DGVMs, SLA is decoupled from the phenology routine in JSBACH and therefore does not play a role in determining productivity.Vcmax 25 (µmol m −2 s −1 ) and Jmax 25 (µmol m −2 s −1 ) are the reference values at 25 • C used for calculation of the Vcmax and Jmax values at the actual temperature.These traits are part of the photosynthesis routine in JSBACH, which is modeled following Farquhar et al. (1980) for C3 plants and Collatz et al. (1992) for C4 plants (see Supplement S1 for a more detailed description with equations of the functioning of traits in JSBACH).Vcmax 25 also determines the reference dark respiration at 25 • C (R d,25 ).Leaf carbon assimilation is the lowest of carboxylation and electron transport rate (calculated from actual Vcmax and Jmax) minus dark respiration, and scaled to canopy level with the use of the leaf area index (LAI).Supplement S2 provides a description and a flowchart with the role of the selected traits in JSBACH processes and the way variation in these traits might propagate in the model.
Only georeferenced trait observations from field sampling and field experiments were used and selected from the TRY database (Kattge et al., 2011), with additional data for SLA from the database by Van Bodegom et al. (2012) and for Vcmax 25 and Jmax 25 from Domingues et al. (2010) (see Table 1 for all references).From these databases, only Vcmax and Jmax data that could be standardized to a reference temperature of 25 • C were used.Most Vcmax and Jmax values had already been standardized to this temperature via the formulation of the photosynthesis model by Farquhar et al. (1980) used in the context of JSBACH (Kattge and Knorr, 2007).For other records for which the temperature during measurement was recorded, standardization was done according to this formulation.For C4 grasses, variation in PEPcase CO 2 specificity (PEP, mmol m −2 s −1 ) instead of Jmax 25 was determined, following Collatz et al. (1992).Since for C4 grasses insufficient observational PEP and Vcmax 25 data were available, these traits were estimated more indirectly applying equations from Simioni et al. (2004), who determined PEP and Vcmax 25 based on leaf nitrogen (N) content (g m −2 ) (regressions based on two C4 species).Therefore, for C4 grasses, additional information on leaf N was obtained from the TRY database.
Each trait observation was linked to a PFT based on information on growth form (shrub, grass, tree), leaf habit (deciduous/evergreen) and photosynthetic pathway (C3/C4) for the species involved.The climatic domain of a species (tropical, boreal etc.) was determined based on the Köppen-Geiger climate classification (Kottek et al., 2006) and applied to the georeferenced observations.

Trait-climate relationships
The observed natural variation of SLA, Vcmax 25 and Jmax 25 was described for each PFT by the simplest models possible, using linear regressions between individual traits and multiple observational climatic variables.The derived traitclimate relationships likely reflect habitat filtering, which commonly leads to trait convergence at community scales.The relationships were implemented in JSBACH to predict traits based on the internally simulated climate.As a consequence, in our choice of climate variables, we were constrained by the environmental predictors that are modeled in JSBACH.
Global climate data on mean annual precipitation (MAP, mm yr −1 ), mean annual relative humidity (Reh, %), mean annual temperature (MAT, • C), and mean temperature of coldest and warmest month (T min and T max , • C) were collected from a global 10 min gridded dataset of mean monthly climate data based on weather stations from the Climatic Research Unit (CRU) (New et al., 2002).Mean annual net shortwave radiation (NSWR, W m −2 ) was calculated based on distance to sun and percentage sunshine (from CRU) according to Allen et al. (1998) on a global 30 min spatial resolution.Mean annual soil moisture was taken from GLEAM, a methodology that estimates soil moisture and evapotranspiration based on remotely sensed data at a global 15 min spatial resolution (SoilMoist, m 3 m −3 for 1 m depth; Miralles et al., 2011).Given the available modeled climatic drivers in JSBACH, soil moisture was the best available estimate of a drought index.Soil moisture was the only edaphic control of trait variation, since soil N was not modeled in this version of JSBACH.
To determine trait-climate relationships, we tested all possible climate combinations in relation to each trait.Regressions with the highest explained variance (highest R 2 adjusted ) were selected, after checking for significance of climatic drivers and distribution of residuals.Table 2 shows for each PFT the selected climatic drivers of the best model, as well as their directional relationship with the traits (third column), and the corresponding R 2 adjusted .In case of co-linearity between environmental drivers (Pearson's correlation higher than 0.7 or lower than −0.7), the driver with lowest significance was omitted.Due to the low number of entries for Vcmax 25 and Jmax 25 , data for the two tropical tree PFTs were combined, as well as data for the two shrub PFTs.Even though we acknowledge that mean trait values for evergreen and deciduous vegetation types will differ in reality, combining these will give better estimates than modeling these traits on their few individual data points alone.Due to differences in other vegetation characteristics, these PFTs will still differ in their emergent properties, as the results will show.In total there were eight relationships for SLA and six for both Vcmax 25 and Jmax 25 .These relationships were subsequently implemented in JSBACH.
In addition, to allow trait calculations to go beyond climatic regions used in the regressions, but to still maintain traits within ranges as observed in nature, we applied two constraints on simulated trait ranges.First of all, predicted traits were constrained to the 2.5-97.5 % quantiles of all individual observations within a PFT (instead of the community means).Secondly, for Jmax 25 and Vcmax 25, an additional constraint was applied to maintain the strong physiological correlation between them (Medlyn et al., 2002;Wullschleger, 1993), keeping Jmax 25 and Vcmax 25 values within the 95 % confidence interval of the linear regression among these two traits (as determined per PFT).Predicted values outside this range were adjusted to the confidence interval border, based on the shortest distance needed to reach the border.Since this version of JSBACH did not simulate N-cycling, soil N could not be used to parameterize Vcmax 25 , even though Navailability is a strong determinant (see references in Kattge et al., 2009).However, observational data was used to limit the above two constraints on Vcmax 25 , in this way indirectly incorporating N-limitations via observed Vcmax 25 values.

Simulation setups
Three different scenarios were performed: (1) a simulation containing the default parameterization of JSBACH with constant parameter values per PFT, based on Raddatz et al. (2007) and Kattge et al. (2009), and adapted to approximate realistic vegetation functioning within the vegetation dynamics module (Brovkin et al., 2009); hereafter called "de-fault simulation", (2) an "observed traits simulation" again with constant values for SLA, Vcmax 25 and Jmax 25 per PFT, but based on observational data only (the weighted community means for each PFT, from here on called "observed global means"), and (3) a "variable traits simulation" in which traits were allowed to vary depending on local climatic conditions.In the latter simulation, at the beginning of every year, each of the three key traits was reparameterized for each PFT in every terrestrial grid cell of the world, depending on the local simulated climatic conditions in each grid cell.To get vegetation into quasi-equilibrium with the simulated climate (but with CO 2 concentration fixed), the coupled model (JSBACH/ECHAM5) was run for 150 years with vegetation dynamics in an accelerated mode (i.e., vegetation was simulated 3 times for each year of simulated climate).Next, to get the slow soil carbon pools into equilibrium, the uncoupled carbon model of JSBACH ("CBALANCE") was run for 1500 yr.The coupled model was continued for 160 yr with updated carbon pools (with the vegetation module in an accelerated mode for the first 50 yr), until carbon pools and vegetation cover had reached equilibrium.The last 10 yr were averaged and used for further analysis, which in this setup (with prescribed seasonal sea surface temperatures) was sufficient to account for the interannual climate variation.

Data comparison
Model performance was investigated by comparing vegetation distribution, GPP and biomass of the simulations with observations.Vegetation distribution of the dominant PFT (Fig. 5) was compared to the potential (natural) vegetation map of Ramankutty and Foley (1999), which in part is used to initialize simulations in JSBACH.The mismatch in spatial resolution (the potential vegetation maps is at a global 30 min spatial resolution) was solved by counting the number of PFTs according to Ramankutty and Foley (1999) that were present in each JSBACH grid cell.The PFT with the highest occurrence was compared with the simulated dominant PFT of JSBACH.To match the PFT classifications, aggregated PFTs were constructed (Supplement Table S6.1 and Fig. S6.1), as the PFTs of Ramankutty and Foley (1999) could not be directly related to those of JSBACH.Temperate and boreal evergreen trees were merged to match extratropical evergreen trees in JSBACH.The same was done for deciduous trees.Shrubs (marginal in JSBACH) were merged into a single PFT ("shrubs") in both JSBACH and the vegetation map.As savannas consist of grasslands and woodlands, matches between savannas and C4 grasses or tropical broadleaved deciduous trees were both classified as correct.Tundra and mixed forests were omitted from the comparison as no equivalent PFT was available in JSBACH.Model performance was determined with Cohen's kappa (κ) (Cohen, 1960), which is the proportion of agreement between vegetation maps, while accounting for chance agreement.Grid cells with sea or ice as the dominant cover were not taken into account, and neither were cases in which more than 1 PFT shared the highest occurrence (i.e., equal number of grid cells).This resulted in a comparison of 2819 grid cells.
Latitudinal patterns of median GPP were compared with data taken from Beer et al. (2010), who combined observa-tional data (eddy covariance fluxes) with diagnostic models to approximate GPP.
Simulated biomass was compared to estimates from Robinson (2007), which include adjusted estimates of belowground (root) carbon in plant biomass.Estimates based on root/shoot quotients from Mokany et al. (2006) and Grace et al. (2006) were used.As we simulated only natural vegetation and the estimates of Robinson (2007) are current estimates including land use change and crops, comparing global total biomass would result in an overestimation by our simulations.Therefore, comparisons were made per m 2 per (aggregated) PFT.JSBACH does not have a separate savanna-like PFT.Therefore, tropical forests and C4 grasslands were averaged and compared with averages of tropical forests and savanna.Furthermore, extratropical trees were averaged and compared with averaged temperate and boreal trees of Robinson (2007).Mediterranean shrubs and tundra were omitted as there were no comparable PFTs in JSBACH.For the tropics, an additional comparison of biomass was made to Saatchi et al. (2011).

Mismatch between observed trait data and default trait settings
Fig. 1 shows the observed trait ranges as collected from the databases., respectively).Such mismatch between observational trait data and default trait values has also been reported by Kattge et al. (2011), who showed how PFT-specific constant SLA values in DGVMs often differ from the means of observed data, and in some cases are at the low end of the observed data distribution.In JSBACH, SLA default values were always lower than the observed global means, except for the tropical broadleaved deciduous trees, but they almost all fell within the 25 % quartiles of the observed range of SLA, except for deciduous shrubs.However, for both Vmax 25 and Jmax 25 , default values fell outside the 25 % quartiles for more than half of the PFTs (five out of eight PFTs for both), and in some cases even outside the minimum and maximum values (tropical deciduous trees, C4 grasses and C3 grasses for Vcmax 25 and Jmax 25 , and additionally extratropical evergreen trees for Jmax 25 ).In contrast to SLA, there was no clear direction of the trait differences for either Vcmax 25 or Jmax 25 : half of the default values were lower than the observed global means (tropical deciduous trees, extratropical evergreen trees, raingreen shrubs and C4 grasses), the other half higher.These differences point to a strong mismatch between default trait parameters in JSBACH with observed natural trait means.
Figure 1 also shows that within each PFT large trait variation was apparent, as well as large overlap in trait ranges among PFTs, which corresponds to other trait analysis (Laughlin et al., 2010;Wright et al., 2005a).This was most manifested for SLA, where all PFTs overlapped with each other.For Vcmax 25 , the only PFTs of which the 95 % confidence interval did not overlap were both tropical evergreen trees and C4 grasses with C3 grasses, and for Jmax 25 both tropical trees and deciduous shrubs with C3 grasses.

Simulated trait variation based on climatic drivers
The observed variation in PFT-specific community means (as given by Table 1) was related to environmental drivers by applying multiple regression.Table 2 (third column) presents for each trait and each PFT the selected climatic drivers used in the regressions and their directional relationship with the traits.All regressions were significant, except for C3 grasses for Vcmax 25 and for Jmax 25 .Although the selected environmental drivers differed among PFTs, NSWR and MAP were most frequently selected as drivers for SLA, as were MAP, MAT and RH for Vcmax 25 , and MAP for Jmax 25 .R 2 adjusted was up to 0.83 for Vcmax 25 and 0.71 for Jmax 25 .For SLA, more variation remained unexplained, with a maximum value of 0.40 for R 2 adjusted .Fig. 2 presents the predicted trait ranges in JSBACH, as calculated by applying the regressions on simulated climate in JSBACH.As in Fig. 1, grey diamonds and black circles again represent default trait values and observed global means.For each grid cell, only simulated trait values of the dominant PFTs are selected to prevent trait values of marginal PFTs with almost no ecosystem impact to obscure trait ranges most important for model output.Raingreen shrubs were almost never dominant, and predicted variation was low, while deciduous shrubs were never dominant.
For all three traits, a strong mismatch of default traits with the simulated trait variation was apparent: for SLA default values (grey diamonds) fell (at least) outside the 25-75 % interquartile range in four PFTs (vs.one in the observed data), and for Vcmax 25 and Jmax 25 for every PFT.
In contrast to the default and observed trait simulation, where PFTs had one constant value per trait, a large climatedriven trait range was predicted for the variable traits simulation.The higher SLA of deciduous versus evergreen trees observed in natural vegetation was reflected in the simulated SLA ranges.As for the observational trait data (Fig. 1), there was large overlap in predicted trait ranges among PFTs as well (Fig. 2), particularly for Vcmax 25 and Jmax 25 .Simulated SLA ranges were narrower and overlapped less among PFTs compared to observed ranges.For some PFTs (e.g., ex- tratropical deciduous trees and C3 grasses for Vcmax 25 and Jmax 25 ), simulated trait variation for Vcmax 25 and Jmax 25 (Fig. 2) was larger than observed variation (Fig. 1).Predicted traits may have fallen outside the ranges of observed variation when climatic ranges as simulated in JSBACH exceeded observational climate ranges used to derive the trait-climate relationships.However, this did not result in unrealistic trait www.biogeosciences.net/10/5497/2013/Biogeosciences, 10, 5497-5515, 2013 values, as limits were set on maximum trait values based on individual observations.Furthermore, we implemented a trade-off between Vcmax 25 and Jmax 25 , keeping them within realistic ranges (even outside climatic ranges used for determination of the regressions).For SLA, such trade-off was not taken into account, but unrealistic SLA values were not very likely as SLA was based on a much larger dataset, covering a very large range of climatic conditions.

Inclusion of trait variation alters predicted global patterns of productivity
The observed traits and variable traits simulation showed a similar pattern of median GPP along the latitudinal gradient compared to the default simulation (Fig. 3), but in absolute values GPP was much higher in the variable traits simulation in the geographic zone between about 30 • S and 30 • N (∼ 50 % around the equator and twice as much between 25 and 30 • N).GPP of the observed traits simulation also reached higher values (up to twice as much) as the default between 25 and 30 • N, but had lower values than the variable traits simulation between 20 • S and 15 • N, and highest values (64-67 %) at higher latitudes.NPP differences between 30 • S and 30 • N were less profound than for GPP, showing that trait variation affects global patterns of GPP and NPP in different ways.Again, at higher latitudes, the observed traits simulation predicted highest NPP.For the variable traits simulation, the reduction from gross to net primary productivity around the tropics is much stronger than for either other simulation, resulting in smaller NPP differences among simulations compared to GPP differences (even though differences may still go up to 75 % between 25 and 30 • N for both observed and variable traits simulations).The large reduction in differences between NPP and GPP among simulations cannot solely be attributed to a higher respiration of PFTs due to higher Vcmax 25 (see Fig. 4), and is mainly a consequence of net productivity exceeding the maximum sizes of the different carbon pools (e.g., leaves, wood, reserves) as defined in JSBACH, resulting in respiration of excess carbon.This loss of excess carbon is visible in the tropical regions for all simulations, but for the variable traits simulation this resulted in a proportionally larger amount of productivity removed.This was enhanced by the higher SLA in these regions, which partly determines the amount of carbon stored in the living parts of the plants, consisting of leaves, fine roots and sapwood.
These results show that, while reproducing global patterns of productivity, incorporating global trait variation in the model leads to strong changes in predicted productivity.Even though GPP is not only determined by the photosynthetic parameters (Vcmax 25 and Jmax 25 ), as it is affected by, for example, water availability as well, on average, shifts in GPP were generally accompanied by shifts in similar directions by mean Vcmax 25 and Jmax 25 (weighted by fractional coverage of PFTs), compared to the default simulation (Fig. 4).Exceptions occurred, as for example, the drop in Vcmax 25 and Jmax 25 above 40 • N did not lead to a coinciding drop in GPP.As expected, as SLA does not play a role in productivity, changes in GPP could not be related to similar shifts in SLA.

Major shifts in vegetation distribution
Figure 5 and Table 3 show how the global distribution of dominant vegetation types as predicted by the simulations strongly changes when incorporating observed global trait values or including climate-driven trait variation.A PFT was considered dominant if it had the highest fractional coverage in a grid cell; this ranged from coverage of almost 100 % (mostly in tropical regions) to only 30 % in some areas at higher latitudes (see Supplement S4 for fractional coverage of the dominant PFTs).Predicted dominant PFTs differed from the default simulation in 35.4 % of the terrestrial grid cells for the variable traits simulation and in 50.5 % of the grid cells for the observed traits simulation.
In the observed traits simulation, tropical evergreen trees, dominant in the default simulation, were taken over by deciduous trees in Africa, South America and Australia and extratropical deciduous trees were replaced by extratropi- cal evergreen trees as the dominant PFT, resulting in less spatially heterogeneous dominant vegetation.In the variable traits simulation, these shifts occurred as well (see Fig. 5 and Table 3).However, in contrast to the observed traits simulation, both changes in dominant tree cover only occurred in limited areas, which resulted in more spatial variation in vegetation in the areas where trees were dominant.The shifts from tropical evergreen to tropical deciduous trees cannot be explained by Vcmax 25 and Jmax 25 , since these tropical PFTs were parameterized with the same values for these traits.The most profound difference between these tropical PFTs seems to be their leaf turnover rate, which is higher for the deciduous than for the evergreen trees.As a consequence, tropical deciduous trees had somewhat lower leaf area index, which meant lower productivity in favorable periods, but also less carbon loss in more stressful circumstances (e.g., drier periods).In some areas, this could have resulted in a higher total yearly NPP for tropical deciduous trees, thereby outcompeting evergreen trees.Another shift in predicted dominant vegetation in the variable traits simulation was an increase in C4 grasses (from 4.8 % in the default simulation to 8.4 %).This occurred mostly in Africa and Australia at the expense of tropical trees and raingreen shrubs.This expansion of C4 grasses below the Sahara coincided with higher fractions of burned area, which promoted the expansion of grasses at the cost of trees.
In the variable traits simulation, bare ground increased as the dominant cover type (from 15.4 % in the default simulation to 19.5 %) in the southwest of the United States (and Mexico), northern Canada and northeast of Siberia at the expense of C3 grasses and deciduous trees, resulting in a shift of the boreal treeline toward lower latitudes.These shifts often coincided with a decrease in Vcmax 25 (see Supplement S5), suggesting lower productivity and consequently less expansion of these PFTs.

Modulation of climate by traits
Predicted climate differed among simulations.In the variable traits simulation enhanced climate-vegetation feedbacks were possible both through changes in traits (modifying PFT properties), as well as through shifts in vegetation distribution (see Supplement S2 describing such interactions).However, differences in climate cannot be related unambiguously to either traits or vegetation shifts.As differences in temperature in unvegetated areas like Greenland and Antarctica show, there does not have to be a direct spatial relation between climate and vegetation.Moreover, some vegetation shifts (compared to the default) do occur in places where climates differ among simulations (e.g., northern Canada), whereas in other areas no such effect occurred (e.g., in tropical zones).
Precipitation differences were local and showed no clear spatial pattern (Fig. 6c and d), but in the variable traits and observed traits simulations it was drier in Canada, Asia and Australia, as well as in large parts of the Amazon rainforests compared to the default simulation.
Compared to the default simulation, in the observed traits simulation mean annual surface-air temperatures were profoundly higher (over 1 • C) in eastern Siberia, Alaska, the US and Australia, and lower (up to 1 • C) in large parts of Europe and Russia, South Africa and South America (Fig. 6a), meaning that (at least) temperature is very sensitive to parameterization of traits.Temperature differences were less profound between default and the variable traits simulation (Fig. 6b), but still went up to around 1 • C. Changes in temperature did not correlate with clear changes in traits or vegetation shifts (e.g., tree-grass shifts), but in the Southern Hemisphere (Australia, Africa, South America) corresponded to differences in transpiration, where cooler areas coincided with higher transpiration.This could be related to the higher Vcmax 25 in these areas, resulting in higher GPP and consequently an increase in transpiration.
The differences between the observed traits simulation and variable traits simulation indicate that by allowing traits to vary and respond to environmental conditions (as in the variable traits simulation), feedbacks between climate and traits result in more moderate temperature shifts, showing the significant effect of adaptive traits on climate.

Comparison of model output with observational data
Cohen's κ, indicating the correspondence of the global map of potential (natural) vegetation of Ramankutty and Foley (1999) with simulated vegetation distribution, was 0.289, 0.282 and 0.334 for the default, observed traits and variable traits simulations, respectively.These values were somewhat lower than the κ's of other DGVMs, e.g., κ = 0.40 for the Lund-Potsdam-Jena DGVM (LPJ) or κ = 0.42 for LPJ with implementation of plant hydraulic architecture, both with 18 PFTs (Hickler et al., 2006), or κ = 0.42 for a consensus map of multiple DGVMs with 6 PFTs (Cramer et al, 2001).However, these models were compared to different vegetation maps and had a different number of classes, making comparisons difficult.Our simulations with 7 vegetation classes performed less well, but this might partly depend on the chosen vegetation map (Ramankutty and Foley, 1999).
A low κ means that for either simulation, a large proportion of the grid cells did not match the potential vegetation map (Supplement Fig. S6.2).However, there is a substantial increase in similarity to observed vegetation from the default and observed traits simulation toward the variable traits simulation.
Mismatches occurred in large parts of the US and Canada, where the simulations predicted mostly C3 grasslands, while according to the potential vegetation map also forests should be present.The same holds for large parts of Europe.The potential vegetation map shows less bare ground than any simulation, resulting in mismatches in the US, but also other parts of the world.Furthermore, almost the whole continent of Australia did not correspond to this map; shrubs and savanna are dominant according to the map, and even though the models did predict C4 grasses there, it was in different areas.Where the default and observed traits simulation had low correspondence with the vegetation map in Africa and South America, the variable traits simulation performed better, mainly with respect to the tropical trees.Even though differences in performance are small, the variable traits simulation matched the potential vegetation map most closely.
Comparing latitudinal patterns of median GPP with estimates from Beer et al. (2010) (thin blue line in Fig. 3), each simulation produces substantial higher GPP at most latitudes (on average leading to 2, 2.4 and 2.6 times higher GPP for the default, observed traits and variable traits simulations, respectively), with the default simulation in general having the smallest differences.variable traits simulation, GPP particularly in the tropics was overestimated.As nothing has been changed in the photosynthesis schemes, this probably points to a parameter mismatch in the (default) model, as in the variable traits simulation the traits are supposed to be more close to realistic values than in the default model.Chen et al. ( 2012) suggest that global GPP estimates based on remotely sensed LAI are underestimated by 9 % when leaf clumping is not taken into account, as this would result in an underestimation of the contribution of shaded leaves to GPP, with the strongest underestimation occurring in the tropics.This implies that the estimates of GPP by Beer et al. (2010) might be too low, as remotely sensed LAI (and fAPAR) data were used to extrapolate GPP estimates from flux-towers to global maps, and thus the actual differences with our simulations may be less.Also many other DGVMs show higher GPP in the tropical areas than the observed median GPP by Beer et al. (2010), and the GPP estimates of the variable traits simulation still do fall within the upper range of GPP predictions by other DGVMs.
Comparisons of simulated biomass per m 2 per (aggregated) PFT with current biomass estimates by Robinson (2007) show that for the combined tropical trees and savannas the variable traits simulation (13.12 kg C m −2 ) most closely matched biomass estimates (13.14 kg C m −2 ) (Ta-ble 4).For extratropical (temperate and boreal) forest and temperate grasslands, the default simulation underestimates carbon in vegetation (7.08 kg C m −2 ).Both the observed traits and variable traits simulation are closer (9.64 and 9.65 kg C m −2 ) to the global estimates for forests (8.83 kg C m −2 ), although they overestimate biomass.For temperate grasslands, either simulation underestimates biomass, with the default simulation deviating the least from global estimates (0.16 vs. 0.67 kg C m −2 ), even though differences among simulations are small.Overall, this implies that of the three simulations, the variable traits simulation is closest to global biomass estimates per PFT per m 2 .
The largest differences in GPP were found in the tropical areas.Comparing simulated biomass of the tropical zone with estimates of Saatchi et al. (2011), either simulation overestimates biomass in most areas (see Supplement S7), especially in the south of Brazil, and large parts of southern Africa (areas with dry forests and savannas).In general, for either simulation the areas with high biomass (> 200 Mg C ha −1 ) are much more extensive than according to Saatchi et al. (2011).However, it has to be noted that for these areas with high biomass, the uncertainties in biomass estimates of Saatchi et al. (2011) are large, ranging from 25 % to over 50 %, meaning differences might be smaller in reality.

Discussion
The aim of this study is to identify the impacts of climatedriven trait variation within PFTs through vegetation responses and vegetation-atmosphere feedbacks on DGVM model behavior.We determined observed trait variation within PFTs as comprehensively as possible based on relationships between measured traits and climate and soil moisture, representing major assembly processes by the abiotic environment.As model intercomparisons have shown that large uncertainties exist in projections of land carbon uptake by DGVMs (Cramer et al., 2001;Friedlingstein et al., 2006, Sitch et al., 2008), incorporation of variation in vegetation responses is important to allow more feedbacks between vegetation and climate and to increase plausibility of model predictions, especially under strong climate change.
Here, we incorporated variation in plant responses based on relationships between observed trait and climate data into the JSBACH DGVM, which revealed profound effects on carbon fluxes and vegetation distribution.

Challenges of modeling trait variation based on trait-climate relationships
In our approach, we used trait-climate relationships to describe the observed natural trait variation and implemented these in JSBACH.These relationships identify and capture multiple (approximate) drivers of natural trait variation, and are thought to reflect abiotic assembly processes.They integrate multiple vegetation responses at different temporal and spatial scales, including acclimation, adaptation of species and species replacement into a spatially and temporally varying trait mean.
One of the issues that might affect the reliability of model performance is the fact that we derived climate-driven trait variation for each trait independently, not explicitly accounting for trait trade-offs.The advantage of a changed trait often implies a cost in another trait (e.g., the leaf economics spec-trum of Wright et al., 2004), possibly dampening the positive effect of the other trait and modifying plant performance.Implementing trait trade-offs helps to restrict possible trait combinations or plant responses, especially when simulated traits are modeled outside the observational climatic ranges used to derive the trait-climate relationships.However, tradeoffs can only be included if the traits involved are explicitly represented in the model and if they are PFT specific.This complicates including and evaluating the importance of trade-offs in DGVMs.In our situation, we took care of the trade-off among Vcmax 25 and Jmax 25 by constraining the values of one trait by the other.Moreover, in JSBACH, changes in these traits are coupled to respiration and transpiration processes as well, likely leading to consistent responses.Trade-offs between SLA and other traits were not accounted for, given the specific way SLA in JSBACH is used in relation to productivity, independent of LAI and phenology.This decoupling of SLA expresses that early phenology is strongly driven by remobilization of carbon from previous year reserves and not from current year productivity.Instead, in JSBACH phenology is determined by the environmental drivers temperature and soil moisture.Therefore, the commonly observed strong trade-off between SLA and leaf life span (LLS) (Wright et al., 2004) is not expressed in JSBACH, in contrast to many other DGVMs.In fact, by incorporating environment-driven variation in SLA, the decoupling between SLA and LLS has been partly diminished compared to the default model, more closely representing natural plant strategies.Trait trade-offs thus need to be evaluated depending on the model formulation.
We updated trait values once per year.This time step is a balance between computational efficiency and ecological realism.By this approach we avoid evaluating ontogenic impacts on trait values and whether environmental impacts differ for different parts of the growing season, for which currently insufficient information is available.Despite presumed differences in plasticity among species and PFTs (although L. M. Verheijen et al.: Trait variation impacts on global vegetation model performance we are not aware of studies analyzing this), current analyses on within-species variation of leaf economics traits (Kattge et al. 2011;Messier et al. 2010) suggest that all PFTs will be sufficiently plastic to adjust their leaf economic traits at yearly timescales to the extents forced by year-to-year differences in climatic drivers.The selected leaf traits in this study are rather plastic and can vary within a year (Dubey et al., 2011;Misson et al., 2006).Other traits, like wood traits, will have longer adaptation times than a year, as changes would involve individual or species turnover.SLA may be somewhat less plastic than the two photosynthetic traits, but a yearly reparameterization is reasonable, as leaves of deciduous species are shed every year and vary in relation to environmental changes (Ma et al. 2011).The maximum lifespan for leaves of evergreen PFTs is three years in JSBACH, meaning yearly SLA shifts occur in (on average) a third of the leaves, resulting in a slight overestimation of SLA variation in evergreens.However, the yearly shift in leaf trait values may not only reflect acclimation, but genetic adaptation and species replacements may also contribute.In contrast to a mechanistic approach, in our approach the impacts and (unknown) timescales of those processes leading to trait shifts do not have to be differentiated.

Implementing variation in PFTs: comparing approaches
Determining trait-climate relationships is broadly accepted and applied in ecology (Ordoñez et al., 2010;Wright et al., 2005b), but using these relationships in combination with observational trait data to calibrate trait variation and implement these relationships in DGVMs is novel.Recently, there are various attempts to implement more (trait) variation in PFTs or calibrate parameters on observational data.For example, Alton (2011) allowed for more variation in a number of traits in the land surface model JULES-SF via tuning of traits to ranges set by observational data (e.g., eddycovariance fluxes, satellite data), but this trait variation was modeled as stochastic processes.Some trait variation based on mechanistic principles has been incorporated in O-CN (Zaehle and Friend, 2010); there trait variation was not the primary goal, but a means to enable feedbacks in the modeling of nutrient cycles.Some DGVMs also implement the concept of environmental filtering, like the JeDi-DGVM (Pavlick et al., 2012).This DGVM models functional diversity by creating a continuum of PFTs or plant growth strategies (PGSs), which are determined by trait trade-offs and habitat filtering.However, this model does not include measurable traits (Reu et al., 2011) and as such does not allow for linking traits to observational trait data or evaluation of trait combinations.Also aDGVM (Scheiter et al., 2013) generates trait variation, and viable trait combinations are selected and inherited via species performance.This means that environmental filtering only acts on trait values through the next gen-eration.aDGVM has not been validated with observational data, nor does it include trait trade-offs.
In contrast to these approaches, we aim to comprehensively model trait variation in PFTs by identifying multiple drivers of trait variation and calibrating these relationships on observational trait data.Our proposed method is correlational and does not explain mechanistically the adaptation, acclimation or turnover processes behind trait variation (Pavlick et al. 2012).Nor does it fully account for constraints by biotic interactions, trait trade-offs or dispersal limitation on trait values (similar to most DGVMs).Still, in our opinion, it is an important and necessary step as it reflects the observed correlation between traits and climatic drivers (Niinemets, 2001;StPaul et al., 2012;Wright et al., 2005b) as changing with changing climate conditions.Importantly, it has the advantage that it does identify and quantify multiple abiotic drivers of trait means and in this way captures a large part of observed trait variation, as shown by the substantial R 2 adjusted of most regressions.

Implications of incorporating observation-based trait variation
The observed global mean trait values of natural vegetation as used in the observed traits simulation strongly deviated from trait values in the default simulation, indicating a mismatch between default PFT trait means and means obtained from natural vegetation.Moreover, either set of constant values contrasts strongly with the large trait variation observed in natural vegetation (Fig. 1).While we applied the most comprehensive database available today, we are aware that estimates of observed trait variation are still uncertain (Table 2) and need to be improved in future applications.Nevertheless, the wide range of observed trait values illustrates how simulations with constant traits do not reflect natural trait variation.In contrast, this variation was reflected by the variable traits simulation where trait variation represented abiotic assembly processes (Fig. 2).
To investigate the effects of trait variation on vegetationclimate feedbacks, it was essential to incorporate vegetation dynamics, to allow trait shifts to alter vegetation distribution and in this way modulate productivity and climate.In contrast to the simulations with constant traits (both default and observed traits simulations), the variable traits simulation enabled enhanced interactions between vegetation and climate to occur, via shifting traits.Such changes in traits did not have a direct effect on climate, but propagated indirectly to climate by modifying different plant properties and fluxes (see Supplement S2 for these pathways).These trait-climate interactions ultimately resulted in more spatial variation in dominant vegetation compared to the other two simulations.As such, predicted vegetation distribution is more a result of temporal dynamics in vegetation properties than is the case in the other simulations where these vegetation properties were prescribed.
In the current model setup, direct and indirect effects (via climate, changes in vegetation properties or distribution) of traits are not easily disentangled, which makes it difficult to pinpoint how and to what degree traits directly modulated model output.Our main aim however, is to identify the integrated impact of trait variation on model behavior via vegetation responses and feedbacks with climate.So, we were more interested in whether and to what extent trait variation alters model performance, than to know the exact pathways, as we acknowledged in advance this is difficult with a coupled model setup.Offline simulations with prescribed climate could help to reduce the number of possible variables causing the changes and may clarify some patterns.However, for our model, such offline simulations were for technical reasons not yet feasible.
Provided the strong effect exerted by climate on traits, major differences among the simulations in predicted vegetation distribution and productivity were expected, the latter especially when parameters that affect assimilation rate are concerned, as sensitivity analyses of DGVMs have shown (White et al., 2000;Zaehle et al., 2005).Indeed, for the observed traits simulation and variable traits simulation, this resulted in large differences in the new equilibrium state, both compared to the default and to each other.Not only were vegetation properties affected, but also climate changed and mean temperatures were altered by up to more than 1 • C.This confirms that JSBACH (and other DGVMs), are quite sensitive to changes in a few key parameters and this may cast some doubts on the reliability and predictive power of DGVMs in general.
As the simulations in this study provide equilibrium states, and as such do not necessarily correspond to current climate or vegetation composition, model comparisons with observational data must be interpreted with care, although it does provide insights in the realism of the simulations.DGVMs are parameterized to produce approximately realistic results, and therefore our simulations were not expected to approach observations better than the default simulation.Even though the variable traits simulations produced high GPP for tropical areas, its biomass estimates and vegetation distribution more closely resembled observational data than the default simulation.
The large differences in model output and the mismatch between default trait values in the model and observed trait variation in nature demonstrate that -besides a correct representation of plant physiology (e.g., photosynthesis, transpiration) -integration of ecological theory will be a step forward to help improve vegetation representation in DGVMs and ESMs.Allowing for variable vegetation responses to climate and soil will have important consequences for predictions by vegetation models, for vegetation distribution and productivity as well as for global current and future climate.A model intercomparison with this approach under elevated CO 2 projections, for which large uncertainties in predictions by current DGVMs exist (Friedlingstein et al., 2006;Sitch et al., 2008), should therefore be one of the next steps.

Conclusions
In this study we identified the impacts of modeling climatedriven trait variation in PFTs, calibrated on observational data, on JSBACH model behavior.The current mismatch of constant trait values in JSBACH with observed natural trait variation and the impact of incorporation of trait variation on model behavior with respect to vegetation distribution, productivity and global climate together emphasize the need for implementation of more observation-based trait variation and concomitant ecological concepts.The suggested approach, based on such data and concepts, reflects vegetation acclimation and adaptation to the environment, and will help enable more reliable modeling of vegetation behavior under unknown climates.

Fig. 2 .
Fig. 2. Simulated trait ranges of the dominant PFTs of each grid cell in the variable traits simulation: (a) SLA (m 2 kg −1 carbon), (b) Vcmax 25 (µmol m −2 s −1 ), (c) Jmax 25 (µmol m −2 s −1 ).For C4 grasses, PEPcase CO 2 specificity (PEP, in mmol m −2 s −1 ) instead of Jmax 25 is modeled.Box plots and symbols as in Fig. 1.PFTs as in Fig. 1 but without cold/deciduous shrubs (DSh), because this PFT was never dominant.Observed global means are added again for illustrative purposes.Note that direct comparison with simulated data is not appropriate here as the range of climatic conditions do not overlap for observed and simulated trait data.

Fig. 6 .
Fig. 6.Difference in mean annual temperature ( • C) and annual precipitation (mm yr −1 ).Observed traits simulation minus default simulation for (a) temperature and (c) precipitation, and variable traits simulation minus default simulation for (b) temperature and (d) precipitation.

Table 2 .
Properties of the selected trait-climate relationships used to parameterize traits.

Verheijen et al.: Trait variation impacts on global vegetation model performance the
effects of replacing default trait constants by observationbased constants from the effects of adding trait variation in the model.As this simulation mainly serves as a control, we will focus our results and discussion on the variable traits simulation.

Table 3 .
Dominant PFT coverage (%) of vegetated grid cells for the three different simulations.

Table 4 .
Biomass estimates fromRobinson (2007)compared to simulated biomass of the three simulations.