Individual-Based Modeling of Amazon Forests Suggests That Climate Controls Productivity While Traits Control Demography

Climate, species composition, and soils are thought to control carbon cycling and forest structure in Amazonian forests. Here, we add a demographics scheme (tree recruitment, growth, and mortality) to a recently developed non-demographic model—the Trait-based Forest Simulator (TFS)—to explore the roles of climate and plant traits in controlling forest productivity and structure. We compared two sites with differing climates (seasonal vs. aseasonal precipitation) and plant traits. Through an initial validation simulation, we assessed whether the model converges on observed forest properties (productivity, demographic and structural variables) using datasets of functional traits, structure, and climate to model the carbon cycle at the two sites. In a second set of simulations, we tested the relative importance of climate and plant traits for forest properties within the TFS framework using the climate from the two sites with hypothetical trait distributions representing two axes of functional variation (“fast” vs. “slow” leaf traits, and high vs. low wood density). The adapted model with demographics reproduced observed variation in gross (GPP) and net (NPP) primary production, and respiration. However, NPP and respiration at the level of plant organs (leaf, stem, and root) were poorly simulated. Mortality and recruitment rates were underestimated. The equilibrium forest structure differed from observations of stem numbers suggesting either that the forests are not currently at equilibrium or that mechanisms are missing from the model. Findings from the second set of simulations demonstrated that differences in productivity were driven by climate, rather than plant traits. Contrary to expectation, varying leaf traits had no inﬂuence on GPP. Drivers of simulated forest structure were complex, with a key role for wood density mediated by its link to tree mortality. Modeled mortality and recruitment rates were linked to plant traits alone, drought-related mortality was not accounted for. In future, model development should focus on improving allocation, mortality, organ respiration, simulation of understory trees and adding hydraulic traits. This type of model that incorporates diverse tree strategies, detailed forest structure and realistic physiology is necessary if we are to be able to simulate tropical forest responses to global change scenarios.


