Analysing Amazonian forest productivity using a new individual and trait-based model ( TFS v . 1 )

Repeated long-term censuses have revealed largescale spatial patterns in Amazon basin forest structure and dynamism, with some forests in the west of the basin having up to a twice as high rate of aboveground biomass production and tree recruitment as forests in the east. Possible causes for this variation could be the climatic and edaphic gradients across the basin and/or the spatial distribution of tree species composition. To help understand causes of this variation a new individual-based model of tropical forest growth, designed to take full advantage of the forest census data available from the Amazonian Forest Inventory Network (RAINFOR), has been developed. The model allows for withinstand variations in tree size distribution and key functional traits and between-stand differences in climate and soil physical and chemical properties. It runs at the stand level with four functional traits – leaf dry mass per area ( Ma), leaf nitrogen (NL) and phosphorus (P L) content and wood density (DW) varying from tree to tree – in a way that replicates the observed continua found within each stand. We first applied the model to validate canopy-level water fluxes at three eddy covariance flux measurement sites. For all three sites the canopy-level water fluxes were adequately simulated. We then applied the model at seven plots, where intensive measurements of carbon allocation are available. Treeby-tree multi-annual growth rates generally agreed well with observations for small trees, but with deviations identified for larger trees. At the stand level, simulations at 40 plots were used to explore the influence of climate and soil nutrient availability on the gross ( 5G) and net ( 5N) primary production rates as well as the carbon use efficiency (C U). Simulated5G,5N and CU were not associated with temperature. On the other hand, all three measures of stand level Published by Copernicus Publications on behalf of the European Geosciences Union. 1252 N. M. Fyllas et al.: Analysing Amazonian forest productivity using a new individual and trait-based model productivity were positively related to both mean annual precipitation and soil nutrient status. Sensitivity studies showed a clear importance of an accurate parameterisation of withinand between-stand trait variability on the fidelity of model predictions. For example, when functional tree diversity was not included in the model (i.e. with just a single plant functional type with mean basin-wide trait values) the predictive ability of the model was reduced. This was also the case when basin-wide (as opposed to site-specific) trait distributions were applied within each stand. We conclude that models of tropical forest carbon, energy and water cycling should strive to accurately represent observed variations in functionally important traits across the range of relevant scales.


