# A new approach to optimal discretization of plant functional types in a process-based ecosystem model with forest management: a case study for temperate conifers

Marc Peaucelle, Valentin Bellassen, Philippe Ciais, Josep Peñuelas, Nicolas Viovy

► To cite this version:

Marc Peaucelle, Valentin Bellassen, Philippe Ciais, Josep Peñuelas, Nicolas Viovy. A new approach to optimal discretization of plant functional types in a process-based ecosystem model with forest management: a case study for temperate conifers. Global Ecology and Biogeography, 2017, 26 (4), pp.486 - 499. 10.1111/geb.12557. hal-01584210

HAL Id: hal-01584210

https://hal.science/hal-01584210v1

Submitted on 15 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
This is the **accepted version** of the article:

Peaucelle, Marc; Bellassen, Valentin; Ciais, Philippe; [et al.]. «A new approach to optimal discretization of plant functional types in a process-based ecosystem model with forest management: a case study for temperate conifers». Global Ecology and Biogeography, Vol. 26, Issue 4 (April 2017), p. 486-499. DOI 10.1111/geb.12557

This version is available at https://ddd.uab.cat/record/218310

under the terms of the © IN COPYRIGHT license
1 A new approach for optimal discretization of plant functional types in a
2 process-based ecosystem model with forest management: a case study of
3 temperate conifers.

4 Marc Peaucelle¹, Valentin Bellassen², Philippe Ciais¹, Josep Peñuelas³,⁴, Nicolas Viovy¹

5 ¹Laboratoire des Sciences du Climat et de l'Environnement, CEA CNRS UVSQ, Gif-sur-Yvette, France

6 ²Centre d'Economie et de Sociologie appliquées à l'Agriculture et aux Espaces Ruraux, INRA, Dijon, France

7 ³CREAF, Cerdanyola del Vallès, Barcelona 08193, Catalonia, Spain

8 ⁴CSIC, Global Ecology Unit CREAF-CSIC-UAB, Cerdanyola del Vallès, Barcelona 11 08193, Catalonia, Spain

9

10 Marc Peaucelle: [marc.peaucelle@lsce.ipsl.fr](mailto:marc.peaucelle@lsce.ipsl.fr) (corresponding author)

11 Valentin Bellassen: [valentin.bellassen@dijon.inra.fr](mailto:valentin.bellassen@dijon.inra.fr)

12 Philippe Ciais: [philippe.ciais@lsce.ipsl.fr](mailto:philippe.ciais@lsce.ipsl.fr)

13 JosepPeñuelas: [josep.penuelas@uab.cat](mailto:josep.penuelas@uab.cat)

14 Nicolas Viovy: [nicolas.viovy@lsce.ipsl.fr](mailto:nicolas.viovy@lsce.ipsl.fr)

15

16 Running title: Defining the optimal number of PFTs using a trait-based approach.

17

18 Key words: PFT classification, conifers, forest management, functional traits, hierarchical
classification.

20

21 **Word Counts**: Abstract – 289, Main body – 5263, References – 1631, Biosketch – 57,
22 Figures – 368, Tables – 90.

23 References: 42

24 Figures: 5

25 Tables: 3

26

27

28

29

Final manuscript for Global Ecology and Biogeography
November 2, 2016
30 **ABSTRACT**

31 **Aim**

32 Dynamic global vegetation models (DGVM) use a discretization of forest vegetation based on
33 plant functional types (PFT). The physiological and ecological parameters used to model a
34 given PFT are usually fixed, being defined from point-based observations while model
35 applications are often grid-based. This rigid approach is causing spatial biases in the results of
36 DGVM-simulated productivity and biomass-related variables. We aim to overcome this
37 limitation with a new approach that uses a hierarchical classification of forest PFT parameters
38 from traits retrieved from the literature and from the TRY global database of plant traits. This
39 approach is applied to temperate conifers in the DGVM ORCHIDEE-FM, which has previously
40 been shown to produce systematic biases in the simulation of biomass and biomass increments.

41 **Location**

42 Temperate coniferous forests; France

43 **Time period**

44 2005-2012

45 **Major taxa studied**

46 Gymnosperms

47 **Methods**

48 The five major coniferous species in France, *Abies alba*, *Picea abies*, *Pinus pinaster*, *Pinus sylvestris* and *Pseudotsuga menziesii*, were grouped objectively into PFTs within the
49 ORCHIDEE-FM DGVM using a hierarchical classification based on 12 key attributes related
50 to photosynthesis, phenology and allometric relationships.
51
## 52 Results

We show that the single PFT covering all temperate coniferous forest that was used by default in ORCHIDEE-FM could be replaced by two representative sub-categories defined by grouping species-level data, without necessarily having to adopt a set of parameters for each species. The definition of new temperate conifer PFTs with this approach allows us to reduce the spatial heterogeneity by 40% in average in model-measurement misfit for stand volume, growth and stand density at the regional scale.

## 59 Main conclusions

The proposed approach to improve the representation of plant functional types in DGVMs, while keeping the number of different PFTs manageable, is promising for application to regions where a single PFT can correspond to a number of different species.
63 INTRODUCTION

64 Dynamic global vegetation models (DGVMs) simulate the biogeochemical and hydrological
65 processes involved in ecosystem functioning. They model the controls on the fluxes of carbon,
66 energy and water from local up to regional and global scale. DGVMs also simulate the response
67 of processes and state variables to environmental drivers. Most of these models use a discretized
68 representation of vegetation using Plant Functional Types (PFT), with generally a dozen
69 different PFTs representing the world's vegetation. By definition, each PFT represents a
70 "mean" type of vegetation defined by a set of traits (translated to model parameters). While the
71 use of a limited number of PFTs is practical to for grid-based simulations at large scale, this
72 rather rigid representation of the vegetation does not allow a robust analysis of regional
73 gradients in productivity and biomass, and leads to over-simplify processes controlling
74 vegetation dynamics. It is therefore a source of uncertainty in model predictions (Alton, 2011).

75 The research question of how to deal with the uncertainties resulting from the discrete and rigid
76 representation of PFTs in DGVMs has been addressed by increasing the number of PFTs, or by
77 introducing a continuous spatial variation of vegetation parameters, e.g. using plant trait
78 datasets (Kleidon & Mooney, 2000; Reu et al., 2011; Van Bodegom et al., 2012; Wullschleger
79 et al., 2014; Verheijen et al., 2015). A continuous parameterization of parameters allows spatial
80 variability of parameters within a given PFT (Van Bodegom et al., 2012; Pavlick et al., 2012;
81 Verheijen et al., 2013). However, this method needs a large amount of data for model
82 calibration. Although plant trait databases are building up an increasingly large bank of relevant
83 data, the quality and spatial representativeness of these data might not yet be sufficient to derive
84 robust relationships at the global scale, and even less so at regional scale, such as in our case
85 for temperate coniferous forest in Europe.

86 In this context, increasing the number of PFTs seems a reasonable compromise. Naudt et al.
87 (2015) defined one PFT by dominant European tree species. This strategy is relevant for forests
in Europe where only a few dominant tree species prevail, but it is impractical in higher
biodiversity systems such as grasslands or tropical forests, in which intermediate levels of PFT
representation are needed.

Past studies mostly focused on biogeographical or climatic envelopes to define new PFTs for
DGVMs (Nemani & Running, 1996; Bonan et al., 2002; Harrison et al., 2010) and often sought
to find an optimal and/or minimal subset of traits to give the best classification (Pillar, 1999;
Pausas & Lavorel, 2003; Pillar & Sosinski, 2003; Bernhardt-Römermann et al., 2008).
Boulangeat et al. (2012) already argued for the necessity to improve PFTs in DGVMs by
regrouping representative species with a classification based on key biological characteristics.