INTRODUCTION
The Amazon forest is the most extensive tropical forest on Earth.Storing an estimated ∼100 Pg C in aboveground woody biomass (Feldpausch et al., 2012), this region accounts for approximately 25% of the global aboveground forest carbon stock (Pan et al., 2011).The Amazon is also hugely diverse, containing many thousands of tree species (ter Steege et al., 2013;Cardoso et al., 2017).Understanding the sensitivity of Amazonia to climate change and anthropogenic disturbance is therefore a research priority, which requires a mechanistic understanding of the processes driving the carbon cycle.
A wide variety of processes act at multiple spatial and temporal scales to influence the Amazon forest carbon cycle.Plotbased measurements of woody biomass, mortality, and species composition across the basin indicate that there are remarkably covarying spatial patterns in forest structure, dynamics and biodiversity.There is a north-east to south-west gradient of declining biomass (Malhi et al., 2006;Quesada et al., 2012;Mitchard et al., 2014;Johnson et al., 2016) coincident in space with declining community mean wood density (WD, Baker et al., 2004;Malhi et al., 2006;ter Steege et al., 2006;Quesada et al., 2012) and increasing tree mortality rates (Phillips et al., 2004;Quesada et al., 2012;Johnson et al., 2016).Soil physical properties are thought to underlie these patterns; the younger, shallower soils of the south-west promoting a composition of short-lived, fast-growing, low WD species due to their poorer structural properties limiting root anchorage (Quesada et al., 2012).Stem mortality rates are not only an important driver of spatial patterns of forest biomass (Johnson et al., 2016), but also temporal patterns of carbon uptake (Brienen et al., 2015).At the tree-level, mortality is linked negatively with WD (Muller-Landau, 2004;King et al., 2006;Kraft et al., 2010;Aubry-Kientz et al., 2013) and positively with growth rate (Chao et al., 2008).Stand-level woody productivity is only weakly related to biomass (Johnson et al., 2016) but may be linked to soil fertility, especially phosphorus (Aragão et al., 2009;Quesada et al., 2012).Drought events can cause forests to shift temporarily from acting as carbon sinks to sources through episodic impacts on tree growth and mortality (Phillips et al., 2009;Gatti et al., 2014;Doughty et al., 2015;Feldpausch et al., 2016) and may be behind a longer-term weakening of the sink (Brienen et al., 2015).Additionally, Malhi et al. (2015) show that variation in gross primary productivity (GPP) between sites is largely driven by variation in the seasonality of precipitation, but this does not follow through to a relationship between climate and net primary productivity (NPP), nor woody productivity, due to differences in autotrophic respiration (Ra) and allocation between sites.
The ecosystem properties of a forest (e.g., biomass, GPP) are the aggregated sum of the individual trees in that forest, so by understanding the life cycle of individual trees and harnessing physical understanding of plant physiology, it is possible to build individual based models from first principles to simulate ecosystem processes (DeAngelis, 1988;Grimm and Berger, 2016;Shugart et al., 2018).In forests, demographic processes (growth, mortality) are linked to size (Lines et al., 2010;Herault et al., 2011;Aubry-Kientz et al., 2013;Bigler, 2016), and hence the sizestructured nature of individual-based models allow for inclusion of size-dependent variation in processes.In diverse tropical forests, species show differences in performance linked to their functional traits (e.g., Herault et al., 2011;Ruger et al., 2012).In recognition of these results, recent advances in vegetation modeling have developed new trait-based and individual-based mechanistic models (aDGVM2 Scheiter et al., 2013;TFS Fyllas et al., 2014;LPJmL-FIT Sakschewski et al., 2015) which enable tree size and functional traits to be explicitly combined for the simulation of forest processes.This is particularly useful for tropical forests where the enormous species diversity precludes the use of species-specific models.Models with the inclusion of trait-continua can represent diverse plant strategies, and move beyond the traditional approach of plant functional types (PFTs) which use a limited number of prescribed trait combinations to describe plant strategies.
The Trait-based Forest Simulator (TFS, Fyllas et al., 2014) differs from aDGVM2 and LPJmL-FIT by including locally observed distributions of key functional traits accounting for variation in tree life cycle and productivity [WD, leaf mass per area (LMA), leaf nitrogen (N), and phosphorus (P) concentrations].This allows TFS to indirectly account for the influence of environmental factors on species composition and leaf nutrients.The use of observed distributions also allows the extent of local co-variation between traits to be accounted for, without assuming that co-variation occurs, nor assuming any underlying theoretical distribution.Trees with high WD typically grow more slowly, live longer and store more biomass for a given size than those with low WD.LMA describes the resource investment in the leaf and is a crucial component of the leaf economics spectrum (Wright et al., 2004), while leaf nutrient concentrations in conjunction with LMA are linked to photosynthetic capacity (Domingues et al., 2010); carbon assimilation per leaf area is expected to increase with LMA and leaf nutrient concentrations.TFS is a non-demographic model that uses inputs of forest structure (tree diameter distribution), trait distributions, local or gridded climate data, and soil type to simulate carbon and water fluxes for specific tropical forest plots.With this array of inputs, the model can account for many of the processes thought to influence Amazonian forest carbon cycling, but currently lacks a dynamic scheme (growth, mortality and recruitment).Here, we add demographics to this model to produce a fully dynamic version of TFS (TFSv.2, Figure 1; Figure S1).
To add growth to the model, we include allocation of NPP.Some typical model approaches to allocation are to maintain the originally specified allometries of the tree (e.g., ED, Moorcroft et al., 2001), to use fixed allocation coefficients (e.g., TROLL, Marechaux and Chave, 2017) that can vary with PFT (e.g., ED2, Medvigy et al., 2009), or to consider allocation itself as a trait (e.g., aDGVM2).In TFSv.2 we add the important trade-off between LMA and leaf lifespan (Wright et al., 2004;TROLL, LPJmL-FIT), where trees with higher LMA also have a longer leaf lifespan to account for the greater carbon investment.A consequence of differing leaf turnover rates is that allocation to leaf biomass may also be different.Hence, we allocate carbon between organs based on the losses from these organs, which includes variation due to LMA for leaves.Most individual-based models represent mortality as a function of PFT specific longevity or maximum size (e.g., SEIB-DVGM, Sato et al., 2007;Sato, 2009;LPJ-GUESS, Smith et al., 2001;ED2;FORMIND, Köhler and Huth, 1998;Fischer et al., 2016), and growth rate or carbon balance (e.g., SEIB-DGVM, LPJ-GUESS, ED2, FORMIND, TROLL, SORTIE, Pacala et al., 1996).Here, we include a mortality scheme based on an empirical relationship predicting the probability of mortality from growth rate (an outcome of traits and environment), wood density (to represent life history) and tree size as a fraction of maximum potential size parameterized from a large dataset of Amazon forest plots.For the establishment of new trees, models often produce seeds or saplings from a reproductive biomass pool (e.g., LPJ-GUESS, ED, aDVGM2).PFT specific minimum reproductive threshold diameters are occasionally used (e.g., TROLL, FORMIND) and production of propagules can include PFT specific dispersal distances (e.g., FORMIND, SORTIE).Sapling establishment is typically dependent on light availability, sometimes with PFT specific light requirements (e.g., FORMIND, LPJ-GUESS, SEIB-DGVM).For this version of TFS (where the minimum size of simulated trees is 10 cm diameter) we propose a simple recruitment scheme where all adult trees have the same chance of reproducing and the number of saplings is dependent on light availability.
We use this updated fully dynamic model (Figure 1) to simulate the carbon cycle, forest structure and population dynamics of two contrasting lowland western Amazon sites, Tambopata (TAM, Malhi et al., 2014) and Allpahuayo (ALP, del Aquila-Pasquel et al., 2014).TAM has a much more seasonal climate than ALP (Figure 2), while soils and leaves in TAM have higher nutrient concentrations than ALP (Table 1).These differences allow us to test the relative importance of climate seasonality and plant traits for Amazon ecosystem function.Furthermore, the sites are home to extensive data collection efforts in forest dynamics and detailed measurements of the carbon cycle, enabling far-reaching evaluation of the model against data.We are investigating, through the mechanisms proposed within TFS v.2, the interplay of traits and climate in driving carbon cycling and forest structure, and the reasons differing climates and traits can influence modeled outcomes.The aims of the paper are to assess (i) whether TFSv.2 converges on forest properties (productivity, demographic and structural variables) within the range of observations, (ii) the impact of climate seasonality and plant traits on forest processes with the TFSv.2 framework, and (iii) the relative importance of leaf and stem traits for forest properties.For the latter two aims, we perform experiments using hypothetical trait distributions extracted from the observed distributions, and the different climates of the two study sites.While no model can be perfect, it reflects our current understanding.For aspects that are well simulated, the model may suggest the mechanisms, and for aspects that are not well simulated, the model highlights where improvement in understanding is needed.

Site Description
Within each site, Allpahuayo and Tambopata, two 1-ha forest inventory plots (ALP-01, ALP-30, TAM-05 and TAM-06) have been established (Table 1).The climate in ALP (mean annual precipitation 2,689 mm) has little seasonality with no dry season and relatively constant temperature, whilst in TAM (mean annual precipitation 1,900 mm) there is a reduction in rainfall below 100 mm during 5 months of the year which coincides with a drop in temperature of approximately 2 • C (Figure 2).Substrate varies between and within sites.Soils in TAM are more fertile than ALP, with particularly high total phosphorus in TAM-06, hence the higher leaf P concentrations in TAM-06 (Table 1).Soils in ALP have a high sand content, while those in TAM are more clay rich.Although ALP-30 has a very high sand content and low fertility, fast-draining substrate, tree roots may be able to access deeper, higher fertility soil below (Malhi et al., 2015).To our knowledge, the sites have not been subject to any recent human impacts.