Introduction
The Amazon basin, encompassing one of the planet's largest forest areas and hosting one quarter of the Earth's biodiversity, constitutes a large reservoir of living biomass (Malhi and Phillips, 2005).Amazon forests also have a substantial influence on regional and global climates (Shukla et al., 1990;Spracklen et al., 2012).These forests are, however, under strong human pressure through logging and forest-to-pasture conversion, and face at present a warming and more variable climate and changing atmospheric composition (Lewis et al., 2004;Gloor et al., 2013).Due to the enormous area of forest within the Amazon basin, these factors have the potential to modify global atmospheric greenhouse concentrations, regional and global climate, and the overall biodiversity of the planet (Cramer et al., 2004).
Traditionally, two approaches have been followed to understand current and future state of the Amazon forests.First, dynamic global vegetation models (DGVMs) have been used to simulate vegetation patterns and carbon fluxes across Amazonia (Moorcroft et al., 2001;Galbraith et al., 2010) with some predicting substantial carbon losses under scenarios of global change (White et al., 1999;Cox et al., 2004) but with others less so (Cramer et al., 2004), or even gains (Huntingford et al., 2013).A second approach to understand Amazonian forests dynamics is through the analysis of longterm field observations of patterns of tree growth and mortality as they relate to climatic and edaphic variations across the basin (e.g.Phillips et al., 2004;Quesada et al., 2012).
Analyses of Amazon forest inventory data, and particularly those of the Amazon Forest Inventory Network (RAIN-FOR) (Malhi et al., 2002), have revealed large-scale temporal trends in biomass and species composition as well as intriguing spatial patterns in many stand properties (Phillips et al., 1998;Baker et al., 2004;Phillips et al., 2009).Specifically, there is systematic spatial variation in species composition, biomass, growth and turnover rates, with western forests exhibiting higher wood productivity, faster turnover time and lower stand wood density compared to eastern forests (Baker et al., 2004;Malhi et al., 2006).This macroecological variation may possibly be explained by the basin-wide observed climate and soil nutrient availability gradients (Ter Steege et al., 2006;Quesada et al., 2012).The climatic gradient comprises a southeast to northwest increase in annual precipitation and decrease in dry season length (Sombroek, 2001), with aboveground wood productivity positively related to precipitation (Malhi and Wright, 2004).On the other hand, a soil age/nutritional axis spans from the northeastern part of the basin to southwestern Amazonia, with generally younger and richer soils in the west and highly weathered nutrient poor soils in the east (Sombroek, 2000;Quesada et al., 2011), although at regional and local scales the patterns are often more complicated than this macro-gradient might imply (Higgins et al., 2011).Soil physical properties (such as rooting depth, drainage and water holding capacity and soil structure) are similarly related to soil age and parental material (Quesada et al., 2010).Poor physical (for example soil depth) conditions (less weathered soils) are often associated with higher soil nutrient availability (Walker and Syers, 1976;Vitousek and Farrington, 1997), leading to increased nutrient concentrations at the leaf level (Fyllas et al., 2009) and thus a potential for higher photosynthetic rates (Reich et al., 1994;Raaimakers et al., 1995).In addition, increased disturbanceassociated mortality rates in soils of poor physical properties lend towards more dynamic stands where faster growing species dominate (Chao et al., 2009;Quesada et al., 2012).This positive feedback mechanism could explain the higher aboveground productivity and turnover rates observed for western forests (Quesada et al., 2012).
The simplistic ways by which plant functional diversity is currently reflected in DGVMs is an important shortcoming in predicting ecosystem response to environmental gradients and their vulnerability to global change (Lavorel et al., 2007).Some of the widely applied DGVMs represent Amazonian plant diversity with only few plant functional types (PFTs), for example the LPJ model uses only two tropicaloriented PFTs (Sitch et al., 2003) and the JULES model only one (Clark et al., 2011).The mean values of key model parameters like photosynthetic capacity, wood density and leaf turnover times are selected to describe an a priori PFT definition (Fyllas et al., 2012).This means that many processes are controlled by a set of fixed parameters that describe viable plant strategies within very limited boundaries.Such PFT implementation has important drawbacks.It is usually based on the average value of a plant trait recorded from different field studies and different species.But recent studies have shown that key traits present a wide variation, dependent upon species identity and site growing conditions (Sultan, 2000;Fyllas et al., 2009;Baraloto et al., 2010a).Thus any given species has the potential to exhibit site-dependent shifts in its trait value; this being in addition to the inter-specific trait variability expected at any given site.Ignoring this plasticity could potentially bias modelling through an underestimation of the PFT's resilience by projecting dramatic but artificial switches in vegetation state caused by the limited and discrete (step-wise) nature of PFT descriptions.
Such unaccounted variability could be particularly important when modelling Amazonian forest dynamics, where environmental heterogeneity and plant functional diversity comprise key components of the ecosystem (Townsend et al., 2008).For example, the variation in leaf mass per area (M a ) recorded within Amazon forests covers an approximately similar range to the one identified in global data sets, ranging from 30 to 300 g m −2 (Fyllas et al., 2009).Similarly, there are large contrasts in soil physical and chemical conditions (Quesada et al., 2010).These important ecosystem flux drivers have now been better quantified, with Amazonwide climate (Malhi and Wright, 2004), soil (Quesada et al., 2011) and functional trait data sets also having been obtained (Baker et al., 2009;Fyllas et al., 2009;Patiño et al., 2009;Patiño et al., 2012).This is in addition to continually expanding long-term forest inventory data in which tree growth, mortality and species composition data are regularly being recorded (Keeling et al., 2008;Chao et al., 2009).
We here introduce a vegetation dynamics model developed as a tool to better analyse these observed Amazonian large-scale productivity patterns.This is achieved through specific incorporations of observed environmental and the biotic variations into the model formulation.Specifically we focus (a) on the architectural variability, expressed through the size-class distribution of a stand, and (b) on the functional variability, expressed through simulated distributions of four important functional traits which are allowed to vary from tree to tree within individual plots.Following a continuum approach, we replace the use of a discrete number of PFTs, with distributions of a functional trait "quartet", the withinstand distributions of which also vary from plot to plot in accordance with observation.
Two axes of functional variation/strategy are represented in the model: the leaf economic and the tree architecture spectra.The four functional traits include leaf mass per area (M a ), leaf nitrogen and phosphorous dry mass concentration (N Lm and P Lm respectively) and wood density (D W ). The first three traits express one component of the leaf economic spectrum (Reich et al., 1997;Wright et al., 2004), i.e. a global photosynthetic tissue trade-off between inexpensive, short-lived and fast payback leaves vs. costly, longlived and slow payback leaves, although we emphasise that other factors such as leaf cation concentrations may be important in this respect (Fyllas et al., 2012;Patiño et al., 2012).Low M a and high nutrient content leaves are associated with comparably short longevity and usually have high (massbased) gas exchange rates (Reich et al., 1994;Raaimakers et al., 1995).Lately the role of P Lm has been highlighted, as it expresses alternative limitations of the photosynthetic efficiency of tropical tree species (Domingues et al., 2010).The fourth trait, D W , is used to represent a tree architectural axis with denser wood species supporting an overall higher aboveground biomass and thus having a higher maintenance respiration (Chave et al., 2005;Mori et al., 2010, although see Larjavaara andMuller-Landau, 2012).These two dimensions capture essentially a growth vs. survival trade-off.There is mixed evidence for a coordination between leaf and stem traits, i.e. a correlation between slow return related leaf traits and denser wood (Chave et al., 2009), with Baraloto et al. (2010b) suggesting that these two axes are independent, but with Patiño et al. (2012) showing some important correlations with foliar traits such as P Lm .For the purpose of this study we consider leaf and stem dimensions as independent axes of tree functional variation, with no predefined interrelationship between the representative traits.However, the observed among-stand variability of these four characters is used to express how growing conditions control plant processes, while the within-stand trait variation represents a range of ecological strategies found under the same growing conditions.
The model is initialised with site-specific tree diameter and functional traits data, and forced with daily climate data.We first test the ability of the model to estimate stand-level water fluxes at three eddy-flux tower sites.For a subset of seven RAINFOR plots where site-specific carbon allocation coefficients are known, a tree-level test of stem growth rates is applied.We further validate the ability of the model to simulate the spatial patterns of aboveground biomass productivity at 40 RAINFOR plots, and subsequently explore the variation of Gross Primary Productivity ( G ), Net Primary Productivity ( N ) and Carbon Use Efficiency (C U ) along established Amazonian climatic and edaphic gradients.