Here, we propose an approach to refine the PFT representation in typical DGVMs through the
case study of improving the representation of the spatial variability of the parameters of the
PFT "temperate needleleaf forest" (TNF) across Europe, using the ORCHIDEE-FM DGVM.
The ORCHIDEE-FM DGVM ("Organising Carbon and Hydrology In Dynamic Ecosystems –
Forest Management"; Krinner et al., 2005; Bellassen et al., 2010) simulates the growth, and the
competition between individuals leading to density-driven mortality within a forest stand, and
optionally an idealized management practice of plantation forests, assumed to be even-aged.
The results from this DGVM was evaluated for temperate broadleaf and needleleaf forest in
Europe against, among other data sources, biometric data obtained from numerous plots from
national forest inventories by Bellassen et al., (2011). Although simulations of biomass and
growth were satisfactory for broadleaf forests at the regional scale, the model was not able to
reproduce regional patterns for temperate needleleaf forests. Bellassen et al. (2011) found
regional biases in conifer forests standing volume and growth as a function of forest age, which
could not be corrected by simply optimizing parameters, unlike for broadleaf forests
(Supplementary Fig.S1). They concluded that the main reason for this model systematic error
was the use of a unique conifer PFT associated with a single set of parameters. This suggests
113 that the model suffers more weakness due to its too coarse discretization of vegetation, than
114 from its internal structure and parameter values.

115 Therefore, in this study, we seek to establish if TNF species can be objectively sub-divided into
116 sub-groups with coherent set of parameters for inclusion in the ORCHIDEE FM DGVM. We
117 search for coherent parameters related to eco-physiological attributes, namely photosynthesis,
118 phenology, allometry and differences in management practice.

## 119 METHODS

120 We proposed a new approach for improving the representation of PFTs in ORCHIDEE-FM
121 based on literature review and field data analysis. This method combines three steps:

122 1) We first identified 12 key sensitive eco-physiological parameters, i.e., parameters that
123 have important influence on the model output. These parameters can be related to
124 photosynthesis, phenology, allometry and management parameters for TNF forests. The
125 definition of sensitive is based on previous work that optimized parameters against net
126 ecosystem exchange (NEE) CO2 and evaporation flux measurements at TNF sites with
127 eddy-covariance technique.

128 2) Next, we retrieved a range of values for each sensitive parameter obtained in step-1,
129 based on a literature review and the TRY global trait database (Kattge et al., 2011). This
130 search was restricted to the five most abundant TNF species present in France (Vallet
131 et al., 2006): Silver fir (Abies alba), Norway spruce (Picea abies), Maritime pine (Pinus
132 pinaster), Scots pine (Pinus sylvestris) and Douglas fir (Pseudotsuga menziesii).

133 3) Based on the results of 1) and 2), we showed that TNF species can be objectively
134 regrouped using a hierarchical classification into two homogeneous clusters that have
135 parameter similarities. Implemented in the ORCHIDEE model, we showed that the new
136 PFTs improve the spatial representation of conifer stands. The consistency between
137 PFTs based on parameter values from the literature and PFTs used by ORCHIDEE was
138 assessed by comparing the parameter values of each cluster defined with independent
139 photosynthesis and phenological data from the literature to the values obtained from the
140 calibration of ORCHIDEE using eddy covariance (EC) evaporation and NEE
141 measurements (Santaren et al., 2007; Kuppel et al., 2014).

142 **Description of the ORCHIDEE-FM model**

143 ORCHIDEE is a process-based vegetation model that can be run point-based or grid-based
144 (Krinner et al., 2005). In the standard version of ORCHIDEE (v.1.9.6) an average forest is
145 simulated as a mature stand in near equilibrium between growth and mortality. For given
146 vegetation, soil type and climatic condition, ORCHIDEE simulates the carbon, water and
147 energy budgets on a half-hourly time step (Fig. 1).

148 Bellassen *et al.* (2010) added a forest management module (FM; Fig. 1) to ORCHIDEE
149 (v.1.9.6) to simulate self-thinning and management. The FM module replaces an “average tree”
150 representation with a statistically distributed population of individual trees and includes rules
151 for both natural mortality and human interventions. This forest management version of the
152 model is called ORCHIDEE-FM (v.1.9.6-rev.30). It simulates carbon stocks and detailed stand
153 structure of forests of varying ages (Bellassen *et al.*, 2011) which can be directly compared with
154 forest inventory data or ecological site biometric observations. The FM equations are fully
155 described by Bellassen *et al.* (2010).

156

157 **Selection of sensitive eco-physiological parameters**

158 The starting point for improving the description of TNF forests is the determination of the
159 variability of so called sensitive or influential eco-physiological parameters among the TNF
species. Among all the parameters of ORCHIDEE-FM, we defined as influential those whose variations have a significant effect on the modelled stand biomass. The selection of influential parameters was carried out by comparing information from previous sensitivity tests made by Bellassen *et al.* (2010) to the results of an optimization of ORCHIDEE eco-physiological parameters made by Kuppel (2012) using eddy-covariance measurements. Both methods allow the identification of sensitive parameters, given different “target” state variables. Bellassen *et al.* (2010) identified “slow” parameters that control stand density and biomass dynamics on time scales of decades (hereafter B-parameters), while Kuppel (2012) identified mainly “fast” parameters that control physiological processes controlling fluxes (hereafter P-parameters), on time scales of days to seasons. Both fast and slow parameters interact with each other in determining ecosystem carbon state variables in the model equations (Fig. 1). This interaction is non-symmetrical, however. Slow parameters controlling long-term forest growth and carbon stocks have a small impact on hourly to seasonal fluxes, whereas fast parameters do influence decadal net primary productivity (NPP) and thus a change from their default value is expected to modify the decadal stand-scale growth, competition and mortality, and thus the long-term changes of biomass.

The method used to select each parameter is fully described in Appendix S2 (see Table 1 for parameters names). The results are summarized here. Three B-parameters were selected as being influential: allometric coefficients that govern the tree biomass-to-circumference ratios ($b_{bc}$), self-thinning parameters ($\ln(\alpha_{st})$ and $\beta_{st}$) and the threshold of stand density at which forest management thinning is set to occur (*rdi_target*) for managed forests. Six P-parameters were found influential: the maximum leaf-level limited carboxylation rate normalized at 25°C ($V_{cmax_opt}$), the optimal temperature for photosynthesis ($T_{opt}$), the slope of the Ball-Berry linear relationship between stomatal conductance ($g_{slope}$) and net assimilation rate, the specific leaf area (SLA) and the critical leaf age for leaf senescence ($L_{age,crit}$).
Processes in which the different parameters identified in this study are involved are illustrated in Fig. 1 and listed in Table 1. The default values of these parameters for the TNF PFT in ORCHIDEE are as defined by Krinner *et al.* (2005) and by Bellassen *et al.* (2010) and listed in Table 2. Optimized values obtained from inversion of the fast parameters against EC measurements are from Santaren *et al.* (2007) and Kuppel *et al.* (2014).

## 190 Literature review of key parameters

We gathered species-specific data for the values of each sensitive parameter listed in Table 1. A search of the literature (up to 2012) was conducted for the predominant TNF species found in France. In addition of this literature search, a number of values were obtained through the TRY database (Kattge *et al.*, 2011), especially for SLA. The number of values found for each parameter is indicated in Fig. 2 and Table 2. The literature was first searched for each species and each parameter; the search was then restricted to temperate regions (Köppen classification; Peel *et al.*, 2007). TRY values that were also in the literature were eliminated to avoid duplication; and when multiple measurements were available at a given site, the median value of the reported values was used. Care was taken to ensure the physiological meaning of the data found in the TRY database corresponded to the model parameter definitions. The use of mean values reduced the size of the dataset and prevented any bias from over-sampling a few better-measured sites.

For all parameters identified as influential, we tried to minimize the uncertainty due to different measurement methods by filtering or correcting data whenever possible. Values were not corrected according to environmental variability, nor species provenance or age. All references were searched for each parameter described in the following sections and a detailed description of the methodology used to homogenize values for each parameter can be found in Appendices B and C in the Supporting Information.
Even though it is a sensitive parameter according to Kuppel et al. (2012), the number of data reviewed for g<sub>slope</sub> was insufficient to perform a separate analysis per species within the TNF PFT (Fig. 2). We thus decided to extract another parameter from the literature, the maximal stomatal conductance (g<sub>smax</sub> – even though it is not directly used in ORCHIDEE). We defined g<sub>smax</sub> as the maximum conductance measured at one site at any measurement period. The model is sensitive to the value of g<sub>slope</sub>, but because we have no evidence that g<sub>smax</sub> is also a key parameter, the following analyses were performed either including or excluding g<sub>smax</sub>.

216

## Clustering of species based on parameters values

For all parameters, except rdi_target, T<sub>opt</sub> (see Table 1 for definitions) and g<sub>slope</sub>, because of insufficient data, the inter-species and intra-species parameter variabilities were calculated using Kruskal–Wallis one-way analysis of variance and the Mahalanobis distance (distance between species taking into account the parameter variance; Mahalanobis, 1936):