Observations
We use published data (del Aquila-Pasquel et al., 2014;Malhi et al., 2014Malhi et al., , 2015) ) from intensive carbon cycling measurements following the GEM (Global Ecosystems Monitoring) protocol (Marthews et al., 2014) in which wood (NPP W , including coarse roots), leaf (NPP L ) and fine root NPP (NPP FR ) (measured using dendrometer bands, leaf litter production and root ingrowth cores, respectively) and Ra (R W , R L , R FR ; measured using gas analysers on stem surfaces, individual leaves, and soil surfaces and then scaled to the whole plot) are monitored.R FR is separated from soil respiration by comparing measurements from cores where roots are excluded to those with roots included.The NPP and Ra from each organ is summed to estimate the plot-level GPP (GPP = NPP + Ra).All are presented as Mg C ha −1 yr −1 .Full methods used for the carbon cycling measurements can be found in the appropriate references (ALP: del Aquila-Pasquel et al., 2014;TAM: Malhi et al., 2014).Because forest demographic rates can vary over small spatial and temporal scales, we calculated stem mortality, biomass mortality and recruitment rates (% yr −1 ), woody productivity rates (Mg C ha −1 yr −1 , differing from NPP W due to exclusion of coarse roots) and mean diameter growth rates (cm yr −1 ) from a larger dataset from inventory plots located within the same sites as the intensive plots in order to include a broader array of observed dynamic rates (Table S1).We refer to these larger plot datasets as clusters.Inventory data were collected following the RAINFOR protocol (Phillips et al., 2016), whereby all stems ≥ 10cm diameter at breast height (D, cm) are botanically identified and have D measured every 1 to 5 years over several decades.These data are collated in the RAINFOR ForestPlots.netdatabase (Lopez-Gonzalez et al., 2009, 2011).We excluded nearby plots on waterlogged soils (ALP-40, MSH-01, TAM-03) from the clusters.Error statistics for the field data are described in Supplementary Methods S1.

Model Description and Development
Non-demographic Version TFS v.1 TFS uses traits and tree size in the estimation of carbon balance of individual trees.The observed diameter distribution of a forest plot is used to initialize the model, and each input diameter becomes an individual tree.Each tree is randomly assigned a suite of four traits [WD (g cm −3 ), LMA (g m −2 ), N (mg g −1 ), P (mg g −1 )] from the observed distribution (retaining any covariation between traits in the observed distribution).Allometry is used to  define the individuals' dimensions and biomass pool sizes based on diameter and traits.Trees compete for light using the Perfect Plasticity Approximation (Purves et al., 2007) whereby crown area and height determine the canopy layer in which the tree is situated and the amount of light received.Soil water availability is estimated daily from climate data with a simple bucket model using parameters based on soil type.The leaf traits are used to estimate the photosynthetic parameters Vc max (the maximum rate of carboxylation, µmol m −2 s −1 ) and J max (the maximum rate of electron transport, µmol m −2 s −1 ) (Domingues et al., 2010), which together with the light received and a water stress function determine the photosynthetic rate (Farquhar et al., 1980) and GPP.Soil moisture and respiration are simulated on a daily time-step, photosynthesis is simulated on an hourly timestep.Plot-level fluxes are determined by summing the fluxes of each individual tree within the plot.For full details of TFS v.1 see Fyllas et al. (2014).

Developments
A small number of alterations made to improve the realism of the model are presented in the Supplementary Methods S2, along with all new equations and parameter values.Briefly, the spatial scale on which canopy light competition is performed was reduced to 20 × 20 m subplots to improve representation of within plot heterogeneity in canopy height (Bohlman and Pacala, 2012).Leaf lifespan (LL, days) was allowed to vary and assigned to individuals based on their LMA following the known relationship from the leaf economics spectrum (Wright et al., 2004).Further, LL was increased for subcanopy trees based on observations.R L (modeled as a function of Vc max ) was reduced for shaded trees.

Tree mortality
The annual probability of mortality (P m ) was estimated as a logistic function of WD, diameter growth rate (Gr, cm yr −1 ), and D relative to maximum diameter (D max ) (Equation 1, Table S2).
This approach was used to determine mortality dependent both on the actual circumstances of the tree (low growth rate producing a high probability of mortality, long-standing negative carbon balance resulting in death), the life history strategy of the tree (high WD being associated with low mortality), and how close the tree is to its D max (trees closer to their D max being more likely to die due to senescence).
Data from 69 forest plots located in the western Amazon from the RAINFOR ForestPlots.netdatabase (Lopez-Gonzalez et al., 2011) were used to fit the model to data using maximum likelihood applied to a binomial process, see Supplementary Methods S3.D max (95th percentile diameter) was estimated for all species in the dataset with at least 20 individuals, following the method of King et al. (2006).We used data for all individual trees present (alive or dead) in at least three censuses and for which a D max value was available.Maximization of the likelihood function was performed using simulated annealing, following similar methods to Chao et al. (2008), Lines et al. (2010), andAubry-Kientz et al. (2013).Different models with all possible combinations of the three predictors were compared using AIC, showing that the full model with all predictors had the greatest support (Supplementary Methods S3).WD was the best single predictor, and WD combined with Gr was the best model with two predictor variables.Relationships between the variables and the probability of mortality are shown in Figure S2.
To apply the equation in TFS v.2 a D max value was assigned to each simulated tree.From observations (Figure S3), there was a weak relationship between species LMA and D max (and no relationship with other traits).Whilst the relationship between LMA and D max is not strong, the (log-normal) distributions of D max binned by LMA shift to the right (i.e., increase) at higher LMA limits (Figure S3).To maintain the distribution of D max values found for Amazonian species and its weak relationship to LMA within the simulated communities, a D max value for each tree was taken from the distribution of observed D max values in 5 classes of LMA (Figure S3) with LMA from Fyllas et al. (2009).In the model, mortality is performed on an annual timescale.After calculating P m for each simulated tree using equation 1, death was determined by a simulated random coin toss with a uniform probability function.In addition, trees with a negative Gr for five consecutive years were killed.

Tree recruitment
Each year a maximum number of saplings (N sap ) are established in each 20 × 20 m subplot.N sap is controlled by the total leaf area index (LAI, m 2 m −2 ) of the subplot with lower sapling densities at higher LAI following a negative exponential function (Fyllas et al., 2010;equation 2), where A is the subplot area (m 2 ).Parameter a represents the maximum potential number of saplings in one m 2 .Trait assignment to saplings follows the same method of random vector generation used at initialization (Fyllas et al., 2014), but based on the traits distributions of the mature trees of the stand rather than the input distributions.The saplings are maintained in a pool of potential recruits for 10 years with a 90% annual chance of survival (within the range of sapling survivorship observed over 3 years in Panama, Welden et al., 1991) before being established as mature trees with D = 10 cm.