Model description
"Traits-based Forest Simulator" (TFS) is an individual-based forest model, i.e. it simulates water and carbon fluxes for each tree in a stand.In the current version of the model, stand structure is prescribed in terms of the number of trees and their diameter at breast height (d).This is thus a "snapshot" version of the model, which does not take into account tree recruitment and mortality.In this version of TFS, each individual is fully described through d, with allometric equations used to estimate other attributes of interest like tree height (H ), crown area (C a ), total leaf area (L a ) and treelevel leaf area index (L).Whole tree biomass is then partitioned to leaf (B L ), stem (B S ), coarse root (B CR ) and fine root (B FR ) biomass using established allometric equations.Allocation of assimilated carbon to different plant components is static, i.e. it does not change with size or resource availability, but rather implements field-derived allocation coefficients (Aragão et al., 2009).The general architecture of the model is presented in Fig. 1.Tree functional diversity is expressed through four traits (M a , N Lm , P Lm , D W ), which are randomly assigned to each tree: these pseudo-data being generated from local observations using a random vector generation algorithm.Leaf photosynthesis is calculated using a modified version of the Farquhar biochemical model (Farquhar et al., 1980), that incorporates leaf chemical and soil moisture effects.The maximum photosynthetic rate is regulated by N L or P L through the co-limitation model of Domingues et al. (2010).In contrast to most ecosystem fluxes models, where photosynthetic rates are directly regulated by water availability (Scheiter and Higgins, 2009;Clark et al., 2011), we couple water "stress" to reduction of canopy conductance by estimating a daily fractional available soil water content for each tree in the stand.Carbon fluxes are simulated on an hourly and water fluxes on a daily time step.
Light competition is based on the assumption of a perfect canopy tessellation.The flat-top version of the perfect plasticity model (Purves et al., 2007) has been used in the current version of TFS to characterise canopy and sub-canopy trees, by assuming that all of a tree's foliage is found at the top of its stem (S1, Canopy Architecture and Radiation Environment).A canopy height Z* is estimated for each forest stand, defining canopy and sub-canopy trees.By summing up the crown area (C a ) of all trees in the stand, Z* is estimated as the height of the last tree that enters to the sum before the cumulative crown area is equal to the plot area.Canopy trees are absorbing a mean daily amount of short-wave solar radiation equal to the sum of mean beam, diffuse and scattered daily radiation in correspondence to the sun-shade model of de Pury and Farquhar (1997).The direct and diffuse fraction of solar radiation is estimated using the Spitters et al. (1986) approximation.The functional configuration of a tree (i.e. the values of the quartet of traits) does not affect its light competitive status, as tree height and crown area are not directly associated with any of the four traits.Future versions of the model will incorporate such effects.
Soil water balance is approximated through a simple bucket model, with soil water content affecting leaf conductance and thus photosynthetic rates.Competition for soil water is approximated through a size hierarchy, i.e. bigger trees, with a more extensive root system are assumed to have access to deeper water (S1, Water Balance Algorithm).By assuming that a tree with a higher leaf biomass (B L ) requires a higher fine root biomass (B FR ), we indirectly implement a M a effect on water competition (S1, Definition, Allometry and Stoichiometry of Individual Trees in TFS).In particular, between two trees of the same size, the higher M a tree will be more competitive in terms of acquiring soil water.
TFS is coded in Java and it is fully described in S1.The main effects of including functional diversity are realised through trait-driven effects on photosynthesis and respiration (Reich et al., 2008(Reich et al., , 2009)).Model components that are linked with any of the four base traits are described in the following paragraphs.All statistical analyses and graphs were made with R (R Development Core Team, 2013).

Within-stand functional diversity
As noted above, TFS employs neither species nor PFT descriptions, but rather a different discrete combination of each the four key functional traits M a , N Lm , P Lm and D W is assigned to each individual tree along with a diameter-based allometry.To achieve this, the four functional characters assigned are generated using a procedure based on the actual values recorded within each plot.This is achieved using a random vector generation algorithm (Taylor and Thompson, 1986).This algorithm, appropriate for generating nonrepeated pseudo-observations from a relatively small sample of observations, was originally developed to provide for a realistic probabilistic representation of shrapnel projectile distributions in military battlefield simulations in the face of only a limited amount of available data (due to the cost and difficulty of undertaking the appropriate experiments).This "ballistic method" is notable in that it was specifically designed to short-circuit the usual step of multivariate density in the generation a pseudorandom population with approximately the same moments as the original sample.The ballistic method is readily programmable as follows (with the underlying rationale as discussed in Taylor andThompson, 1986 andThompson, 1989) and with the following description based on Visual Numerics (2014).
First take a vector X with n multivariate observations (x 1 , . . ., x N ).To generate a pseudo data set from x, one observation (x j ) is first chosen at random and its nearest m neighbours, x j 1 , x j 2 , x j m are then determined and with the mean x j of those nearest neighbours subsequently calculated.Next, a random sample u 1 , u 2 , . . ., u m is generated from a uniform distribution with lower bound 1 m − 3(m−1) m 2 , and upper bound 1 m + 3(m−1) m 2 .The random variate is then estimated as m l=1 u 1 (x j l − x j ) + x j and the process then repeated as required.Somewhat subjective here is the selection of the appropriate value of the number of nearest neighbours (m) although the nature of the simulations is not strongly dependent upon that value (Taylor and Thompson, 1986).Thus, following their recommendation and as in the Visual Numerics (2014) default, we have taken here m = 5.
In our case, applying this procedure resulted in a coordinated trait quartet for each tree in a stand being generated on the basis on the smaller observational trait quartets sampled from trees in the same stand (Baker et al., 2009;Fyllas et al., 2009;Patiño et al., 2012) and without any assumptions having to be made about their underlying statistical distributions.Thus no single functional trait "average stand" value is used (or even required).Further, between-stand differences in the trait distributions and their covariances are also intrinsically taken into account.This is because each stand is characterised by its own multivariate trait sample and size distribution.More fertile plots have an overall lower M a and higher N Lm and P Lm compared to infertile plots (Fyllas et al., 2009), with this being reflected in the photosynthetic capacity of individual trees, as described in the next paragraph.

Photosynthesis
A tree-level leaf area index (L), estimated as the ratio of L a to C a , is used to compute the energy, carbon and water fluxes for each tree in a stand.The net photosynthetic rate (µmol m −2 s −1 ) is given by with C α the atmospheric CO 2 mixing ratio (µmol mol −1 ), C c the CO 2 mixing ratio inside the chloroplast and g S the CO 2 stomatal conductance (mol m −2 s −1 ) calculated from Medlyn et al. ( 2011) and modulated by a soil moisture term.The leaf-level photosynthetic rate A N is scaled up to the tree-level by multiplying with the C a of the tree.
The co-limitation equation suggested by Domingues et al. (2010), whereby the leaf level photosynthetic capacity (area basis) is potentially limited by either nitrogen or phosphorus, is used in TFS to estimate the leaf maximum carboxylation and electron transport rates: ).The canopy-level photosynthetic capacity V Cmax (µmol m −2 s −1 ) is estimated using the tree-level leaf area index L, taking into account within-canopy gradients in light and photosynthetic capacity based on Lloyd et al. (2010).Nutrient optimisation is approximated using equations in Lloyd et al. (2010), with M a also dependent on the height of each tree (H i ) and the mean canopy height (H S ): with a H an empirical coefficient.