$$D(x) = \sqrt{(x - \mu)^T \Sigma^{-1} (x - \mu)} \quad (1)$$

where x is a vector of data (here different parameters) with p columns, μ the mean of the distribution of length p, Σ the covariance matrix (p*p) of the distribution. These methods have the advantage of providing a statistical test of similarity. Their drawback is the impossibility of obtaining a unique set of species clusters based on the variability of all the parameters.

As all parameters are quantitative, we used a clustering algorithm to objectively define homogeneous groups of species based on Ward's method (1963). In this method, all parameter median values are centred and scaled to calculate dissimilarity (distance, eq.2).
$$D(C1, C2) = \sqrt{\frac{n1n2}{n1+n2}} d(\mu1, \mu2) \quad (2)$$

229 Where $D(C1, C2)$ is the distance between clusters $C1$ and $C2$, $n1$ and $n2$ the number of data, and
230 $d(\mu1, \mu2)$ the distance between their respective centroids.

231 Data treatment was entirely performed with R software (R Team 2011).

232 **ORCHIDEE-FM Simulations**

233 First, we analysed previous regional simulations from Bellassen *et al.* (2011) in the light of the
234 new PFT classification by comparing model/measurement errors from Bellassen *et al.*, (2011)
235 to the distribution of the newly defined PFTs.

236 In addition, four sets of simulations were performed to assess the gain (or loss) of information
237 obtained by different sub-divisions of TNF species, each based on coherent parameter values.

238 The first set (REF) uses the default parameter values and processes of ORCHIDEE-FM as
239 described by Bellassen *et al.* (2010), that is parameters are the same for all TNF forests. The
240 second configuration (S1) use one PFT to represent TNF forest but this PFT is parameterized
241 according to results of the literature review. The third set (S2) was performed using the new
242 PFTs corresponding to species groups obtained from the hierarchical clustering of parameters.

243 The fourth simulation (S3) was performed with the parameters of each TNF species (5 in total),
244 each species being modeled as a separate PFT. All the simulations are evaluated against the
245 French NFI observations interpolated at a 0.25° spatial resolution. We simulated 50 years old
246 stands. To minimize bias between simulations and observations, we averaged the simulation
247 outputs over a 10 year period (between years 45 and 55) and only trees with a diameter higher
248 than 7.5 cm were used.

249 The same protocol and pedo-climatic inputs as in Bellassen *et al.* (2011) were used for all the
250 simulations. No historical change in species distribution or forest area was taken into account
251 in the simulations. The performance of each simulation was assessed against observations by
252 using the Spearman's rank correlation coefficient for three different variables: the average stand
253 aboveground volume, volume increment and stand density.

254

255 **RESULTS**

256 **Synthesis of parameter values**

257 We found differences in the quantity of data collected for the different parameters. For instance,
258 336 values were collected for SLA while only one value per species could be obtained for T<sub>opt</sub>
259 and rdi_target (Table 3). There is also a range in the number of data available per species. The
260 values for *Pinus sylvestris* and *Picea abies* are more numerous than for *Abies alba*, *Pinus*
261 *pinaster* and *Pseudotsuga menziesii* in particular for V<sub>cmax</sub>, SLA and g<sub>smax</sub>. The mean coefficient
262 of variation (CV) for parameters V<sub>cmax</sub>, SLA and g<sub>smax</sub> are relatively high (0.35, 0.43 and 0.71,
263 respectively). This spread of parameter distributions could be the result of environmental
264 conditions, bias in measurement method or representation.

265 **Differences between species**

266 Species was found to be a statistically significant explanatory variable for differences in P-
267 parameters. According to the Kruskal-Wallis one-way analysis of variance, the values of
268 V<sub>cmax25</sub>, J<sub>max</sub>/V<sub>cmax</sub>, L<sub>age,crit</sub> and SLA were found to be significantly different between species (p-
269 value<0.05). Variance analysis shows a significant inter-species variability for most P-
270 parameters, whereas B-parameter differences are statistically not significant between species.
271 For parameters V<sub>cmax</sub>, J<sub>max</sub>/V<sub>cmax</sub>, SLA and Longevity, the Mahalanobis distance between species
272 is found to be higher than for g<sub>smax</sub> (Fig.S2 in Appendix S1). However for B-parameters, the
273 distance between species is relatively lower than for P-parameters. All *Pinus* species fall into
274 the same group for P-parameters, which supports grouping these species into one group for
275 modelling carbon stocks and fluxes.

276 **Differences in parameters between literature values and eddy covariance optimized values**
277

278 We observe discrepancies between literature values and values optimized from eddy covariance
279 data for several parameters. $V_{cmax25}$ values observed in the literature for *Pinus pinaster* and *Pinus
280 sylvestris* (90.31 and 106.61µmol m⁻² s⁻¹) are higher than optimized values obtained by Kuppel
281 *et al.* (2014) from 20 conifer sites (viz., 41.7 and 31.9 µmol m⁻² s⁻¹; Table 3). By contrast, for
282 non-pine species, $V_{cmax25}$ estimated from literature data (at leaf scale) are found to be within the
283 range of the results from Kuppel *et al.* (2014). For $T_{opt}$, only values obtained for *Pinus sylvestris*
284 are close to the values optimized by Kuppel *et al.* (2014), but the median value considering all
285 TNF species is twofold higher in the literature data than in Kuppel *et al.* (2014) (32.3°C
286 compared to 17.5°C). SLA literature values are half those of Kuppel *et al.* (2014) (10 mm² mg⁻¹),
287 excepted for *Pseudotsuga menziesii* for which we found an SLA of 7.08 mm² mg⁻¹. For
288 $L_{age,crit}$, only the literature value for *Pinus sylvestris* (34.0 months) is close to Kuppel *et al.*
289 (2014) (35.6 months). Too few data were available to compare $g_{slope}$ for all species, however
290 the single literature value reviewed for *Pinus pinaster* (10.05) is in the range of the inversion
291 results of Kuppel *et al.* (2014) (8.8) and Santaren *et al.* (2007) (8.1). The values of $g_{slope}$ obtained
292 from the literature for *Picea abies* (5.94) and *Pinus sylvestris* (5.59) are smaller than values
293 found by the Bayesian calibration of Kuppel *et al.* (2014). Results from the Santaren *et al.*
294 (2007) and Kuppel *et al.* (2014) optimizations are summarized in Table 3.
## Differences between literature and values derived from forest inventory data for B-parameters

In addition to the literature search, B-parameters linked to biomass dynamics ($b_{bc}$, $\ln(\alpha_{st})$ and $\beta_{st}$) were re-determined by fitting to the French national forest inventory (NFI) data (IFN, 2012; refer to Appendix S2 for a complete description). Fitted values of $b_{bc}$ are from this study, while values of $\ln(\alpha_{st})$ and $\beta_{st}$ are from Charru et al. (2012) (also obtained from French NFI). For these three B-parameters, most of the values fitted using the French NFI data were found to be out of the range of the literature values (blue dots in Fig. 2). Values for $b_{bc}$ and $\beta_{st}$ are larger for *Pinus* compared to literature values, but the opposite is found for $\ln(\alpha_{st})$. For *Picea abies*, values of $b_{bc}$ and $\beta_{st}$ are smaller than literature values and larger for $\ln(\alpha_{st})$. Exceptions are *Pseudotsuga menziesii* for which $\ln(\alpha_{st})$, $\beta_{st}$, and $b_{bc}$ are in the same range between the two approaches, and *Abies alba* for parameters $\ln(\alpha_{st})$ and $\beta_{st}$.

## Species clustering

The first clustering of species, based on all parameters, is called global classification (Fig. 3a). It results in a first group composed of *Abies alba*, *Picea abies* and *Pseudotsuga menziesii* and a second group with *Pinus sylvestris* and *Pinus pinaster*. Two potential biases were tested for:
1) parameter redundancy – the inclusion of different parameters related to the same eco-physiological process may give too much weight to this process in the classification – and 2) the type of parameter used for the classification, because P- and B-parameters are related to different processes and may have uneven weights in the classification. To test for these possible biases, three other classifications were performed. A classification based on P-parameters only was found to be equivalent to the global classification, yet with a reduced distance between groups and a longer distance between *Pseudotsuga menziesii* and *Picea abies*, *Abies alba* (Fig. 3b). Another classification based on B-parameters resulted in grouping of *Pseudotsuga*
319 *menziesii* and *Picea abies* on one side, *Pinus sylvestris*, *Abies alba* and *Pinus pinaster* on the
320 other side (Fig. 3c).