Closing the carbon cycle
The biomass pool sizes of leaves, stem, fine, and coarse roots were updated annually based on the NPP produced over the year and losses (leaf, root and branch shedding) removed over the year.
Losses of coarse and fine roots were not included in TFS v.1 and so were added into the new dynamic version.Fine root lifetime was set at 6 months (Jiménez et al., 2009).Coarse root losses were set equivalent to branch fall with the assumption that coarse root biomass also contains a "labile" component equivalent to the branch component of stems.Given that some individuals would have higher leaf losses than others due to differing LL, it follows that the allocation to leaf biomass should also vary, depending on the tradeoff between LL and LMA.To account for this, allocation was carried out in two stages.First, all losses are replaced using the NPP.The surplus NPP available after replacement of loses is then allocated to wood, root, and leaf biomass pools using static coefficients.Because the majority of losses are derived from roots and leaves, we hypothesized that most of the surplus NPP is allocated to wood (96%) with 2% going each to fine roots and leaves.With the updated biomass values, plant dimensions (D, height, crown area, and leaf area) were updated.Individuals with NPP less than total losses allocate the available NPP to biomass pools in proportion to the amount of losses from each pool, but the pools are not fully replaced and the tree has a negative Gr.

Simulations
Validation Simulations The forest inventory data from the first census of each plot were used to initialize the model (each diameter becoming an individual tree).Functional traits were taken from Fyllas et al. (2009) and the GEM-TRAIT project (Supplementary Methods S1, Figure S4).Soil depth and texture were taken from Quesada et al. (2010).Soil depths are similar (2.2-2.4 m) in TAM and ALP-30, but shallower in ALP-01.We ran the model for 500 years to allow ample time to reach an equilibrium forest structure, beginning in 1983 when the first forest inventory data were collected.Gridded climate data with a daily time-step for the period 1983-2012 were taken from Sheffield et al. (2006).The climate data is comparable with local data from the ALP site (del Aquila-Pasquel et al., 2014; Figure 2A); however, it has higher precipitation than the local data at the TAM site (Malhi et al., 2014; Figure 2B), with a shorter dry season.For simulation years post-2012, the 1983-2012 climate was recycled by randomly selecting a year of data.Atmospheric CO 2 was kept constant at 380 ppm.Due to inherent stochasticity in the model, 10 runs were performed.The mean of the last 100 years of runs was used to compare with the observed data, with errors calculated as standard deviation from the mean across the 10 runs, accounting for all stochastic processes.

Datasets for Model Validation
We compared observed and simulated GPP, NPP, Ra, and NPP and Ra by organ at the plot-level, and calculated the normalized mean bias (NMB, Equation 3) to evaluate the datamodel agreement where S i is the simulated and O i the observed value for n plots.In Malhi et al. (2015) canopy NPP is separated into four components (leaf, twig, fruits and flowers, and herbivory); here we combined leaf production and herbivory for NPP L .For NPP W we include NPP stem of trees > 10 cm D, NPP branch, and NPP coarse root.We compared simulated demographic rates and forest structure properties against the observed forest inventory data.Observed and simulated distributions (e.g., size class distributions) were compared with Wilcoxon rank sum tests.

Experimental Simulations
In a second set of simulations to test the impacts of trait distributions and climate on forest processes we ran the model under a range of climates and trait distributions.We use the climate data for the two sites (Figure 2).As the trait distributions from the plots overlap (Table 1; Figure S4) we generated four different trait distributions that varied systematically to use as in model experiments (Figure S5).To do this, we combined all the trait data from all plots and separated it in two ways.First, the samples were split between those with N higher or lower than 20 mg g-1.This produced two opposing traits distributions; higher fertility, lower LMA leaves-"fast" leaf traits, as opposed to lower fertility, higher LMA leaves-"slow" leaf traits.Due to its link with LMA, LL was shorter in simulations with "fast" leaf traits and longer in simulations with "slow" leaf traits.Mean WD was higher for the "slow" leaf traits distributions than the "fast" leaf traits distributions but still showed wide variation and overlapped between the two distributions.Second, the samples were split between those with WD value higher or lower than 0.6 g cm −3 .The N, P and LMA values were similar between high and low WD distributions.Simulations were performed using the ALP and TAM climates and each trait distribution (2 climates x 4 trait distributions = 8 combinations), with 10 runs of each.The simulations were initialized with 200, 10 cm D trees to examine the influence of traits and climate on stand development as well as equilibrium forest structure and carbon cycling.We used two-way ANOVA to determine the percentage of variance in simulated variables explained (ω 2 ) by altering the input climate (ALP vs. TAM) or the input traits ("fast" vs. "slow" leaf traits, or high vs. low WD).
All statistical analyses were performed using R 2.15.1 (R Core Team, 2012).