Respiration
Tree respiration includes a growth and a maintenance component, both computed daily.Growth respiration is considered as a constant fraction (0.25) of daily photosynthesis (Cannell and Thornley, 2000).Three different maintenance respiration formulations are allowed in TFS (S1, Respiration), but in this study we use the one described below.Leaf maintenance respiration R mL is estimated as a fraction of V Cmax (Scheiter and Higgins, 2009): Stem maintenance respiration is estimated from the sapwood volume (V S ) of a tree: with δ = 39.6 (µmol m −3 s −1 ) as reported in Ryan et al. (1994) for tropical trees.Sapwood volume is estimated by inversion of the pipe model and assuming that the ratio of leaf area to sapwood area ( LS ) increases with the height and the wood density for tropical trees (following Calvo-Alvarado et al., 2008;Meinzer et al., 2008): with λ 1 = 0.066 m 2 cm −2 , λ 2 = 0.017 m cm −2 , δ 1 = −0.18m 2 cm −2 and δ 2 = 1.6 cm 3 g −1 .
Sapwood area (m 2 ) and volume (m 3 ) are then calculated from with L a the total leaf area of the tree (m 2 ) and with C D the crown depth (m).
Coarse-root maintenance respiration R mCR is estimated as in Scheiter and Higgins (2009): where CN is the root C : N ratio estimated on the basis of the simulated N R assuming a dry weight carbon fraction of 0.5.Fine-root maintenance respiration R mFR is assumed to be equal to leaf respiration.
All respiratory components are corrected with the temperature dependence function of Tjoelker et al. (2001).The total maintenance respiration R m is then (11)

Stomatal conductance
Initially, a maximum (no water stress) stomatal conductance, g s,max is calculated following Medlyn et al. (2011Medlyn et al. ( , 2012)): with g 0 (mol m −2 s −1 ) the minimum stomatal conductance, g 1 (−) an empirical coefficient that represents the water use efficiency of the plant, and D C the leaf-to-atmosphere vapour pressure difference.Values of g 0 and g 1 that lead to the best model performance were different between sites, as indicated by the model calibration procedure.For the basin-wide simulations constant values of g 0 = 0.020 (mol m −2 s −1 ) and g 1 = 5.0 (−) were used, close to the estimates of Domingues et al. (2014).In future versions of the model, we anticipate that g 0 and g 1 will be related to other functional traits.The maximum stomatal conductance is subsequently reduced to the actual g S by multiplying the second term of Eq. ( 8) with a water stress coefficient.
In contrast to most ecosystem flux models, where photosynthetic rates are directly regulated by water availability (Scheiter and Higgins, 2009;Clark et al., 2011), we couple soil water deficit to canopy conductance by estimating a daily fractional available soil water content ϑ i for each i tree in the stand (S1, Water Balance and Soil Water Stress).This term is then used to estimate the water stress γ i that has a direct effect on stomatal conductance, as also described in Keenan et al. (2010).

Study sites and simulations setup
Three sets of site data were used to explore the behaviour of the model.These include a set of three eddy flux measurement (EFM) sites, seven plots with intensive carbon balance and allocation measurements (IMs), and 40 permanent measurement (PM) plots.