321 To decrease redundancy, a second classification was based on all parameters except ln(αst), Optj
322 and Jmax/Vmax to minimize redundancy of these parameters with βst, Optv and Vcmax25 respectively
323 (Figs 3d and 3e). The two main clusters of this second classification are found to be the same
324 as in the (first) global classification, and appear even more clearly since the distances between
325 species within each cluster are relatively smaller than in the global classification. If we separate
326 again between the reduced parameter set for the P-and B-parameters (hereafter rP and rB), the
327 classification result using rB-parameters only is equivalent to the one using rP-parameters only.
328 From this we conclude that two new PFTs describing temperate conifers can be defined, namely
329 PFT1 with *Abies alba/Picea abies/Pseudotsuga menziesii* and PFT2 with *Pinus pinaster/Pinus
330 sylvestris*.

331 **Model evaluation**

332 We reprocessed the results of Bellassen et al. (2011) with the new classification. The spatial
333 distribution of the new PFTs within the TNF group is illustrated in the Supplementary
334 Information (Fig. S3). It allows us to explain qualitatively the spatial distribution of differences
335 between observed and simulated volumes and volume increments (Fig. 4). From the
336 observations available in the French NFI dataset we separated pixels where only one PFT was
337 present from pixels where the two PFTs coexist. Pixels which show an overestimation of the
338 volume increment correspond to a predominance of PFT2 (Fig. 4a); and the opposite, pixels
339 showing an underestimation of the volume increment correspond to a predominance of PFT1.
340 This distinction is also visible on pixels where the two PFTs are present (Fig. 4b).

341 Figure 5 illustrates the spatial correlation between observed and simulated volume, volume
342 increment and stand density with 1, 2 or 5 PFTs. Considering TNF as a single PFT with
parameters values retrieved from the literature (S1) degrades the spatial representation of volume and volume increment compared to the reference simulation (REF). Increasing the number of PFTs improves the spatial correlation (r) between observed and simulated total aboveground volume ($r_{REF}=0.27$, $r_{S2}=0.73$), volume increment ($r_{REF}=0.44$, $r_{S2}=0.56$) and average stand density ($r_{REF}=-0.16$, $r_{S2}=0.61$). However, further increasing the number of PFTs up to one PFT per species in S3 does not improve the results compared to S2, excepted for the volume increment ($r_{S2}=0.56$, $r_{S3}=0.7$).

## 350 DISCUSSION

### 351 The species clustering method

352 The main result of this study is that the temperate coniferous species planted in French forests 353 can be grouped into a few distinct PFTs from our clustering analysis method based on *sensitive* 354 parameter similarity (Fig.3). The use of a reduced set of influential parameters, designed to 355 prevent redundancy of information across parameters, leads to a similar grouping of species 356 when considering either P- or B-parameters, or both, to group species. This suggests that the 357 parameter classification method proposed here can determine groups of species. The significant 358 inter-species variability observed for most P-parameters, but not for B-parameters could be 359 partly explained by the fact that some B-parameters related to forestry (allometric parameters 360 or *rdi*) are determined from managed stand data, where management strategy can mask species 361 effects. By contrast, species effects are present in P-parameters.

362 Moreover, all eco-physiological parameters are driven by environment and evolutionary 363 selection. Inter-specific variability should be statistically distinguishable from environmentally 364 induced variability. For example here, too few data are available to screen environmental effects 365 or to perform robust statistical tests on all parameters. Thus co-linearity between environment 366 and trait values can potentially bias the classification. However, most PFTs in DGVMs are
currently parameterized without screening out possible environmental effects, and their
parameterization is usually based on far fewer references than the 156 data used here (Zaehle
et al., 2005; Alton, 2011; Wullschleger et al., 2014). The method used in this study is thus
relevant for use with DGVMs other than ORCHIDEE-FM.

Past studies sought to define PFTs by finding an optimal set of traits (Pillar & Sosinski, 2003;
Bernhardt-Römermann et al., 2008; Boulangeat et al., 2012). In contrast, in this study we started
from the modelling viewpoint and sought to define PFTs by selecting key functional traits
related to parameter values. We managed to define two different PFTs from TNF: one group is
formed by pines: *Pinus pinaster* and *Pinus sylvestris*, and the other by spruces and firs: *Abies
alba*, *Picea abies* and *Pseudotsuga menziesii*. We can see that the classification results in the
separation of species with *pinus* trees clustered together, arguing for a genus-based
classification. Liston *et al*. (1999) already demonstrated clear phylogenetic differences between
pines and spruces. However, all eco-physiological parameters are driven by environment and
strong selections have been operating for a long time on these species. Thus, the grouping of
pines in a same cluster may be the result of parallel acclimatization of species sharing the same
bio-climatic conditions.

This proposed re-classification of the TNF of France into two PFTs seems to explain the spatial distribution of simulated misfits of the volume and volume increment in Bellassen *et al.* (2011) (Fig.4, Sup.Fig.3), even within pixels where the two PFTs are present (Fig. 4b). Moreover, the use of two PFTs improves the spatial correlation between observed and simulated stand variables by more than 40 %, while the use of one PFT by species does not bring more precision.
This result tends to confirm our method and that the differences between simulated and observed stands are mostly attributable to the PFT representation rather than a potential bias linked to local soil-climate conditions. Past studies highlighted the discrepancies between the leaf and the canopy scale, but also the structural uncertainties of current DGVMs, which
therefore cannot be parameterized with leaf-scale observations (Bonan et al., 2012). Difficulties in parameterizing one unique PFT can also arise from the heterogeneity of species within this PFT.

## Inconsistencies between parameter values from eddy-covariance and literature data

The inconsistency between the literature review and eddy-covariance (EC) based optimization results from Kuppel et al. (2014) implies that the Bayesian optimization based on variables like fluxes (which depend on many parameters) is inadequate because it does not correctly represent physiological processes. The optimization method may benefit from the growing availability of trait observations in the TRY database. However, as highlighted by Bonan et al. (2012), it also suggests that some fundamental processes in ORCHIDEE-FM are lacking or poorly represented, thus leading to compensations in the parameterization.

We have no information about the correlation of optimized parameter errors related to temperate needleleaf forest in Kuppel et al. (2014). However, Kuppel (2012) showed strong posterior error correlations between the carboxylation rate, the needle lifespan and the optimal temperature for example when considering all multi-site optimizations and all PFTs (see Fig.S13 of Kuppel (2012) for the correlation matrix of parameters errors). These correlations occur because of the inability (or equifinality, i.e. different combinations of parameter values can lead to the same result) of EC data to discriminate the value of each single parameter. Traits optimized in Kuppel et al. (2014) can be biased due to equifinality or a poor representation of some processes, thus leading to an over-parameterization when constrained by EC data.

In the other hand, parameters optimized by Kuppel et al. (2014) represent traits integrated at the canopy level, while measurements gathered in this study were performed on individual leaves. Therefore, the comparison between optimized and measured traits should be made in the light of the different biases specific to each method and scale.
These differences between leaf-scale measurements of parameters from literature data and the values optimized using ecosystem-level observations of fluxes of NEE, H, and E can be attributed to the lack of physiological consideration in the mathematical procedure used to optimize parameters from ecosystem-scale flux data. Despite the constraint imposed by setting a realistic prior error value for each parameter in the inversion (based on expert judgment), the minimization of a Bayesian cost function produces a posterior value based on purely mathematical criteria. Such a value may work well locally, but might be disconnected from eco-physiological reality (Pelikan 2005; Ward et al., 2010). This discrepancy can arise from the erroneous representation of eco-physiological processes in ORCHIDEE, but also compensations or co-variations between parameters dealing with same eco-physiological processes that cannot be separated by NEE measurements alone. Another point revealed by this study is the potential bias introduced by giving equal weight to all the geographic sites used for optimization. Most of the sites (20 in total) optimized in Kuppel (2012) were *Picea abies* and *Pinus sylvestris* whereas the only site optimized by Santaren et al. (2007) was *Pinus pinaster*. In the end, parameter values of optimizations were allocated to the same PFT, ignoring the different species. However we show here differences between species for key parameters inside the same plant functional type. The problem highlighted here is an over-representation of one species in the optimization. Unless this over-representation reflects the real species distribution on the ground, it will bias the final value retained for the global PFT. Especially in the Northern Hemisphere, where we expect an increase in conifer plantations to meet the objectives of climate change mitigation, the non-representation of the different capacity of species will lead to the wrong estimation of carbon stocks and fluxes (Vallet et al., 2009; Matthews et al., 2014).