Validation Simulations of Carbon Cycling, Forest Structure, and Demographics
Simulations of the carbon cycle show good agreement with GPP (Figure 3A; Table 2, NMB = −2.6%)and Ra (Figure 3C; Table 2, NMB = 0.6%) observations, where GPP and Ra were higher in ALP than TAM.The NMB of simulated NPP was a little higher (Table 2, NMB = −8.5%),and the simulated values were higher for ALP than TAM (due to a slight underestimation of GPP for TAM-05 and slight overestimation of Ra in TAM-06), whereas the observations were similar for both sites (Figure 3B).
Examining NPP by organ (Figures 3D-F; Table 2), NPP FR was overestimated (NMB = 53.4%)whilst NPP W and NPP L were underestimated (NMB = −29.5% and −18.6%, respectively).The observed data show variation in allocation of NPP to different organs between plots, which is not reflected in the model simulations (Figure S6).Simulated respiration by plant organ shows similar patterns to NPP by plant organ, but with higher biases and overestimation of R L (Figures 3G-I; Table 2).
The final steady-state forest structure at the end of simulations is the result of 500 years of drift from the initial (observed) structure due to recruitment, growth and mortality.The final simulated size-class distributions (combined from all runs) for all plots were significantly different from the initial distribution (Figure 4, Wilcoxon rank sum test p < 0.0001).The percentage of large stems (>30 cm diameter) was higher in simulations (33.6-41.9%)compared to observations (13.9-20.1%).
For all plots, stem numbers initially decreased and reached a steady-state after c. 100-200 years (Figures 5A-D), with stem numbers much lower than observations (NMB = −47.6%,Table 2).Simulated BA initially increased and then reduced in all plots reaching a steady state which was below the cluster mean for TAM plots and similar to the cluster mean for ALP plots (Figures 5E-H; Table 2).Equilibrium above-ground biomass (AGB) was simulated more accurately (NMB = 5.5%, Table 2).Observed temporal patterns of stem numbers and basal area are shown in Figure S7.
The distributions of individual tree growth rates were significantly different between simulations and observations for all plots with growth somewhat overestimated in simulations (Figure 6, Wilcoxon rank sum test, p < 0.0001); however, growth rates were within realistic ranges.A higher proportion of trees were simulated to grow between 0 and 0.1 cm yr −1 than in the data, but the data also showed a number of trees that shrank (from reductions in diameter of stressed trees or measurement error).The simulations produced a higher proportion of trees (14.3-26.1%)growing in 0.5-1 cm yr −1 than in the data (8.2-11.7%),especially in the ALP plots.
Observed distributions of woody productivity were not significantly different between the ALP and TAM clusters (Figure S10, Wilcoxon rank sum test, p = 0.13), whilst the simulated distributions were significantly different between plots (Wilcoxon rank sum test, p < 0.0001), with higher woody productivity in ALP plots (Figure S10).Compared to the observations, distributions of woody productivity for ALP plots were not significantly different (Wilcoxon rank sum test, p > 0.06), while the simulated productivity rates of TAM plots were underestimated (Wilcoxon rank sum test, p < 0.0001).

Climate and Trait Influence on Experimental Simulations
The influence of climate and traits on simulated GPP, NPP, Ra, stem number, BA, AGB, woody productivity, LAI and demographic rates was tested.Climate was the predominant driver of simulated GPP, NPP, Ra and stem number in both trait experiments ("fast" vs. "slow" leaf traits [leaf traits experiment] and high vs. low wood density [WD experiment], Table 4; Figure 7; Table S3) with higher GPP, NPP and Ra, and lower stem number under the ALP climate.LAI was markedly influenced by climate for the leaf traits experiment (proportion of variance explained ω 2 = 0.91), with higher LAI under the ALP climate, but traits had a stronger role in the WD experiment (ω 2 = 0.23), with higher LAI under the high WD traits distribution.Although climate was the stronger driver of NPP in both experiments, the importance of the traits distribution was higher for the leaf traits experiment with higher NPP with "slow" leaf traits.AGB was driven mostly by climate, but with a role for traits in both experiments (ω 2 = 0.20).The results for BA in the leaf traits experiment was similar to AGB, but in the WD experiment traits became much more important (ω 2 = 0.75, Figures 7e-l; Table 4).Interestingly, in the WD experiment the high WD distribution lead to lower BA but higher AGB compared with the low WD distribution, while in the leaf traits experiment the "slow" distribution was associated with higher BA and AGB compared with the "fast" distribution.Climate was the major driver of woody productivity for the leaf traits runs, but traits also had a role in the WD runs (ω 2 = 0.22).Demographic rates (stem mortality and recruitment) were driven by trait distributions in both cases, but with a much stronger effect for the WD experiment (Table 4, higher rates with lower WD).Interactions were rarely significant, and explained very little variance (ω 2 ≤ 0.02).
The simulations showed similar variation in stand development under all experimental conditions.In all circumstances, stem numbers increased quickly, then declined to an equilibrium level (Figures 7a-d).Despite equilibrium stem numbers not being strongly influenced by trait distributions (Table 4), the peak stem numbers reached early during succession were higher for simulations using "slow" leaf traits than "fast" leaf traits (Figures 7a-d).This is likely due to the higher demographic rates resulting from the "fast" leaf trait distributions.WD increased during all simulations with a significantly different mean WD between the first and last years (Table S4; Figure S12).

Simulation of Carbon Fluxes
Our findings suggest that differences in GPP between sites are due to climate.Simulations using the ALP climate had consistently higher GPP even under a range of different trait distributions.This occurs in the model due to a reduction in stomatal conductance as a result of seasonal soil moisture stress in TAM (Figure S13), which limits photosynthesis (Fyllas et al., 2014).Is it realistic that climate is the main driver of GPP differences between these sites?Malhi et al. (2015) found a strong relationship with seasonal water deficit and GPP across a range of sites, suggesting that seasonal water limitation reduces photosynthesis.This is supported by GPP observations from flux towers showing a seasonal reduction in photosynthesis at southern Amazon sites (Restrepo-Coupe et al., 2013).Remote sensing data also suggests a seasonal reduction in photosynthesis is expected at TAM (Guan et al., 2015).However, the role of climate as the sole driver of photosynthesis seasonality has been questioned, with evidence that leaf demography (leaf dropping, flushing, and age related changes in photosynthetic efficiency) determines seasonal changes in canopy photosynthetic capacity, driving seasonality of GPP (Hutyra et al., 2007;Restrepo-Coupe et al., 2013;Wu et al., 2016).In TFS v.2, photosynthetic capacity is constant throughout the year.Hence, the mechanism by which climate influences GPP may be more complex than we propose.