Eddy flux sites
Daily climate and energy flux data from three EFM sites (Caxiuanã [1.72  W]) were used to assess the ability of the model to estimate canopy-level water fluxes.Data were obtained from the Large Scale Biosphere-Atmosphere Experiment in Amazonia (LBA) project (http://daac.ornl.gov/LBA/lba.shtml).In particular mean daily climate parameters including incoming radiation, temperature, precipitation, relative humidity and wind speed were used to force the model.Latent heat flux (λE in W m −2 ) was used to estimate a daily mean canopy conductance defined as G C = λE D C .The EFM data cover a period from 2001 to 2008 for Caxiuanã, from 2000 to 2005 for Manaus and from 2002 to 2004 for Tapajós.G C was only estimated for days with a complete diurnal record of λE.At each one of the EFM sites the mean daily G C (mol m −2 s −1 ) was compared between observations and simulations.The model was initialised with size-class distribution and functional traits data from RAINFOR permanent plots located near the eddy flux towers.Specifically, CAX-06 inventory data were used for Caxiuanã, BNT-04 for Manaus, and TAP-55 for Tapajós.We note that the EFM sites are mainly found at the eastern part of Amazonia (Fig. 2) growing on low nutrient status soils.
The model was initially calibrated to the site-specific values for g 0 and g 1 of Eq. ( 8) that gave the best performance.A standardised major axis (SMA) regression, forced through zero, was used to verify the ability of the model to simulate G C , with a regression slope close to one indicating a good model performance.

Intensive measurement (IM) sites
The ability of the model to realistically simulate carbon fluxes at the tree level is evaluated using data from the seven intensive measurement plots (Aragão et al., 2009;Malhi et al., 2009).These sites are amongst the intensively surveyed plots within the RAINFOR network (Fig. 2), where measurements of all major components of the C cycle are recorded (Malhi et al., 2009).At these plots, a detailed assessment of the carbon stocks is applied, and N allocation coefficients to different plant components are estimated (Aragão et al., 2009;Malhi et al., 2011;Doughty et al., 2014).These sitespecific coefficients are used to calculate the amount of simulated N that is allocated to stems N,s (kg C yr −1 ).
The IM sites of interest include two plots at Agua Pudre in Colombia (AGP-01 and AGP-02), one (ALP-30) at Allpahuayo/Peru, one (BNT-04) at Manaus/Brazil, one in Caxiuanã/Brazil (CAX-06), one in Tambopata/Peru (TAM-05) and one in (TAP-55) Tapajós/Brazil.Based on data from Quesada et al. (2011), AGP-01, AGP-02, TAM-05 can be considered to be located on fertile soils, with the other four plots on infertile ones.Available soil depth data (Quesada et al., 2011) and functional traits data (Fyllas et al., 2009) were used for site-specific simulations.For all seven sites we estimated the observed average multi-annual growth rate (2000)(2001)(2002)(2003)(2004)(2005)(2006) of each tree from forest census data, in order to compare it with the simulated N,s .
The daily climate was extracted from the Princeton Global Meteorological Forcing Data Set (Sheffield et al., 2006).These simulations are used to validate the ability of the model to accurately estimate tree-level stem growth, under a given stand structure, a given climatic and soil profile and functional traits configuration of the established trees.Average observed stem growth rate (per 10 cm d bins), expressed in carbon units (i.e.kg C yr −1 ), is compared with simulated N,s using the York method of best straight line, which holds when both x and y observations are subject to correlated errors that vary from point to point (York et al., 2004).

Permanent measurement (PM) sites
Inventory data from 40 RAINFOR permanent measurement plots (Fig. 2), including tree diameter and multi-annual growth for all trees greater than 10 cm curated/managed in ForestPlots.net(Lopez-Gonzalez et al., 2011) Climate data for the same period were used here with the first year again used as a spin-up period (Sheffield et al., 2006).For those 40 PM plots, sample distributions of the traits quartet are available (Fyllas et al., 2009) as well as a description of soil chemical and physical properties (Quesada et al., 2011).
At the PM sites the simulated stand-level aboveground N was compared with observed rates of aboveground growth B ABG [kg C m −2 yr −1 ] for trees that survived during the 2000-2006 time period using a SMA regression.A second step was to explore the way that G , N and C U vary across an Amazon climatic and soil nutrient availability gradient (Quesada et al., 2010).The site scores of a principal components analysis (PCA) on the soil properties of the 40 PM plots (see Fyllas et al., 2009) are used to categorise plots along a nutrient availability gradient ( 1 ), while the key climatic variables used were the annual mean temperature T a and annual total precipitation P a .A Kendall correlation coefficient (τ ) was used to identify potential relationships of G , N and C U with T a , P a and 1 , as in most cases non-linear associations were observed.

Randomisation exercise
In order to explore (a) the importance of including trait variability and thus functional diversity in our simulations and (b) the importance of including constraints that are known to control the large-scale patterns of Amazonian forest dynamics, we conducted a randomisation exercise with the model being run under four alternative setups at the 40 permanent RAINFOR plots.The first setup denoted as var-tr is the variable-trait simulation with trait initialisation based on the observed stand-level trait distribution as described in the previous paragraphs (default setup).The second setup, denoted as fix-tr, is a fixed-trait simulation with all trees having the same (data set mean) values for each trait: this thus representing a single PFT case.The third setup (rand-tr) is a variable-trait simulation with trait initialisation based on random values of the trait quartet as recorded in any individual along the 40 permanent plots.This setup thus ignores any potential patterns of functional trait biogeography, i.e. traits are not related to the environmental or edaphic conditions under which a tree is growing.The fourth setup (rand-tr-N) is a variable-trait simulation in which the photosynthetic capacity of an individual is only defined by its leaf N content and thus the NP co-limitation constraint is removed.These alternative setups were compared by considering both the slope and the R 2 of SMA regressions between the predicted and the observed N,S .

Canopy conductance simulations at the EFM sites
Values of best model performance for g 0 and g 1 were different between sites, with g 0 = 0.035 (mol m −2 s −1 ) and g 1 = 7.5 at Caxiuanã, g 0 = 0.035 and g 1 = 7.0 at Manaus with g 0 = 0.01 and g 1 = 2.5 these being somewhat lower than the estimates of Domingues et al. (2013) at Tapajós.Simulated G C was underestimated for Caxiuana (α = 0.85 ± 0.05) and Manaus (α = 0.90 ± 0.02), with the model overestimating G C in Tapajós (α = 1.28 ± 0.04), but exhibiting an overall adequate performance (Fig. 3).For simulations at the IM and the PM sites, constant values of g 0 = 0.02 (mol m −2 s −1 ) and g 1 = 5(−) were used, which are found within the range of values in the EFM sites and reported estimates (Medlyn et al., 2012;Domingues et al., 2013).

Stem growth rate simulations at the IM sites
The mean simulated stem growth rate N,s of each tree in the seven IM plots was compared with the observed aboveground biomass gains ( B ABG ) for the 2000-2006 period.An accurate simulation of N,s can be seen for small size classes, but with greater differences between the observed and the simulated multi-annual growth found for bigger trees (Fig. 4).At infertile ALP-30, the estimate slope of the York model indicated an overestimation of aboveground production (α = 1.18±0.06),driven mainly by an overestimation of the mid-size classes.At BNT-04 the model underestimated the overall growth (α = 0.82 ± 0.03).Aboveground growth was overestimated in CAX-06 (1.11±0.07).At TAP-55 (α = 1.44 ± 0.15) the model underestimated aboveground production (0.90 ± 0.06).At fertile AGP-01 (α = 1.36 ± 0.08) and AGP-02 (α = 1.25±0.05)an overestimation of aboveground productivity was observed although with simulations of most size classes falling within the observed ranges.At TAM-05 (α = 0.79±0.07)though, the simulated aboveground growth was underestimated with the overall slope driven by divergences in smaller size classes.The range and distribution of N allocation to stem growth is adequately captured by TFS as summarised in Fig. S2.1 in the Supplement.

G , N and CU simulations at the PM sites
Simulated stand-level aboveground net primary productivity N,A was positively associated with observed changes in aboveground biomass of trees that survived in the PM plots over the 2000-2006 period B ABG , with an R 2 = 0.42, suggesting an adequate model behaviour (Fig. 5).A summary of simulated stand-level G , N and C U relationship to key environmental drivers is given in Table 1 (see also Fig. S2.2 in the Supplement).G and N and C U were not associated with temperature.However, all three measures of stand-level productivity were positively related to annual precipitation and soil nutrient availability.

Randomisation exercise simulations
Results from the randomisation exercise (Fig. 6) found the fully constrained default setup (var-tr) to have the best predictive performance (R 2 = 0.42 with a SMA slope a = 0.92).This is as compared to the fixed-trait simulation (fix-tr) single PFT parameterisation with a decreased predictive ability TFS (R 2 = 0.29, a = 0.82) and an overall higher mean predicted aboveground productivity.Not accounting for the site-specific distribution of the trait quartet, i.e. bypassing potential biogeographic patterns of functional diversity and/or environmental trait interactions (rand-tr), also reduced the predictive ability of the model (R 2 = 0.29, a = 0.74).Finally, the random trait no-NP co-limitation setup (rand-tr-N) similarly led to an inferior model performance (R 2 = 0.33, a = 0.88) and with the highest mean simulated aboveground productivity.

Discussion
We report here on the core components of an individualbased model that has been developed in order to help better understand the patterns revealed by recent integrated measurements of climate, soils, functional diversity and stand dynamics for a wide range of forests across the Amazon basin.Table 1.Kendall correlation coefficients (τ ) and associated significance levels (p) between simulated gross primary productivity ( G ), net primary productivity ( N ), carbon use efficiency (C U ) and key environmental factors.Significant associations are indicated with bold.In its current setup the model does not explicitly simulate regeneration and mortality dynamics but rather uses the observed size distribution of trees at the study sites, thus taking into account stand structure and functional trait variability as observed along the main climatic and edaphic axes of the Amazon basin.With the current setup we were able to reproduce the tree-and stand-level N patterns found across Amazonia and to explore for potential environmental controls over stand-level G , N and C U .

Scientific outcomes
Our simulations found no association of stand-level gross primary productivity ( G ) with temperature, probably due to the relatively small range of variation of temperature across our plots.G decreased until an annual temperature of approximately 26 • C but remained relative constant above this point (Table 1, Fig. S2.2 in the Supplement).However, our simulations suggest that a strong association of G with the annual precipitation and soil nutrient availability of the plots.
G was positively related to annual precipitation over the entire range observed in the 40 PM plots.The association of G with the nutrient availability axis is in agreement with fertilisation experiments showing an increase with nutrient supply (Giardina et al., 2003).In our basin-wide examination of G the soil nutrient availability and stand structure gradients are not, however, independent (Quesada et al., 2012), as in the RAINFOR network permanent plots it has been observed that bigger/older trees are more abundant on eastern infertile forests, where soil physical conditions can support a bigger tree size (Baker et al., 2009) with a lower risk of trees being uprooted (Chao et al., 2009).Bigger trees generally support a greater foliage area and thus could significantly contribute to the overall carbon assimilation of the stand.However, bigger trees on infertile plots are generally characterised by lower leaf nutrient concentrations (Fyllas et al., 2009) and thus slower assimilation rates (Reich et al., 1994;Domingues et al., 2010).On the other hand a higher abundance of smaller trees with higher gas exchange rates is observed on more dynamic, fertile plots.Ultimately this indicates that stand structure should be specifically taken into account when simulating G in tropical forests, and thus individual-based models could significantly contribute towards a deeper understanding of the functioning and sensitivity of these ecosystems.
In our simulations stand-level net primary productivity ( N ) showed no significant association to annual temperature but increased with soil nutrient availability and annual precipitation (Table 1, Fig. S2.2 in the Supplement).Our N simulations are in agreement with field observations of increasing aboveground wood productivity with precipitation (Quesada et al., 2012).Based on TFS parameterisation, photosynthetic rates are expected to be higher at a greater soil nutrient availability due to associated higher leaf N and P concentration (Fyllas et al., 2009;Domingues et al., 2010).Using a similar parameterisation for a "sun and shade" big leaf model, Mercado et al. (2011) found an increase in net canopy assimilation rate with leaf P content in agreement with our positive association between N and soil nutrient availability.Their simulated G accounted for approximately 0.30 of the observed wood productivity in 33 study plots, and thus the R 2 = 0.42 between simulated N and aboveground growth found here suggests a marginally improved model behaviour.It should be noted that our definition of soil nutrient availability ( 1 ), based on the PCA analysis in Quesada et al. (2010), directly relates to soil P content.As shown first in the analysis of Quesada et al. (2012), where data from almost 60 plots were considered, aboveground N is positively related to soil P content in lowland tropical forest.The increased N in fertile environments (apart from the higher G ) seems to be enhanced by the greater abundance of small trees there.As tree size increases maintenance respiration likely "consumes" an increasing proportion of assimilated carbon, and thus at large size classes the proportion of trees which have enough carbon to allocate to growth decreases (Givnish, 1988;Cavaleri et al., 2008).This is in line with the negative relationship between coarse wood production and maximum height documented for some Amazonian trees (Baker et al., 2009).
In our simulations carbon use efficiency (C U ) ranged from 0.43 to 0.54.Recent research suggests that the C U is not as constant as had been previously suggested (DeLucia et al., 2007;Zhang et al., 2009).For example the meta-analysis of DeLucia et al. (2007) found that C U varies from 0.23 to 0.83 in different forest types.Our average estimate of C U = 0.51 is, however, above the range reported in Malhi (2012).Zhang et al. (2009) identified a negative trend of the N / G ratio with temperature at the range of 20-30 • C, as also simulated here especially above 26 • C (Fig. S2.2 in the Supplement).Simulated C U increased with soil nutrient availability, being marginally lower at infertile (0.48) compared to fertile (0.50) plots.This is attributable to smaller size class trees (with lower relative respiratory costs) constituting a greater proportion of the total stand biomass on higher nutrient status soils.One factor relating to soil nutrient availability but not included in the current version is an implicit consideration of the respiratory costs of plant nutrient uptake (Lambers et al., 2008) either directly, or through other processes such as organic acid exudation (Jones et al., 2009) or the symbiotic associations (Duponnois et al., 2012).One would expect these costs to be proportionally higher for a stand of low nutrient status, especially with regard to P (Quesada et al., 2012).

Practical implications
The modelling of tropical forest carbon fluxes and stand dynamics has traditionally involved approaches aimed at a balance between simplicity, computational economy, and complexity.On one hand, the enormous biological and biogeochemical heterogeneity of tropical forests (Townsend et al., 2008) places special importance on how modellers prioritise both the amount and the detail of processes that should be included to capture the main controls and feedbacks.On the other hand, the finding that Amazonia is dominated by just 227 tree species (Ter Steege et al., 2013) implies that most biogeochemical cycling in the world's largest tropical forest is performed by a tiny sliver of its diversity.At one end of the complexity spectrum are individual-based models which are able to properly simulate population dynamics and thus lags due to demography.Individual-based models of tropical forests have traditionally focused on realistically representing the light environment (TROLL - Chave, 1999) or grouping tree species on the basis of their different responses to environmental resources as suggested by field observations (FORMIND - Köhler andHuth, 1998, LPJ-GUESS -Hély et al., 2006).At the other end of the complexity spectrum are DGVMs which simulate population dynamics more simplistically (but see Moorcroft et al., 2001;Scheiter and Higgins, 2009).Using a DGVM model, Verheijen et al. (2013) allowed for within-PFT climate-driven trait variation to occur and achieved an improvement of the predicted vegetative biomass and PFT distribution patterns.A similar rationale was followed in Wang et al. (2012) where it was shown that the inclusion of multi-trait covariance in DGVM can be used to constrain model parameters and reduce uncertainties in simulated ecosystem productivity.Fisher et al. (2010) applied the individual-based Ecosystem Demography model (Moorcroft et al., 2001), and showed that by varying traits related to demographic processes, forest and biomass dynamics exhibited a wide range of responses to climate forcing.
Most of the above approaches have used discrete PFTs to represent tree species and functional diversity.These studies suggest that by allowing for within-PFT trait variability a more plastic and realistic response to the relevant environmental drivers is observed.In contrast to the above, TFS replaces the use of PFTs with trait distributions, following a different model philosophy and architecture using the concept of multidimensional trait continua.In particular, considering functional diversity to be expressed by a multidimensional trait space, the use of PFTs selects a number of clusters where the central vector defines the average trait values of each PFT (Fyllas et al., 2012).Recent studies (Verheijen et al., 2013;Wang et al., 2012) allow for the average trait values to be shifted based on empirical climatic and/or trait inter-correlation functions.In contrast, the use of trait continua does not cluster the multidimensional trait space but rather allows any realistic trait combination (as suggested by the limited sampling of the actual population) to be simulated.Successful trait combinations under given environmental conditions may then be expected to emerge as a byproduct of model dynamics (Higgins et al., 2014).A similar to TFS representation of functional diversity has been implemented in the DGVM model (Scheiter and Higgins, 2009;Scheiter et al., 2013) where the importance of including trait variability in simulations of vegetation dynamics has also been highlighted.In TFS, variable-trait (R 2 = 0.42) simulations led to a better model performance compared to fixedtrait (R 2 = 0.29) simulations (Fig. 6).Thus including functional diversity in simulations of vegetation dynamics is expected not only to suggest less vulnerable communities under changing climatic conditions (Fauset et al., 2012;Scheiter et al., 2013) but also, it seems, to better describe the current patterns of key ecosystem properties like aboveground productivity.
A few modelling studies that implement a similar traits continua approach have recently been published.Scheiter and Higgins (2009) were the first to develop an individualbased framework that eschews the use of PFTs and allows for plants to allocate carbon as a function of local environmental conditions.Falster et al. (2011) presented a model where they used leaf economic strategy, height, wood density and seed size to scale up from individual-scale processes to landscape predictions.Pavlick et al. (2013) applied an interesting approach where they used 15 traits to incorporate trait diversity within a plant community in a DGVM.The rationale of the above models is that they allow different plant functional strategies to be available in a specific location with given environmental conditions (for example a grid cell), and that by setting up a set of functional trade-offs they "filter out" poorly adapted trait combinations from the community.This is effectively an implementation of ideas arising from the environmental filtering/community assembly theory to predict an optimum plant community at a given location (Keddy, 1992;Scheiter et al., 2013;Fortunel et al., 2014).By contrast, drawing on recent findings on the processes controlling Amazonian forest dynamics, we have here attempted to incorporate within TFS the relevant observed associations between functional trait diversity, stand-structure and soil physical and chemical properties (Fyllas et al., 2009;Quesada et al., 2012).Although there are similarities with some of the more recent models discussed above, to our knowledge this is the first time all these linkages have been represented in a single modelling framework.Our approach has been made possible (and thus differs from others) because of the type and quantity of observational constraints used.For example in any given plot we do not force the model to select some "optimum" trait combination based on the prevailing environmental conditions, but we rather assume that the observed trait distribution reflects that of the evolutionarily stable community structure occurring at each site.Similarly we do not require the model to predict what the optimum treesize class distribution would be.Rather, we initialise simulations with what is observed.We have here employed this implementation as our primarily aim in this first instance has been to validate the predictive ability of the model at some extensively monitored Amazonian plots.
Even with these prescribed constraints, the trait randomisation exercise yielded some interesting outputs regarding the importance of trait variability in simulations of forest dynamics.As already discussed, the default variable-trait (vartr) simulations gave the best TFS performance in terms of predicting patterns of aboveground production at the 40 permanent measurement plots with fixed trait (fix-tr) TFS simulations showing a lower predictive ability and an overall higher mean AN .This pattern of trait variability reducing aboveground biomass is in contrast to a similar simulation by Scheiter et al. (2013), where variable trait simulations gave rise to a higher mean biomass because of an increased chance of selecting a trait combination allowing trees to grow larger.This difference arises from the photosynthesis NP co-limitation constraint hardwired into the current version of TFS, as the use of the Amazon-wide mean N L and P L values leads inevitably to universally phosphorus limited estimates of V Cmax and J max that reduce the overall predictive ability of the model.Indeed, when the NP colimitation is removed, the variable-trait simulations (randtr-N) do actually yield the highest AN estimates.Finally, the random variable trait setup (rand-tr) resulted again in a similarly poor TFS behaviour (R 2 = 0.29), emphasising the importance of potential environment-trait interactions in accounting for between-stand structural differences.In other words, in the modelling tropical forest dynamics it is clear that trait distributions cannot be used without a consideration of how they may be shifted by the local growing conditions.
The four functional traits used in the current version of TFS, i.e. leaf dry mass per area, leaf nitrogen and phosphorous concentrations and wood density, are directly related to the rates of tree photosynthesis and respiration.For that reason they provide a stable basis that should allow alternative ecological strategies based on well-known trade-offs such as the "growth vs. survival" to be implemented in trait-based vegetation dynamics models.These four traits have been extensively studied around 70 plots in the Amazon and their patterns of variation and inter-correlation have been analysed (Baker et al., 2009;Fyllas et al., 2009;Patiño et al., 2009Patiño et al., , 2012)).For the purpose of this study, it is important to know the within-stand variation of the functional traits used, i.e. the trait values at the individual level across different plots.Additional functional traits that were considered but not included as base traits in this version of TFS were the seed size and the ratio of leaf area to sapwood area.Seed size is an important functional trait that expresses a tolerance vs. fecundity trade-off, with seed size trading-off with seed number and with larger seed species being more tolerant at more stressful places (Muller-Landau, 2010).However, data on seed size are usually available at the species level, i.e. intraspecific variation is not usually recorded, and thus this kind of data cannot be included in the current version of TFS.The ratio of leaf area to sapwood area, LS , is an important trait that can be used to constrain the hydraulic architecture of trees (Meinzer et al., 2008).Here LS is expressed as a function of D W and H (Eq. 7) and it is not used as an independent (base) trait.Future versions of TFS will include this aspect of functional variability, but for this first study we have selected just a small set of key traits in order to maintain a relatively simple model structure.
Like most modelling efforts, TFS represents work in progress.We identify three particularly promising avenues for future improvements.Firstly, discrepancies between the observed and simulated stem level growth rates, particularly in larger size classes, could result from the allometric equations used to estimate aboveground biomass and growth not being species or size specific.The allometric equations used here express a generic height (H vs. d relationship for Amazonia, without taking into account habitat and species differences, so a more accurate representation of tree architecture would probably result in better biomass growth estimation.Indeed, H −d relationships do vary significantly among species (King, 1996;Poorter et al., 2006) and across regions (Nogueira et al., 2008;Feldpausch et al., 2011;Goodman et al., 2014).An additional source of bias when estimating stem-level growth rates could be related to the uniform (static) allocation coefficient used in this study.For example, Litton et al. (2007) showed that allocation to aboveground tree biomass components increases with age and the availability of resources.Furthermore, Castanho et al. (2013) improved the predictions of a DGVM by adjusting allocation coefficients based on soil texture.Such ontogenetic and/or resource-based shifts in patterns of carbon allocation could be potentially modelled through the use of dynamic allocation schemes (Friedlingstein et al., 1999;Franklin et al., 2012).
The importance of realistically representing autotrophic respiration processes in models of vegetation dynamics is also highlighted here.Modelling respiration has proven to be a difficult task (Cannell and Thornley, 2000), and accurate representation of this component is of great importance for understanding the global C cycle (Valentini et al., 2000).For example, the way that respiration is represented in DGVMs could have a substantial influence over the way that the dynamics of Amazonian forest under scenarios of climatic change are simulated (Huntingford et al., 2004;Galbraith et al., 2010).Nitrogen content of plant tissue has been proven a good predictor of respiration rates (Reich et al., 2008).However, Mori et al. (2010) suggested a mixed-power equation in which the exponent varies from 1 to 3/4 as size increases.Both the Reich and Mori models are implemented in TFS, but we found that a third method, combining the size and nitrogen control, performed better.Thus we suggest that an amalgamation of those two approaches could provide a better way to estimate respiration fluxes in the new generation of dynamic vegetation models.In addition, leaf phosphorous content seems to constrain respiration rates more strongly than nitrogen content in some tropical forests (Meir et al., 2001;Meir and Grace, 2002), and thus inclusion of a phosphorus constraint in future equations of leaf respiration could increase their realism.
Finally, discrepancies in the observed versus the simulated canopy conductance G C could result from the parameterisation of the stomata conductance model of Medlyn et al. (2011).The estimates for g 0 and g 1 used in the 40 PM plots simulations were taken as constant.However, Medlyn et al. (2011) suggested that g 0 and g 1 could vary with functional group.Thus the Amazon-wide parameterisation used here should be replaced with local-level estimates when appropriate gas exchange data are available, and ultimately with estimates based on linked functional traits as evidenced through recently documented associations between structural characteristics such as wood density and leaf area to sapwood ratio with leaf physiological traits such as M a and leaf 13 C / 13 C ratio (Patiño et al., 2012), although we also note that the extent of such structural/physiological linkages remains the subject of debate (Baraloto et al., 2010b).Alternative stomatal closure equations as a function of soil water availability (Harris et al., 2004) should also be tested along with the conductance model in future versions of the model.

Conclusions
We set out to develop a modelling framework for tropical forests that is relatively simple yet adequately complex to capture the main ecological gradients in the world's most extensive tropical forest.Our study places special emphasis on processes highlighted by recent field studies to strongly influence Amazonian forest dynamics, especially functional trait diversity and its association with multiple soil properties (Fyllas et al., 2009).In summary, TFS is characterised by a relatively simple setup, which is capable of reproducing water and carbon fluxes as observed at both daily and multiannual timescales.TFS represents an important link between inventory data and large-scale models with the incorporation of the continuum of plant strategies, through the inclusion of trait distributions providing a step towards better representing diversity in vegetation modelling (Lavorel et al., 2007), representing important processes and trait variation that cannot be adequately accounted for by a DGVM approach to vegetation modelling.Since TFS is based heavily on measured data, the model is well suited to testing hypotheses related to the present-day Amazon biogeography and biogeochemical fluxes.

Figure 1 .
Figure 1.The five basic components of the model and information flow among them.Tree by tree traits and size initialisation takes place at the beginning of each simulation.Carbon and water fluxes, as well as gross and net primary productivity are estimated daily.

Figure 2 .
Figure 2. Geographic distribution of study sites.Dark grey triangles indicate the three eddy flux tower sites (with local names), light grey circles indicate the seven intensive measurement plots (with plot codes), and crosses indicate the coordinates of the 40 RAINFOR permanent measurement plots.
, are used to (a) validate the ability of the model to accurately simulate stand-level carbon fluxes and (b) explore patterns of G , N and C U along the Amazonian climatic and soil nutrient availability gradient.The size-class distribution within each PM site is used to initialise the stand structure of the model and simulate patterns of productivity for the 2000-2006 period.

Figure 3 .
Figure 3. Simulated against observed mean daily canopy conductance G C for the three sites with eddy flux data.The broken line represents a 1 : 1 relationship and the continuous line illustrates a standardised major axis (SMA) regression.

Figure 4 .
Figure 4. Simulated stem growth rate N,s against observed aboveground biomass change B ABG for different size classes for the 2000-2006 period.Upper panel: infertile plots.Lower panel: fertile plots.The broken line represents a 1 : 1 relationship.The continuous line illustrates the straight line fit using the York method (see text for details).

Figure 5 .
Figure 5. Simulated stand-level aboveground net primary productivity ( AN ) against observed stand-level aboveground biomass growth ( B ABG ) of surviving trees, at the 40 PM plots.The line illustrates a SMA regression of α = 0.92(0.72, . . ., 1.18) and R 2 = 0.42.Red dots indicate high nutrient availability and blue dots indicate low nutrient availability plots.

Figure 6 .
Figure 6.Summary of the randomisation exercise simulations.(a) Simulated stand-level aboveground net primary productivity ( AN ) against observed stand-level aboveground biomass growth ( B ABG ) for the four different setups.The slope of the SMA (a) and the adjusted R 2 are given in parentheses for each setup.Different colours indicate different setups.(b) Simulated Amazon-wide aboveground net primary productivity ( AN ) for the four different setups.