Our results thus highlight the importance of confronting the results of “brute force” parameter calibration methods using mathematics and EC data with a more refined approach based on eco-physiological knowledge before using them for model improvement.
The use of a literature review as performed in this study can provide better optimizations by applying appropriated ranges of variability for each parameter. One of the limitations of a literature search however, is the lack of data harmonization related to differences in environmental conditions and methodologies. We tried to correct for some methodological effects but too few data were available to allow a correction based on biotic and abiotic environmental as well as management factors, which leads to large spread in parameter estimates (Fig. 2). Results obtained for each parameter and the associated uncertainties are debated in detail in the extended discussion (Appendix S4).

## CONCLUSION & PERSPECTIVES

In this study we proposed a new method combining literature review and hierarchical classification to constrain sensitive parameters in a DGVM model. The model used was ORCHIDEE FM and our results are model-specific, but the approach is applicable to any other DGVM with a discrete representation of PFTs. We proposed to group species into coherent groups based on parameter similarity for refining the number of PFTs. This method was applied to temperate conifers in France for 5 dominant species, and resulted in the successful derivation of two PFTs with coherent parameters. The parsimonious choice of two new PFTs to describe TNF improved the spatial representation of conifers stands carbon dynamics by more than 40% compared to the reference simulation where all TNF have the same parameter values. However using a more refined representation with PFT per species did not improve the simulations, and did not degrade them either. Our approach of coherent parameter classification is thus a good compromise between the too rigid simple PFT representation and a more complicated fully trait-based continuous parameterization of vegetation parameters in DGVMs. Thanks to the increasing availability of trait observation in global databases, our method could be applied to
refine other PFTs, for instance the “temperate broadleaf forests” PFT which currently groups 465 together in a model like ORCHIDEE species with different characteristics such as fast growing 466 pioneer species (*Populus* sp.) and slow growing dominant species like oak (*Quercus* sp.) 467

In addition, we showed clear discrepancies between observed parameter values and values 468 optimized against eddy covariance flux measurement. This result highlights the importance of 469 confronting the results of “brute force” parameter calibration methods with a more refined eco- 470 physiological knowledge. We argue that this warrants further research for a better 471 understanding of the definition of parameters, their characterization and the adequacy of point- 472 based field observations to parameterize grid-based models. 473

474

475 Acknowledgements

The authors would like to acknowledge the financial support from the European Research 476 Council Synergy grant ERC-SyG-2013-610028 IMBALANCE-P. We thank the anonymous 477 reviewers and the associate editor for their constructive comments that improved both the 478 readability and the robustness of the manuscript. This work was granted access to the HPC 479 resources of TGCC under the allocation 2015-6328 made by GENCI. The authors thank Jens 480 Kattge, Ian Wright, Peter Reich, Hans Cornelissen, Bill Shipley, Belinda Medlyn, Eric Garnier, 481 Jenny Ordonez and William Cornwell for their permission to use the TRY data in this study. 482

483

484 Data Accessibility

All data and sources used in this study are listed in Appendix S3.

486

487
488 **Biosketch**

489 M. Peaucelle did his PhD at the Climate and Environment Science Laboratory (LSCE, France).

490 His research focuses on improving the representation of temperate forest in the ORCHIDEE

491 DGVM. In particular, he is interested in improving the representation of plant functional types

492 by a more continuous approach exploring different ways to include the plasticity of plant

493 functional traits.
494 Supporting information

495 Appendix S1: Supplementary tables and figures (list below)

496 Table S1a: Optimized parameters from Kuppel (2012).

497 Table S1b: Parameters used in the sensitivity analysis from Bellassen *et al.*, (2010).

498 Table S2: Data sources for French sylvicultural reference systems.

499 Figure S1: Representation of the relative difference between observed and simulated stand volume and volume increment.

501 Figure S2: Representation of the Mahalanobis' distance between each species for photosynthesis and biomass parameters.

503 Figure S3: Spatial distribution of the two new PFTs over France.

504 Figure S4: Evolution of the relative density index against the dominant height over a rotation for a stand of *Picea abies*

506 Appendix S2: Detailed methodology for the literature review, the homogenization of values for each parameter.

508 Appendix S3: List of references retrieved for each parameter.

509 Appendix S4: Extended discussion about uncertainties associated with each parameter.
## 510 References

511 Alton, P.B. (2011) How useful are plant functional types in global simulations of the carbon, water, and energy cycles? *Journal of Geophysical Research: Biogeosciences*, 116.

513 Bellassen, V., Le Maire, G., Dhôte, J.F., Ciais, P. & Viovy, N. (2010) Modelling forest management within a global vegetation model—Part 1: Model structure and general behaviour. *Ecological Modelling*, 221, 2458–2474.

516 Bellassen, V., le Maire, G., Guin, O., Dhôte, J.F., Ciais, P. & Viovy, N. (2011) Modelling forest management within a global vegetation model—Part 2: Model validation from a tree to a continental scale. *Ecological Modelling*, 222, 57–75.

519 Bernhardt-Römermann, M., Römermann, C., Nuske, R., Parth, A., Klotz, S., Schmidt, W. & Stadler, J. (2008) On the identification of the most suitable traits for plant functional trait analyses. *Oikos*, 117, 1533–1541.

522 Van Bodegom, P.M., Douma, J.C., Witte, J.P.M., Ordoñez, J.C., Bartholomeus, R.P. & Aerts, R. (2012) Going beyond limitations of plant functional types when predicting global ecosystem–atmosphere fluxes: exploring the merits of traits-based approaches. *Global Ecology and Biogeography*, 21, 625–636.

526 Bonan, G.B., Levis, S., Kergoat, L. & Oleson, K.W. (2002) Landscapes as patches of plant functional types: An integrating concept for climate and ecosystem models. *Global Biogeochemical Cycles*, 16, 5–23.

529 Bonan, G.B., Oleson, K.W., Fisher, R.A., Lasslop, G. & Reichstein, M. (2012) Reconciling leaf physiological traits and canopy flux data: Use of the TRY and FLUXNET databases in the Community Land Model version 4. *Journal of Geophysical Research: Biogeosciences*, 117.

532 Boulangeat, I., Philippe, P., Abdulhak, S., Douzet, R., Garraud, L., Lavergne, S., Lavorel, S., Van Es, J., Vittoz, P. & Thuiller, W. (2012) Improving plant functional groups for dynamic
models of biodiversity: at the crossroads between functional and community ecology. *Global Change Biology*, 18, 3464–3475.

Charru, M., Seynave, I., Morneau, F., Rivoire, M. & Bontemps, J.-D. (2012) Significant differences and curvilinearity in the self-thinning relationships of 11 temperate tree species assessed from forest inventory data. *Annals of Forest Science*, 69, 195–205.

Harrison, S.P., Prentice, I.C., Barboni, D., Kohfeld, K.E., Ni, J. & Sutra, J.-P. (2010) Ecophysiological and bioclimatic foundations for a global plant functional classification. *Journal of Vegetation Science*, 21, 300–317.

IFN. 2012. “Raw Inventory Data.” http://inventaire-forestier.ign.fr/.

Kalnay, E., Kanamitsu, M., Kistler, R., Collins, W., Deaven, D., Gandin, L., Iredell, M., Saha, S., White, G., Woollen, J., Zhu, Y., Leetmaa, A., Reynolds, R., Chelliah, M., Ebisuzaki, W., Higgins, W., Janowiak, J., Mo, K.C., Ropelewski, C., Wang, J., Jenne, R. & Joseph, D. (1996) The NCEP/NCAR 40-Year Reanalysis Project. *Bulletin of the American Meteorological Society*, **77**, 437–471.