Beyond GPP
Organ-level Ra and NPP are key areas for model improvement and further research.Both Ra and NPP show similar patterns; R W and NPP W are consistently underestimated, while R L , R FR and NPP FR are overestimated.This may be due to the low simulated stem densities and larger tree sizes than observed (Figures 4, 5),    The mean of each variable from the last 100 years was taken for 10 simulations and modeled using anova.Results are presented from the two trait experiments; "fast" vs. "slow" leaf traits and high vs. low WD.
and the lower proportional allocation to wood of large trees.It may also be due to the allocation scheme underestimating NPP W by preferentially replacing lost biomass (mainly leaf litter and root turnover).In any trait-based model it is important to include ecological trade-offs.In TFS v.2, we increased leaf lifespan (LL) with increasing LMA (Supplementary Methods S2; if all trees had the same LL low LMA trees would develop unrealistically dense crowns).Given this difference in LL, it is not possible or desirable to use static allocation coefficients (e.g., those observed from plot-level data e.g., Malhi et al., 2011) for all trees.Trees with a short LL need to allocate more carbon to leaves to maintain their canopy.This brings us to consider allocation as a trait-mediated function that will vary between individuals and species, perhaps as part of a "fast-slow plant economic spectrum" (Reich, 2014).Another aspect of allocation that was not incorporated into the model is a non-structural carbohydrate pool.In some sense, the annual timestep of allocation and growth allows for storage of photosynthate within the year, as during seasonal water stress photosynthesis is reduced but respiration continues.However, there is no consideration of longer-term storage, which has shown to be common across plants globally (Martínez-Vilalta et al., 2016), and could improve the model.Understanding mechanisms for different allocation patterns between species, sites, and under varying conditions remains a challenge, and outcomes from new datasets (e.g., Doughty et al., 2014Doughty et al., , 2015;;Girardin et al., 2016), in addition to organ-level understanding of respiration, should be incorporated into models when available.
The observations of woody productivity rates from forest inventory data (which is different from NPP W as it does not include branch fall or coarse roots, and here observations were based on a larger forest inventory dataset) were not significantly different between plots.In contrast, the simulations matched the observations for ALP but underestimated for TAM.Higher simulated wood production in ALP results from higher GPP due to less seasonality, but this is not supported by the data.For TFS v.2, climate is a stronger driver of woody productivity than traits, but with a small driving role for WD (Table 4).Baker et al. (2009) investigated plot-based observations of woody productivity across the Amazon and found it to be driven by environmental factors rather than differences in WD or maximum size between sites.Data from intensive plots show that while GPP is linked to climate, NPP does not show the same patterns (Malhi et al., 2015).Hence, there is progress to be made in improving our understanding of the linkages between GPP and woody growth.

