Methods

Methods include the calibration of the forest structure and th leaf phenology modules with inventories data and leaf litterfall observation from Tapajos and Paracou (see inventories and fluxes), exploring the variations among repetitions of simulations with similar parameters, and the methods used for simulation outputs evaluation.

Structure calibration

We calibrated the three forest structure parameters (m,  aCR  and bCR) for each site. aCR  and bCR are not independent, and we used the TALLO global database of crown radius (CR) and diameter (dbh) measurements (Jucker et al., 2022) to infer their relationship. To do so, we restricted the TALLO database to observations located within 10 km around sites from which we generated a thousand pairs of (aCR, bCR) values. Each pair of values was determined by randomly drawing 10 individuals per 10-cm diameter class to generate a size-balanced dataset to which the following model was fitted:  log(CR) ~ N[aCR +bCRlog(dbh), 2]. This resulted in the following linear relationship between the two parameters:  bCR =-0.39+0.59aCR +bCR , with bCR the error around the relation. This relationship constrained the exploration of the three-dimensional parameter space, so we only had to calibrate aCR, bCR , and m. Based on preliminary exploratory analyses with the previous version of TROLL, we defined the range of calibration for each parameter and site as follows:  aCR varied from 1.60 to 2.00 in Paracou and from 2.3 to 2.7 in Tapajos with a step of 0.05, bCR from -0.30 to 0.10 in both sites with a step of 0.05, and m from 0.030 to 0.050 in both sites with a step of 0.0025. This resulted in 9 aCR 5 bCR9 m2 site=810 triplets of parameter values. 

For each set of three parameter values, we performed a 600-year simulation from bare ground over a 4-ha area. Simulations were run with an external seed rain uniformly distributed across species, so that the simulated community structure is an emergent property resulting from the community assembly mechanisms embedded in the model. As succession unfolds and the number of mature trees increases in the simulation, internal seed production increases according to the assumed relationships between individual size and fecundity. An alternative to uniform seed rain across species would be to prescribe non-uniform seed rain based on species’ regional abundances. This approach would tend to make the simulated species abundances more closely resemble the observed regional abundances. In contrast, uniform seed rain as simulated here, biases the simulated abundances towards evenness across species, and differences in simulated abundances reflect differences in demographic performance controlled by the model trait-based parameterization rather than prescribed differences in the seed rain. Each simulation was forced each year by randomly drawing a year among the ten years of climatic data. In doing so, we avoided applying a periodic climatic forcing or any potential trend linked to global warming.

To evaluate the forest structure simulated with each triplet of parameter values, we compared simulated to observed total aboveground biomass (AGBtot, Mg.ha-1), total number of stems (Ntot, ha-1), and tree abundances per 5-cm diameter class (Ni, ha-1 for dbh class i) at the end of the 600-year regeneration. The Paracou reference dataset was a 2015 inventory of trees with dbh >10 cm in six 6-ha plots (Derroire et al., 2023). The Tapajos reference dataset was a 1999 inventory of trees with dbh > 10 cm in 19.75 ha along four 1-km transects (Rice et al., 2004). At both sites, we calculated the relative root mean squared error. We extracted the simulation with the lowest RRMSEP at each site and used the corresponding values for m,  aCR  and bCR in all subsequent simulations.

Phenology calibration

The new leaf phenology modules lay on a shedding threshold potential:

\[ \Psi_{T,o}=min(a \times \Psi_{TLP}, -0.01\times h-b) \]

with \(\Psi_{T,o}\) the leaf shedding potential in MPa, \(\Psi_{TLP}\) the leaf water potential at turgor loss point in MPa, and \(h\) the tree height in metre. This equations thus introduces to new parameters which are \(a\) the leaf resistance parameter and \(b\) the tree height susceptibility parameter. Once the leaf shedding potential is reached tree losses their leaf quicker following \(delta\) the base variation in residence time multiplying factor for drought shedding. While this scheme is process-based uncertainty relied on the values of \(a\) , \(b\) and \(delta\) that we calibrated. We all the combinations of the following values:

  • a0: [0.01, 0.025, 0.05, 0.075, 0.1, 0.2, 0.3, 0.4, 0.5]
  • b0: [0.05, 0.01, 0.015, 0.02, 0.04, 0.06, 0.08, 0.10, 0.12, 0.14, 0.16, 0.18, 0.20]
  • delta: [0.1, 0.2, 0.3, 0.4, 0.5]

We first ran a spin-up simulation from bareground during 600 years to establish a mature forest with \(a=0.2\), \(b=0.02\) and \(delta=0.2\). We then ran 10 years of simulation with each combinations of parameters that we burnt, to allow the forest and leaf dynamics to adjust to new parameters, and we kept a second run of 10 years for calibration. For each simulation and couples of parameters, we built simulated observation of litterfall based on the sampling date of the observed litterfall in Paracou and Tapajos (with corresponding varying time lapse). We further computed for each simulation the same two metrics built with observed litterfall:

  • Day of the litterfall pike as the day of the maximum annual value
  • The litterfall rate is the ratio between the pike value (average of the annual maximum value and the previous and following values) divided by the base flow (average values between January and April)

We computed for each simulation the Root Mean Square Error of Prediction (RMSEP) of the euclidean distance between observed and simulated pikes and rates:

\[ RMSEP = \sqrt{| (ratio_o-ratio_s)^2+(pike_o-pike_s)^2 |} \]

We finally selected the set of parameters for each site with the lowest RMSEP to be used in TROLL evaluation.

Variations approach

To assess TROLL simulations variations due to stochasticity we ran ten repeats of 10-years simulations at Paracou and Tapajos with best parameters for a single spin-up of 600-years from bareground for each site. As we were mainly interested in daily fluxes that adjust to change in model parameters over a few years, we used only variations between 10 years to save ressources and computation time. But for a better evaluation of TROLL stochasticity for forest structure and composition we should used repeats of 600 years spin-up simulations from bareground.

Evaluation methods

We run an extra 10-years simulations (run3) with best parameters after calibration for TROLL evaluation. We evaluated fluxes of carbon (growth primary productivity GPP), water (evapotranspiration ET and relative soil water content RSWC) and leaf (Leaf Area Index LAI), forest structure (basal area BA & abundance) and forest composition of species (rank abundance curves) and functional traits (quantiles). For fluxes, we investigated daily variations, 15-days variations, daily variations against known drivers, daily observed versus predicted, and computed several evaluation metrics (R squared of the linear regressions of null intercept R2, Pearson correlation coefficient CC, root mean square error of prediction RMSEP, standard deviation of the error of prediction SD and their relative counterparts RRMSEP and RSD). We evaluated fluxes when available against eddy-flux measures from towers, ground observations (litter traps and terrestrial lidar) or satellite observations.

Finally, to quantify the envelopes of stochastic simulation outputs, we ran ten replicates of 600-year simulations starting from bare ground with the six calibrated parameter values.