Kattge, J., Díaz, S., Lavorel, S., Prentice, I.C., Leadley, P., Bönisch, G., Garnier, E., Westoby, M., Reich, P.B., Wright, I.J., Cornelissen, J.H.C., Violle, C., Harrison, S.P., Van BODEGOM, P.M., Reichstein, M., Enquist, B.J., Soudzilovskaia, N.A., Ackerly, D.D., Anand, M., Atkin, O., Bahn, M., Baker, T.R., Baldocchi, D., Bekker, R., Blanco, C.C., Blonder, B., Bond, W.J., Bradstock, R., Bunker, D.E., Casanoves, F., Cavender-Bares, J., Chambers, J.Q., Chapin Iii, F.S., Chave, J., Coomes, D., Cornwell, W.K., Craine, J.M., Dobrin, B.H., Duarte, L., Durka, W., Elser, J., Esser, G., Estiarte, M., Fagan, W.F., Fang, J., Fernández-Méndez, F., Fidelis, A., Finegan, B., Flores, O., Ford, H., Frank, D., Freschet, G.T., Fyllas, N.M., Gallagher, R.V., Green, W.A., Gutierrez, A.G., Hickler, T., Higgins, S.I., Hodgson, J.G., Jalili, A., Jansen, S., Joly, C.A., Kerkhoff, A.J., Kirkup, D., Kitajima, K., Kleyer, M., Klotz, S., Knops, J.M.H.,
Kramer, K., Kühn, I., Kurokawa, H., Laughlin, D., Lee, T.D., Leishman, M., Lens, F., Lenz, T., Lewis, S.L., Lloyd, J., Llusià, J., Louault, F., Ma, S., Mahecha, M.D., Manning, P., Massad, T., Medlyn, B.E., Messier, J., Moles, A.T., Müller, S.C., Nadrowski, K., Naeem, S., Niinemets, Ü., Nöllert, S., Nüske, A., Ogaya, R., Oleksyn, J., Onipchenko, V.G., Onoda, Y., Ordoñez, J., Overbeck, G., Ozinga, W.A., Patiño, S., Paula, S., Pausas, J.G., Peñuelas, J., Phillips, O.L., Pillar, V., Poorter, H., Poorter, L., Poschlod, P., Prinzing, A., Proulx, R., Rammig, A., Reinsch, S., Reu, B., Sack, L., Salgado-Negret, B., Sardans, J., Shiodera, S., Shipley, B., Siefert, A., Sosinski, E., Soussana, J.-F., Swaine, E., Swenson, N., Thompson, K., Thornton, P., Waldram, M., Weiher, E., White, M., White, S., Wright, S.J., Yguel, B., Zaehle, S., Zanne, A.E. & Wirth, C. (2011) TRY – a global database of plant traits. *Global Change Biology*, 17, 2905–2935.

Kleidon, A. & Mooney, H.A. (2000) A global distribution of biodiversity inferred from climatic constraints: results from a process-based modelling study. *Global Change Biology*, 6, 507–523.

Krinner, G., Viovy, N., Noblet-Ducoudré, N. de, Ogée, J., Polcher, J., Friedlingstein, P., Ciais, P., Sitch, S. & Prentice, I.C. (2005) A dynamic global vegetation model for studies of the coupled atmosphere-biosphere system. *Global Biogeochemical Cycles*, 19.

Kuppel, S. (2012) Assimilation de mesures de flux turbulents d'eau et de carbone dans un modèle de la biosphère continentale, Versailles-St Quentin en Yvelines.

Kuppel, S., Peylin, P., Chevallier, F., Bacour, C., Maignan, F. & Richardson, A.D. (2012) Constraining a global ecosystem model with multi-site eddy-covariance data. *Biogeosciences*, 9 (10), 3757–3776.

Kuppel, S., Peylin, P., Maignan, F., Chevallier, F., Kiely, G., Montagnani, L. & Cescatti, A. (2014) Model-data fusion across ecosystems: from multi-site optimizations to global simulations. *Geoscientific Model Development*, 7 (6), 2581–2597.
Liston, A., Robinson, W.A., Piñero, D. & Alvarez-Buylla, E.R. (1999) Phylogenetics of *Pinus* (*Pinaceae*) based on nuclear ribosomal DNA internal transcribed spacer region sequences. *Molecular Phylogenetics and Evolution*, 11, 95-109.

Mahalanobis, P. (1936) On the generalised distance in statistics. *Proceedings National Institute of Science*, 2, 49-55.

Matthews, R., Sokka, L., Soimakallio, S., Mortimer, N., Rix, J., Schelhaas, M.-J., Jenkins, T., Hogan, G., Mackie, E., Morris, A. & Randle, T. (2014) Review of literature on biogenic carbon and life cycle assessment of forest bioenergy. Final Task 1 report, The Research Agency of the Forestry Commission, Farnham, UK.

Naudts, K., Ryder, J., McGrath, M.J., Otto, J., Chen, Y., Valade, A., Bellasen, V., Berhongaray, G., Bönisch, G. & Campioli, M. (2015) A vertically discretised canopy description for ORCHIDEE (SVN r2290) and the modifications to the energy, water and carbon fluxes. *Geoscientific Model Development*, 8, 2035-2065.

Nemani, R. & Running, S.W. (1996) Implementation of a hierarchical global vegetation classification in ecosystem function models. *Journal of Vegetation Science*, 7, 337-346.

Pausas, J.G. & Lavorel, S. (2003) A hierarchical deductive approach for functional types in disturbed ecosystems. *Journal of Vegetation Science*, 14, 409-416.

Pavlick, R., Drewry, D., Bohn, K., Reu, B. & Kleidon, A. (2012) The Jena Diversity-Dynamic Global Vegetation Model (JeDi-DGVM): a diverse approach to representing terrestrial biogeography and biogeochemistry based on plant functional trade-offs. *Biogeosciences*, 10, 4137-4177.

Peel, M.C., Finlayson, B.L. & McMahon, T.A. (2007) Updated world map of the Köppen-Geiger climate classification. *Hydrology and Earth System Sciences*, 11 (5), 1633-1644.
604 Pelikan, M. (2005) *Hierarchical Bayesian Optimization Algorithm. *Studies in Fuzziness and Soft Computing*, 31–48.

606 Pillar, V.D. (1999) *On the identification of optimal plant functional types*. *Journal of Vegetation Science*, 10, 631–640.

608 Pillar, V.D. & Sosinski, E.E. (2003) An improved method for searching plant functional types by numerical analysis. *Journal of Vegetation Science*, 14, 323–332.

610 Reich, P.B., Rich, R.L., Lu, X., Wang, Y.-P. & Oleksyn, J. (2014) Biogeographic variation in evergreen conifer needle longevity and impacts on boreal forest carbon cycle projections.

611

612 *Proceedings of the National Academy of Sciences of the United States of America*, 111, 13703–13708.

614 Reu, B., Zaehle, S., Proulx, R., Bohn, K., Kleidon, A., Pavlick, R. & Schmidtlein, S. (2011)
615 The role of plant functional trade-offs for biodiversity changes and biome shifts under scenarios
616 of global climatic change. *Biogeosciences*, 8, 1255–1266.

617 Santaren, D., Peylin, P., Viovy, N. & Ciais, P. (2007) Optimizing a process-based ecosystem model with eddy-covariance flux measurements: A pine forest in southern France. *Global Biogeochemical Cycles*, 21.

620 Team, R.D.C. (2011) R: A Language and Environment for Statistical Computing, Vienna, 621 Austria.

622 Vallet, P., Dhôte, J.-F., Moguédec, G.L., Ravart, M. & Pignard, G. (2006) Development of total aboveground volume equations for seven important forest tree species in France. *Forest Ecology and Management*, 229, 98–110.

625 Vallet, P., Meredieu, C., Seynave, I., Bélouard, T. & Dhôte, J.F. (2009) Species substitution for carbon storage: Sessile oak versus Corsican pine in France as a case study. *Forest Ecology and Management*, 257, 1314–1323.
628 Verheijen, L., Brovkin, V., Aerts, R., Bönisch, G., Cornelissen, J., Kattge, J., Reich, P., Wright, 629 I. & Bodegom, P. van (2013) Impacts of trait variation through observed trait-climate 630 relationships on performance of a representative Earth System model: a conceptual analysis. 631 *Biogeosciences*, 10 (8), 5497-5515.

632 Verheijen, L.M., Aerts, R., Brovkin, V., Cavender-Bares, J., Cornelissen, J.H.C., Kattge, J. & 633 van Bodegom, P.M. (2015) Inclusion of ecologically based trait variation in plant functional 634 types reduces the projected land carbon sink in an earth system model. *Global Change Biology*, 635 21, 3074–3086.

636 Vetter, M., Churkina, G., Jung, M., Reichstein, M., Zaehle, S., Bondeau, A., Chen, Y., Ciais, 637 P., Feser, F., Freibauer, A., Geyer, R., Jones, C., Papale, D., Tenhunen, J., Tomelleri, E., 638 Trusilova, K., Viovy, N. & Heimann, M. (2008) Analyzing the causes and spatial pattern of the 639 European 2003 carbon flux anomaly using seven models. *Biogeosciences*, **5**, 561–583.