Simulation of Forest Structure and Demographic Rates
Whilst the total standing stock and volume of vegetation (BA, AGB) is reasonably simulated, there is a lack of model performance with respect to the partitioning within individual stems (stem density, size class distribution, diameter growth rates, Table 2; Figures 4-6).The deviations from the observations were particularly strong in the ALP plots, and in the experimental simulations stem density was controlled only by climate.The higher productivity of the aseasonal ALP climate may enhance the dominance of large stems.The underestimation of stem density does not appear to be linked to the low mortality rates, as the simulation of stem density was unresponsive to altering input traits, while mortality rates increased under the low wood density simulations (Table 4).
To understand this further, a correlation matrix of various simulated stand-level variables was produced to give insight into the internal mechanics of the model relating to stem density.As shown in Table 5 and Figure S14, stem density is tightly negatively correlated to mean tree diameter, a density-size relationship that is commonly observed in forests (e.g., Long and Smith, 1983).Mean tree diameter is strongly positively correlated with GPP and only weakly with mortality rate (Table 5), while GPP itself is related to climate (Table 4), hence stem density being linked to climate.These findings should be taken with caution, as correlation does not mean causation, especially for GPP, and mean stem diameter which are not correlated when only simulations under the same climate are tested.Differences between observed and predicted size class distributions could be due to the observed forests not being in an equilibrium state.However, size class distributions in many forests across the tropics have been studied, and generally show patterns similar to the observed distributions (Enquist and Niklas, 2001).The cause of poor simulation of stand structure may relate to differential growth strategies of understory trees (e.g., Hubau et al., 2019) that may require more differences in parameterization than are already included in TFSv.2.Farrior et al. (2016) show that size class distributions can be simulated by non-physiological demographic models when size dependent rates are included.The experimental simulations produce a mean tree diameter range of 23.2-36.1 cm (Figure S14), higher than the observed range of 19.9-21.9cm from the first census of each plot.The close relationship between simulated stem density and mean tree diameter suggests that to improve our simulation of stem density it is necessary to improve the simulation of mean tree diameter.Over-estimation of large trees (Antonarakis et al., 2011) and growth rates (Medvigy et al., 2009;Medgivy and Moorcroft, 2012) were also a feature of the ED2 model, and for growth were resolved by optimizing poorly constrained parameter values, an approach from which TFSv.2 would also likely benefit.Particular uncertain parameters that may improve model performance if optimized are the allocation coefficients, respiration parameters (including the fraction of V cmax used to estimate leaf respiration at different levels of shading, Equation S4), and the recruitment model parameters (Equation 2).
Simulation of mean tree diameter may be a useful variable for benchmarking the performance of individual and cohort based models (Fisher et al., 2018).
WD had a strong influence on BA and, to a lesser extent, AGB (Table 4).This is through its impact on mortality rates which limit tree lifespans and hence the size that they can reach and the build-up of biomass within an individual.The importance of large trees for stand-level biomass has been observed in many tropical forests (DeWalt and Chave, 2004;Baraloto et al., 2011;Fauset et al., 2015) and our model formulation with individuals of varying size enables this effect to be included.Further, understanding of the link between tree mortality and stand-level biomass has been in development over the past decade (Malhi et al., 2006;Quesada et al., 2012;Galbraith et al., 2013) and has been supported by a large forest inventory dataset (Johnson et al., 2016).As shown by Johnson et al. (2016), typical big-leaf DVGMs are unable to include links between mortality and biomass as a fully described size structure is necessary, lending support for size structured models (Fisher et al., 2018) such as the TFS v.2 model presented here.
Stem mortality was underestimated in all plots (Figure S8; Table 3).The slightly higher biomass mortality rates compared to stem mortality rates is likely due to the fact that simulated tree sizes were larger than observed (Figure 4), so even if fewer trees died than in observations, the biomass content of these simulated trees was higher.Mortality rates were linked more to traits than climate (Table 4), and the higher rates for TAM-06 are due to the high abundance of low WD palm species Iriartea deltoidea and Socratea exorrhiza (Malhi et al., 2014) in the input traits distribution (Table 1).One potential reason for underestimating mortality is an increase in WD values over time in the simulations (Figure S15, a mean difference between the initial and final simulation years of 0.09 g cm −3 across all plots in validation simulations, see also Figure S12), as the probability of mortality (P m ) is strongly linked to WD. Whilst there is a growth advantage for trees of low WD which show faster growth than higher WD trees of similar size in "disturbed" subplots (Figure S11), this is not a sufficient trade-off with increased mortality to maintain the initial prescribed WD distribution throughout the simulations.Hence, the higher mortality of low WD trees causes a shift in the community toward higher WD trees, which itself reduces the community mortality rates as high WD trees have a lower chance of dying.The final WD values in the plot simulations are similar to observed plot mean WD values from the eastcentral Amazon (approaching 0.75 g cm −3 , Quesada et al., 2012).This suggests that the final simulated WD distributions are not unrealistic for forests with low mortality rates.To maintain lower plot mean WD values more stochastic disturbance is required, or a recruitment advantage for low WD trees, which is not currently included in the model.
To determine the importance of the WD shift for mortality, we calculated the mortality rate produced from the model output of the final simulation year for one run using either the existing WD value, or WD reduced by 0.09 g cm −3 .This showed that the increase in WD caused the stem mortality rate to be reduced by c. 12% of that simulated without the WD shift (Figure S16), and can explain c. one third of the total underestimation.In addition to the increase in WD during the simulations, the somewhat Correlations use the mean of the last 100 years of each run.
higher simulated than observed growth rates (Figure 6) may be an additional cause of low simulated stem mortality.
During model development, we initially used a P m equation parameterized with a larger plot dataset from the whole Amazon basin.However, there is a 2-fold difference in stem mortality rates across the basin (Quesada et al., 2012).When we used the basin wide parameterization, the underestimation of mortality rates was more pronounced (on average 0.65% yr −1 across the four plots).The spatial variation in mortality is thought to be linked to differences in species composition that are related to soil structure (Quesada et al., 2012).Thus, by including WD in the estimation of P m , some of the spatial variation in mortality that is related to species composition should be accounted for.Given that we produce a higher mortality rate when parameterizing the equation with only Western Amazon plot data, despite taking into account the species composition of the plots, this suggests that other factors (possibly soils) beyond WD that vary spatially are also important for determining tree mortality, and are areas for future research and model development.
An interesting emergent outcome of the simulations is the linkage between stem mortality and recruitment rates (Figure S9; Table 3).These processes are modeled independently, with recruitment driven by subplot-level LAI, and mortality linked to tree-level growth, size and WD.As for mortality, the recruitment scheme also includes high stochasticity.Given that the simulations run to a steady-state stem number, it is necessary that the demographic rates are similar over the long-term.The potential mechanisms that lead to linkages between mortality and recruitment rates are summarized in Figure 8. Stem mortality itself influences subplot LAI directly by removing trees (and their crowns) and indirectly by influencing tree lifespan and hence how large crowns can grow.These effects may be the reason why WD had an impact on mean plot LAI in the traits experiment (Table 4).Interannual variation in LAI is mostly driven by climatic events, with years with a particularly strong dry season or low incoming radiation resulting in a reduction in LAI due to low NPP, and therefore an increase in recruitment.This low NPP also influences growth rates, and hence more trees die when the climate is less favorable.
There was surprisingly little influence of leaf traits on GPP despite the explicit linkages in the model between LMA, leaf FIGURE 8 | Conceptual diagram of model inputs (filled ovals) and model processes that affect mortality and recruitment rates.1. Climate impacts NPP, with lower NPP under more seasonal climates.NPP then impacts LAI by determining the amount of carbon available for leaf production.2. NPP impacts diameter growth rate by determining the amount of carbon available for wood production.Traits, predominately WD, also impact diameter growth as the biomass-diameter allometry is dependent on WD. 3. Mortality rates are directly influenced by WD and by diameter growth rates.4. Diameter growth directly influences tree size, and mortality indirectly influences tree size via tree longevity (trees with a longer lifespan can reach larger sizes) and directly by removing trees. 5.The stand structure, particularly the presence of large crowns, determines the LAI.The recruitment rates are directly driven by plot LAI.
nutrients and the photosynthetic capacity.Comparing results from the two experiments (high vs. low WD and "fast" vs. "slow" leaf traits), the different WD distributions (with very little differences in other traits) has a stronger influence on model output than the different leaf trait distributions (which still have some difference in WD).The effects seen in the leaf traits experiment may therefore be more to do with the WD differences.WD showed large shifts over all model simulations (as discussed above), but there were only small, and less consistent, shifts in leaf traits (Figure S12) with seemingly little competitive advantage for particular leaf traits in the model.The lack of large shifts in leaf nutrients and influence of leaf traits on GPP suggests that trade-offs are more balanced for these traits.High leaf nutrient concentrations confer high photosynthetic rates but also high leaf respiration rates, while low LMA values confer high leaf area (for a given leaf biomass) but shorter leaf lifespans.Where shifts did occur they largely favored low leaf nutrient concentrations which may be due to selection for lower respiration rates, or to the weak negative association between WD and leaf nutrients.Whilst there is only scant data on shifts in plant traits, a recent study of five old-growth Neotropical forests over 10 years observed increases in WD with no change in leaf nutrients (van der Sande et al., 2016).However, no such shift in WD has yet been detected in much larger long-term plot networks (Lewis et al., 2009;Brienen et al., 2015;Esquivel-Muelbert et al., 2019).