640 Ward, J.H. (1963) Hierarchical Grouping to Optimize an Objective Function. *Journal of the*
641 *American Statistical Association*, **58**, 236–244.

642 Ward, B.A., Friedrichs, M.A.M., Anderson, T.R. & Oschlies, A. (2010) Parameter optimization 643 techniques and the problem of under determination in marine biogeochemical models. *Journal* 644 of Marine Systems, 81, 34–43.

645 Wullschleger, S.D., Epstein, H.E., Box, E.O., Euskirchen, E.S., Goswami, S., Iversen, C.M., 646 Kattge, J., Norby, R.J., van Bodegom, P.M. & Xu, X. (2014) Plant functional types in Earth 647 system models: past experiences and future directions for application of dynamic vegetation 648 models in high-latitude ecosystems. *Annals of Botany*, 114, 1–16.

649 Zaehle, S., Sitch, S., Smith, B. & Hatterman, F. (2005) Effects of parameter uncertainties on 650 the modeling of terrestrial biosphere dynamics. *Global Biogeochemical Cycles*, 19.
<table><thead><tr><th>Parameter</th><th>Description</th><th>Equation used in ORCHIDEE-FM</th><th>Processes involved</th></tr></thead><tbody><tr><td>V<sub>cmax_opt</sub></td><td>Optimal maximum rate of carboxylation (µmol m<sup>-2</sup>s<sup>-1</sup>)</td><td>V<sub>c</sub> = <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">V<sub>cmax_opt</sub> × Ci</div><div style="padding: 0 5px;">Ci + K<sub>c</sub>(1 - <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">Oi</div><div style="padding: 0 5px;">K<sub>o</sub></div></div>)</div><br>V<sub>c</sub>= RUBISCO activity (µmolCO<sub>2</sub> m<sup>-2</sup>s<sup>-1</sup>)<br>Oi= intracellular oxygen concentration<br>K<sub>c</sub> and K<sub>o</sub>=Michaelis-Menten constants of RUBISCO activity</td><td>Impacts the gross primary production and substantially the stomatal conductance</td></tr><tr><td>T<sub>opt</sub></td><td>Optimal temperature of photosynthesis (°C)</td><td>V<sub>c</sub> = V<sub>cmax_opt</sub> × <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">(T<sub>air</sub> - Tmax)(T<sub>air</sub> - Tmin)</div><div style="padding: 0 5px;">(T<sub>air</sub> - Tmax)(T<sub>air</sub> - Tmin) - (T<sub>air</sub> - T<sub>opt</sub>)<sup>2</sup></div></div><br>T<sub>air</sub>= Air temperature (°C)<br>Tmin, Tmax= minimal and maximal temperature of photosynthesis (°C)</td><td>Impacts the temperature dependency of the photosynthesis, and daily and seasonal variations of the GPP.<br>It is the same equation used for V<sub>jmax</sub> (maximal rate of the RuBP regeneration)</td></tr><tr><td>g<sub>slope</sub></td><td>Slope of stomatal conductance in the Ball-Berry model (unitless)</td><td>g<sub>s</sub> = g<sub>min</sub> + gslope × A × <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">hr</div><div style="padding: 0 5px;">Ca</div></div><br>gs= stomatal conductance (mol m<sup>-2</sup> s<sup>-1</sup>)<br>A = net assimilation of the photosynthesis (mol m<sup>-2</sup> s<sup>-1</sup>)<br>hr = relative humidity (unitless)<br>Ca = boundary layer CO<sub>2</sub> concentration (µmol mol<sup>-1</sup>)<br>g<sub>min</sub> = empirical coefficients (mol m<sup>-2</sup> s<sup>-1</sup>)</td><td>Parameter which modifies the slope of the linear relationship between the stomatal conductance and the assimilation. It substantially impacts the estimation of the NEE and E.</td></tr><tr><td>SLA</td><td>Specific length area (mm<sup>2</sup> mgC<sup>-1</sup>)</td><td>LAI = SLA × Bm<sub>leaves</sub><br>LAI= Leaf area index (m<sup>2</sup> leaves m<sup>-2</sup>ground)<br>Bm_leaves= leaves biomass (mgC)</td><td>Parameter which determines the Leaf Area Index and impacts the GPP.</td></tr><tr><td>L<sub>age,crit</sub></td><td>Critical leaf age for leaf senescence (month)</td><td>V<sub>c</sub> = V<sub>cmax_opt</sub> × e<sub>rel</sub> (<div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">af</div><div style="padding: 0 5px;">Leaf<sub>agecrit</sub></div></div>)<br>af= leaf age<br>Leaf_agecrit= mean leaf lifespan<br>e<sub>rel</sub>=relative efficiency</td><td>Parameter which determines the mean leaf lifespan.<br>It impacts the photosynthetic activity and the senescence of leaves.</td></tr></tbody></table>
<table><caption>Table 1: List of parameters, equations and processes involved in ORCHIDEE-FM and selected for this study.</caption><tr><td>Vjmax/Vcmax</td><td>Maximum flux of electron / maximum rate of carboxylation ratio (unitless)</td><td>Vjmax = 2 × Vcmax<br/>Vjmax = maximal rate of the RuBP regeneration (µmol m<sup>-2</sup>s<sup>-1</sup>)</td><td>The parameter impacts the estimation of the GPP and substantially impacts the stomatal conductance.</td></tr><tr><td>a<sub>bc</sub></td><td>Coefficient of biomass-circumference allometry</td><td>M = a<sub>bc</sub> × D<sup>b<sub>bc</sub></sup></td><td rowspan="2">Determines the circumference of trees and impacts thinning events and the exportation of biomass</td></tr><tr><td>b<sub>bc</sub></td><td>Coefficient of biomass-circumference allometry</td><td>M = total aboveground tree dry biomass (kg)<br/>D = diameter at breast height (cm)</td></tr><tr><td>a<sub>st</sub></td><td>Coefficient of self thinning equation</td><td>dens<sub>max</sub> = α<sub>st</sub>/Dg<sup>β<sub>st</sub></sup></td><td rowspan="2">Determines the maximal density of trees in the stand. Impacts thinning events and consequently the biomass distribution among trees</td></tr><tr><td>β<sub>st</sub></td><td>Coefficient of self thinning equation</td><td>dens<sub>max</sub> = stand maximum density in ind ha<sup>-1</sup> (individuals per hectare)<br/>Dg = quadratic mean diameter (meters).</td></tr><tr><td>rdi_target</td><td>Targeted value or relative density index (unitless)</td><td>thinning<sub>if rdi</sub> = rdi<sub>target</sub> ± δrdi<br/>rdi = relative density index<br/>δrdi= tolerable variation of rdi around rdi_target</td><td>Determines the moment when the thinning occurs. Impacts the biomass and exportation</td></tr></table>

652

653

654

655

656
<table><thead><tr><th rowspan="2">Parameter</th><th colspan="5">Number of data from literature / TRY</th><th rowspan="2">Total</th></tr><tr><th>Abies alba</th><th>Picea abies</th><th>Pinus pinaster</th><th>Pinus sylvestris</th><th>Pseudotsuga menziesii</th></tr></thead><tbody><tr><td>V<sub>cmax25</sub></td><td>13 / 0</td><td>30 / 3</td><td>13 / 0</td><td>15 / 2</td><td>20 / 0</td><td>93</td></tr><tr><td>T<sub>opt</sub></td><td>2 / 0</td><td>2 / 0</td><td>2 / 0</td><td>2 / 0</td><td>2 / 0</td><td>10</td></tr><tr><td>g<sub>slope</sub></td><td>0 / 0</td><td>7 / 6</td><td>2 / 0</td><td>8 / 0</td><td>0 / 0</td><td>23</td></tr><tr><td>g<sub>smax</sub></td><td>10 / 0</td><td>20 / 1</td><td>21 / 0</td><td>18 / 1</td><td>15 / 0</td><td>86</td></tr><tr><td>SLA</td><td>3 / 23</td><td>13 / 10</td><td>20 / 3</td><td>24 / 224</td><td>10 / 6</td><td>336</td></tr><tr><td>L<sub>age,crit</sub></td><td>2 / 2</td><td>7 / 2</td><td>7 / 0</td><td>17 / 4</td><td>4 / 2</td><td>47</td></tr><tr><td>J<sub>max</sub>/V<sub>cmax</sub></td><td>19 / 0</td><td>25 / 2</td><td>7 / 0</td><td>11 / 1</td><td>5 / 0</td><td>70</td></tr><tr><td>b<sub>bc</sub></td><td>3 / 0</td><td>39 / 0</td><td>3 / 0</td><td>22 / 0</td><td>12 / 0</td><td>79</td></tr><tr><td>α<sub>st</sub></td><td>2 / 0</td><td>6 / 0</td><td>3 / 0</td><td>8 / 0</td><td>4 / 0</td><td>23</td></tr><tr><td>Rdi_target</td><td>1 / 0</td><td>1 / 0</td><td>1 / 0</td><td>1 / 0</td><td>1 / 0</td><td>5</td></tr></tbody></table>

Table 2: Number of observations retrieved for each parameter from the literature and the TRY database. Refers to Table 1 for the description of each parameter.
<table><caption>Table 3: Default (ORCHIDEE current value), optimized (Santaren et al., 2007; Kuppel et al., 2014) and observed (from the literature and the TRY database) values for each parameter used in this study. TNF correspond to the current Temperate Needleleaf evergreen Forest plant functional type. (-): no data available.</caption><thead><tr><th rowspan="2">Parameter</th><th rowspan="2">Description</th><th>Current value</th><th colspan="2">Santaren et al., 2007</th><th colspan="8">Median value from literature</th><th rowspan="2">Unit</th></tr><tr><th>TNF</th><th>TNF</th><th>TNF</th><th>A.alb</th><th>P.abi</th><th>P.men</th><th>P.syl</th><th>P.pin</th><th>PFT1</th><th>PFT2</th><th>TNF</th></tr></thead><tbody><tr><td>V<sub>cmax25</sub></td><td>Optimal maximum rate of carboxylation</td><td>35</td><td>41.7</td><td>31.94</td><td>46</td><td>46.77</td><td>38.35</td><td>106.61</td><td>90.3</td><td>45.5</td><td>100.2</td><td>61</td><td>µmol m<sup>-2</sup> s<sup>-1</sup></td></tr><tr><td>T<sub>opt</sub></td><td>Optimal temperature of photosynthesis</td><td>25<br>25</td><td>-</td><td>17.49<br>17.49</td><td>36.81<br>33.20</td><td>31.58<br>21.76</td><td>29.94<br>27.3</td><td>27.56<br>19.89</td><td>38.84<br>36.87</td><td>31.78<br>27.42</td><td>33.2<br>28.38</td><td>32.34<br>27.8</td><td>Vmax °C<br>Jmax °C</td></tr><tr><td>g<sub>slope</sub></td><td>Slope of stomatal conductance in the Ball-Berry model</td><td>9</td><td>8.1</td><td>8.841</td><td>-</td><td>5.94</td><td>-</td><td>5.59</td><td>10.05</td><td>-</td><td>-</td><td>-</td><td>unitless</td></tr><tr><td>g<sub>smax</sub></td><td>Maximal stomatal conductance</td><td>-</td><td>-</td><td>-</td><td>3.41</td><td>3.88</td><td>4.12</td><td>3.9</td><td>3</td><td>3.89</td><td>3.78</td><td>3.89</td><td>mm s<sup>-1</sup></td></tr><tr><td>SLA</td><td>Specific leaf area</td><td>4.63</td><td>-</td><td>10</td><td>5.51</td><td>5.4</td><td>7.08</td><td>4.73</td><td>3.1</td><td>5.6</td><td>4.73</td><td>4.85</td><td>mm<sup>2</sup> mg<sup>-1</sup></td></tr><tr><td>L<sub>age,crit</sub></td><td>Critical leaf age for leaf senescence</td><td>29.5</td><td>-</td><td>35.6</td><td>98.7</td><td>102.33</td><td>65.4</td><td>34</td><td>51</td><td>78</td><td>42.5</td><td>51</td><td>Months</td></tr><tr><td>J<sub>max</sub>/V<sub>cmax</sub></td><td>Maximum flux of electron / maximum rate of carboxylation ratio</td><td>2</td><td>-</td><td>-</td><td>2.2</td><td>2.2</td><td>3.01</td><td>2.1</td><td>2.5</td><td>2.2</td><td>2.1</td><td>2.2</td><td>unitless</td></tr><tr><td>a<sub>bc</sub></td><td>Coefficient of biomass-circumference allometry</td><td>7.03b<sub>bc</sub><sup>-4.76</sup></td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>kgDM</td></tr><tr><td>b<sub>bc</sub></td><td>Coefficient of biomass-circumference allometry</td><td>2.295</td><td>-</td><td>-</td><td>2.52</td><td>2.36</td><td>2.4035</td><td>2.2355</td><td>2.13</td><td>2.37</td><td>2.234</td><td>2.355</td><td>ln(kgDM)<br>ln(m)<sup>-1</sup></td></tr><tr><td>α<sub>st</sub></td><td>Coefficient of self thinning equation</td><td>12,2</td><td>-</td><td>-</td><td>12.548</td><td>12.7925</td><td>12.3165</td><td>12.48097</td><td>12.544</td><td>12.4995</td><td>12.544</td><td>12.5</td><td>ln(ind ha<sup>-1</sup>)</td></tr><tr><td>β<sub>st</sub></td><td>Coefficient of self thinning equation</td><td>1.6</td><td>-</td><td>-</td><td>1.7505</td><td>1.664</td><td>1.6255</td><td>1.7915</td><td>1.815</td><td>1.664</td><td>1.815</td><td>1.7295</td><td>ln(ind ha<sup>-1</sup>)<br>ln(m)<sup>-1</sup></td></tr><tr><td>Rdi_target</td><td>Targeted value or relative density index</td><td>0.75</td><td>-</td><td>-</td><td>0.59</td><td>0.54</td><td>0.63</td><td>0.44</td><td>0.59</td><td>0.56</td><td>0.64</td><td>0.56</td><td>unitless</td></tr></tbody></table>
<figure><img src="image_2.png" /><figcaption>Fig. 1: Representation of the ORCHIDEE-FM model and related processes. Round corner boxes correspond to ORCHIDEE modules (bold). The main simulated processes and the corresponding time step (Δt) are listed for each module. In the STOMATE module are listed the different carbon compartments, while the FM module illustrates the management processes simulated by ORCHIDEE. Arrows represent shared variables, processes and the different links between modules. Parameters used in this study are listed in parenthesis (Table 1).</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>Fig.2: Boxplot representation of parameters for each species (median, first and third quartile, first and last deciles, black open circles represent values outside 1.5 times the interquartile range above the upper quartile and below the lower quartile). Refers to Table 1 for the description of each parameter. Orange diamond represent mean values; blue dots represent values fitted on French national forest inventories data. Number of data available for each species is specified on the right of each boxplot. LIT/TRY correspond to the number of data from literature and from the</figcaption></figure>
677 TRY database. The red line represents the current ORCHIDEE value of the TNF PFT. Parameters with significant Kruskal-Wallis test : '*' p-

678 value < 0.05; '**' p-value < 0.01

679

680
<figure><img src="image_4.png" /><figcaption>Fig.2 (contd.)</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Fig.3: Clusters of species (*Abies alba*/*Picea abies*/*Pseudotsuga mensiezii*/*Pinus sylvestris*/*Pinus pinaster*) obtained with different parameter selections. a) with all parameters; b) with photosynthesis-parameters (P); c) with biomass and management-parameters (B); d) with a reduced set of P-parameters (rP) and e) with a reduced set of B-parameters (rB) and f) with significantly different parameters between species.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Fig.4: Distribution of the mean relative difference (“RD %”) between observed and simulated volume increment in Bellassen *et al.*, (2011) according to the new PFT classification for a) pixels where only one PFT is present and b) where both PFTs are present on the same pixel. Blue striped frequencies represent PFT1 (*Abies alba*/*Picea abies*/*Pseudotsuga mensiezii*), red frequencies PFT2 (*Pinus sylvestris*/*Pinus pinaster*). Vertical lines represent the median relative difference for each PFT. The asterisk indicates a significant difference between the two means at *P*<0.05</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Fig.5: Spearman's correlation coefficient between observed and simulated volume, volume increment and stand density for simulations: REF (standard), S1 (1 PFT), S2 (2 PFTs) and S3 (5 PFTs). Results are for 50y old stands over France</figcaption></figure>