Importance of Climate
Our results suggest a high sensitivity of forest productivity to climate.However, a limitation of the model is that we do not account for variation in plant hydraulic traits.In the model, the impact of waters stress is mediated by a simple link between soil moisture and stomatal conductance.Whilst this approach is used in many vegetation and land surface models, its use has been called into question as it performs poorly compared to other formulations (e.g., Egea et al., 2011;Joetzjer et al., 2014).The use of hydraulic traits and simulation of the full soil-plant-air continuum could provide better estimation of the sensitivity of photosynthesis to water stress (e.g., Christoffersen et al., 2016), although care is required in the parameterisation of soil hydraulics (Verhoef and Egea, 2014).The inclusion of plant hydraulic traits would be especially useful in models aiming to incorporate the wide variety of strategies in tropical forests, such as TFS.A variety of water use strategies have been shown (Klein, 2014), and plant hydraulic traits show correlations (Bartlett et al., 2016;Meinzer et al., 2016) potentially allowing inclusion of hydraulic trait strategies from measurements of few parameters or from traits already typically included (Markesteijn et al., 2011;Christoffersen et al., 2016).As well as improving representation of response to seasonal water stress, adding such traits could improve prediction of drought related mortality (Fisher et al., 2018), which has been shown to be linked to hydraulic failure (Rowland et al., 2015).Future shifts in precipitation could cause changes in species composition with consequences for forest dynamics, further supporting the need for the incorporation plant hydraulics into individual-based models that can allow for changing functional composition.
Empirical studies on tropical forests have also attempted to determine the roles of plant traits and/or climate on forest structure and functioning, particularly focusing on above ground carbon.For example, Poorter et al. (2015) found that rainfall was the strongest driver of AGB along a broad rainfall gradient, when considering effects of rainfall, average stem size, diversity and soil fertility.Poorter et al. (2017) assessed the impacts of climate, community mean traits, soil properties and disturbance (logging) impacts on biomass stocks and dynamics, and report findings broadly supported by those presented here.They found that AGB was controlled by both seasonal water stress and WD (c.f.Table 4), and that woody productivity (referred to as biomass growth of survivors) was controlled by mean annual precipitation and WD.In contrast with their findings, here leaf traits had an impact on AGB equal to that of WD (Table 4), while LMA was not important for AGB in the analysis of Poorter et al. (2017).In a study of an Andean altitudinal gradient using a simplified and non-demographic version of TFS v.1, Fyllas et al. (2017) found that functional trait and climate changes along the transect were the key drivers of changes in GPP and NPP with elevation, with a stronger role for traits on NPP than GPP.The results here are similar, with climate the key driver of GPP and NPP, and small role for leaf traits in controlling NPP.However, the main aspect of the climate differs in these two cases, with solar radiation important in the Andean elevation gradient (which is moist throughout) and precipitation important in the lowland Western Amazon.
Through a very detailed evaluation of the model, we have highlighted areas of poor model performance (organ level NPP and R, stem density, stem mortality) and suggested areas of future model development.While these areas of poor performance exist, many of the higher order variables were simulated well (GPP, NPP, Ra, AGB, BA).This gives us more confidence in the experimental results with regard to these well simulated variables, as does the congruence with the empirical results of Poorter et al. (2017).However, the results should be nonetheless be viewed with caution because of the uncertainty in the accuracy of internal mechanisms.

CONCLUSIONS
We analyzed carbon fluxes, forest structure and demographic rates using an individual-based, traits-continua model.The unique model structure allows us to explore linkages between climate, traits, carbon cycling, tree mortality and forest structure.Crucially, we include trait continua, which allow representation of diverse plant strategies.We found that between two contrasting sites, differences in productivity were driven by climate, rather than plant traits.Drivers of simulated forest structure were more complex, with a key role for WD mediated by its link to tree mortality, and surprisingly little effect of varying leaf traits.Our simulations showed consistent shifts in WD during simulations but not leaf traits, and highlighted the importance of including trait trade-offs.Future model development should focus on linking soil physical properties to mortality, improving allocation between organs and organ level respiration, differentially parameterizing the physiology of trees under different levels of shading, and adding hydraulic traits.This type of model that incorporates diverse tree strategies, detailed forest structure and realistic physiology is necessary if we are to be able to simulate tropical forest responses to global change scenarios.
authors contributed data and commented on or approved the manuscript.

FIGURE 1 |
FIGURE 1 | Schematic of the dynamic version of the TFS model.* denotes processes added to the new model version, ∧ denotes alterations to the original model.Thin arrows show where state variables are used in model processes or are used in the description of other state variables, thick arrows show where one model process informs a sequential process or the state variables.A more detailed schematic is included in the Figure S1.

FIGURE 2 |
FIGURE 2 | Mean monthly values of temperature ( • C) (A,B) and precipitation (mm) (C,D) for Allpahuayo and Tambopata taken from Sheffield et al. (2006) gridded climate data for the years 1983-2012 (blue, solid line: mean, dashed line: standard deviation) and from local data (red, points and error bars, del Aquila-Pasquel et al., 2014, Malhi et al., 2014 obtained with Web Plot Digitizer).

FIGURE 3 |
FIGURE 3 | Observed and simulated carbon fluxes for Tambopata and Allpahuayo plots for validation simulations.Total GPP (A), NPP (B) and Ra (C), leaf, wood and fine root NPP (D-F), and leaf, wood and fine root Ra (G-I).All units Mg C ha −1 yr −1 .Observations from Malhi et al. (2015) with error bars as standard error.

FIGURE 5 |
FIGURE 5 | Stem number (A-D, stems ha −1 ) and basal area (E-H, m 2 ha −1 ) from 500 years of simulation (blue-standard deviation range across all simulations), and data (brown-standard error range across all plots within the site cluster, red-standard error range of the simulated plot) for validation simulations.

FIGURE 6 |
FIGURE 6 | Distributions of stem diameter growth rates (cm yr −1 ) from observed data (red) and simulations (blue) for validation simulations.

FIGURE 7 |
FIGURE 7 | Stand development (a-d stem density, e-h basal area, i-l biomass) and carbon fluxes (m-p GPP, q-t NPP, u-x Ra) under different combinations of input trait distributions and climates from experimental simulations.Column 1 and 3-ALP climate, column 2 and 4-TAM climate, column 1 and 2-"fast" (blue) and "slow" (green) traits distributions with split based on leaf N, column 3 and 4-high (magenta) and low (cyan) WD distributions.Polygons show standard deviation ranges from 10 model runs.

TABLE 2 |
Observed and simulated carbon fluxes and forest structure metrics for validation simulations.

TABLE 4 |
Proportion of variance explained (ω 2 ) by traits distribution, climate, and their interaction for various simulated carbon cycling, forest structure and demography variables from experiment simulations.

TABLE 5 |
Correlation matrix of simulated plot-level variables from each run of the experimental simulations (10 runs of 8 climate x trait combinations, 80 simulations).