Modeling the Effects of Global Change on Ecosystem Processes in a Tropical Rainforest

: Research Highlights: Ongoing land ‐ use change and climate change in wet tropical forests can potentially drive shifts in tree species composition, representing a change in individual species within a functional group, tropical evergreen trees. The impacts on the global carbon cycle are potentially large, but unclear. We explored the differential effects of species within this functional group, in comparison with the effects of climate change, using the Century model as a research tool. Simulating effects of individual tree species on biome ‐ level biogeochemical cycles constituted a novel application for Century. Background and Objectives: A unique, long ‐ term, replicated field experiment containing five evergreen tree species in monodominant stands under similar environmental conditions in a Costa Rican wet forest provided data for model evaluation. Our objectives were to gain insights about this forest’s biogeochemical cycles and effects of tree species within this functional group, in comparison with climate change. Materials and Methods: We calibrated Century, using long ‐ term meteorological, soil, and plant data from the field ‐ based experiment. In modeling experiments, we evaluated effects on forest biogeochemistry of eight plant traits that were both observed and modeled. Climate ‐ change simulation experiments represented two climate ‐ change aspects observed in this region. Results: Model calibration revealed that unmodeled soil processes would be required to sustain observed P budgets. In species ‐ traits experiments, three separate plant traits (leaf death rate, leaf C:N, and allocation to fine roots) resulted in modeled biomass C stock changes of >50%, compared with a maximum 21% change in the climate ‐ change experiments. Conclusions: Modeled ecosystem properties and processes in Century were sensitive to changes in plant traits and nutrient limitations to productivity. Realistic model output was attainable for some species, but unusual plant traits thwarted predictions for one species. Including more plant traits and soil processes could increase realism, but less ‐ complex models provide an accessible means for exploring plant ‐ soil ‐ atmosphere interactions.


Introduction
Interactions between plants and soil are fundamental to terrestrial ecosystem processes, which are strongly influenced by changes in climate and land use [1]. These changes influence biological activity, which in turn alters C cycling, with feedbacks to climate change. As engines of Earth systems, tropical systems are particularly important because they regulate global cycles of C and water. The tropics also house the majority of the world's population and highest biological diversity on Earth, with many species interacting in complicated and multifaceted ways. Climates without modern precedents such as warmer air temperatures are predicted to emerge earliest in the tropics, creating a sense of urgency for understanding climate-change effects on tropical ecosystems that are expected to have large impacts on ecological and social systems [2]. The tropics have also been undergoing land-use change, with timber harvesting occurring in at least 20%, nearly 4 million km 2 , of this biome in 2000-2005 [3]. Gross C emissions from tropical deforestation were reported to be 2.2 Pg C year −1 for the period 2000 to 2010 [4,5]. This loss represents a large potential for replacement of that forest C through tree regrowth and replenishment of soil C stocks. Considering the expansion of forest plantations [6], often containing a single planted species, a better understanding of the biogeochemistry of shifts in individual species will inform modeling approaches and management decisions.
Process-based models, which are based on vegetation characteristics and their relationships to biogeochemistry and hydrology, can provide a research tool for testing and understanding the mechanisms by which plants, soil, and climate influence ecosystem functioning, and for predicting future scenarios under global change [1,[7][8][9][10][11]. Modeling approaches differ in their complexity and scale of resolution, ranging from global context models that are focused at the 2° grid cell scale to the plot or micro-site scale [1]. All models have inherent trade-offs between complexity, realism, interpretability of output, and user accessibility [12]. Thus, the choice of model depends on the purpose of modeling, questions to be addressed, and constraints in model use, including data available for calibration and evaluation.
Our purpose was to gain insight into effects on 'within-biome' properties and processes in tropical rainforests that could influence biome-level biogeochemical cycles, potentially with consequences for global cycles. More specifically, we sought to capture effects at the micro-scale of individual species within the functional group of tropical evergreen trees, in comparison with the effects of two regional climate-change aspects, and simulate their effects on biogeochemical cycles at the biome scale. In this case, if a single type of model can be used flexibly, it can bridge two very different yet potentially linked system conceptualizations, one at the global/macro/extracted scale and the other at the micro/realistic scale, and thereby provide insights across scales. Thus, our other purpose was to evaluate the Century model in terms of this flexibility. Century is a process-based, biome-scale ecosystem model because its linking of C, N, and phosphorus (P) cycling allows for analysis at the regional level [13][14][15], and thus identification of strong feedbacks to global cycles. At the microscale, the forest version of Century contains over 50 tree-trait parameters, which pertain mostly to biogeochemistry rather than physiology, and to which our field-based observations had been purposefully aligned. Century has not been used previously to compare the effects of plant species within a functional group and relate them to biome-level biogeochemistry.
The tropics contain at least 40,000 species of trees, compared with 124 in continental Europe [16]. Thus, modeling the effects of individual tropical tree species would be daunting without the organizing principle of plant functional traits (PFTs). The concept is that species composition could influence ecosystem dynamics because individual species can differ in their acquisition, allocation, use, and storage of resources. Differences among tree species in traits that influence the quantity, chemistry, and location of detrital production [17,18] could drive changes in C and nutrient cycling, with feedbacks to ecosystems, and ultimately to global cycles, given the extent of changes in species composition resulting from land-use changes such as reforestation following deforestation. Our fieldbased studies indicated that species within the tropical evergreen tree functional group differed in multiple traits, and in their effects on ecosystem processes. For example, C and N cycling differed two-to eight-fold among five species of trees planted in experimental plantations in a single site in wet tropical forest [18][19][20]. Reforestation in the humid tropics in 1990-1997 was estimated at 10,600 km 2 /yr based on satellite imagery, with 60% of it attributed to plantations [21]. A change in species composition and decreased diversity usually accompanies plantation establishment, which may also occur during natural regeneration if seed sources of native species are lacking. The consequences for rainforest properties and processes of shifts in tree species composition and diversity at a regional scale are not well understood, however.
Despite differences among process-based models, they share general biophysical concepts and incorporate basic biological constraints, trade-offs, and physical processes that occur across all plant-soil-atmosphere systems [22]. As such, the model equations or parameter values may be viewed as explicitly stated hypotheses to be tested across ecosystems, i.e., parameters are regarded as variables, as a way to evaluate our general concepts of biogeochemical and hydrological cycling. Using this principle, rather than simply adjusting model parameters for each site to attain expected outputs, one can elucidate processes that need clarification through field-based studies. The process of calibrating a process-based model for a specific site can reveal gaps in our understanding of water and nutrient cycling or the existence of previously unrecognized pathways. For example, in calibrating the Century model for Hawaiian ecosystems undergoing long-term soil development, researchers found that modeled P budgets could not be balanced without allochthonous inputs of P. This led to the search for the 'missing' P and discovery that P inputs in dust from Asia contributed substantively to ecosystem-level P budgets [23].
For this modeling study, our objectives were to evaluate the effects on tropical rainforest biogeochemistry of: (1) two aspects of climate change that have been already occurring in the region and (2) variation in plant traits within a single functional group, tropical evergreen trees. We examined whether input parameters regarding plant traits in Century were sufficient to capture realistic effects of individual species on ecosystem dynamics. We integrated data from a unique, longterm field study in a replicated experimental setting to examine whether richer observed data about the vegetation can inform regional-level models about ecosystem dynamics. Field-based measurements of various plant traits and ecosystem properties and processes in five evergreen tree species were conducted at La Selva Biological Station in the Caribbean lowlands of Costa Rica, in a site in which climate, soil type, and land-use history are similar among replicated mono-dominant plantations, described by Fisher [24]. We first calibrated Century for this site, and then used the model in an experimental mode to address three topics and questions: 1. Model adaptation for the sub-biome. What insights about the ecosystem development and properties can be gained from calibrating Century for this site? We predicted that La Selva P input parameters would not need adjustment because wind patterns in this region are not likely to bring P dust inputs. We predicted that because La Selva has such high rates of N cycling [19] relative to Puerto Rican rainforests [25], N fixation parameters would need to be adjusted higher than estimates for other tropical rainforests. We expected that high rainfall at this site could limit aboveground NPP (ANPP) via limitations of light on cloudy days. 2. Sensitivity to climate variables. Which aspects of climate change would have the most effect on ecosystem processes? Tropical wet forests are already experiencing measurable climate change, manifested in Central America as increased minimum air temperatures and drier dry seasons [26,27]. We predicted that higher night-time temperature would have a greater effect in this tropical wet forest, given the negative effects on aboveground ANPP found by Clark et al. [28]. 3. Sensitivity to plant traits and model realism. Is this biome-based simulation model sensitive to variation in traits within a functional group of species? If so, which modeled traits have the most effect on ecosystems, and what is the extent of their effects? We were most interested in insights about the biological underpinnings of the ecosystem responses to functional plant traits (PFTs), specifically for traits associated with C and N dynamics. Previous field-based studies indicated that eight traits that are already included in Century had differed significantly among five experimental tree species in the same functional group; the species also differed in their effects on ecosystem properties and processes [18,19,29,30]. We also addressed the question of whether Century would generate realistic ecosystem properties under the five experimental tree species, and what types of model adjustments would be required to achieve this.
The model is composed of a soil organic matter/decomposition submodel, a water budget model, a grassland/crop submodel, a forest production submodel, and management and events scheduling functions [14]. It computes the flow of C, N, P, and sulfur through the model's compartments. Century does not explicitly incorporate the effects of plant species per se, but the input files define numerous variables regarding the biogeochemical and phenological characteristics that vary among plant species. Century contains a file that allows the user to specify sequences of plants and management, including grazing, crop planting and harvest, fertilizer application, tree removal, and burning. For the studies described here, we used the forest submodel of Century, version 4.5 (https://www.nrel.colostate.edu/projects/century/).

Data Available
La Selva, in the Atlantic lowlands of Costa Rica (10° 26 N, 83° 59' W), had a mean annual precipitation (MAP) of 396 cm and mean annual temperature of 25.8 °C prior to 1992 [40]. Input data for monthly means of precipitation ) and temperature  were obtained from La Selva's meteorological station (www.ots.ac.cr/meteoro/default.php?pestacion=2), for which MAP was 420 cm (Table S1).
Input variables regarding vegetation, e.g., tissue chemistry and leaf death rate (Table S2), were based on measurements in the five experimental species from a 12-ha site within La Selva. In this randomized block design that contained several tree species, a single species of evergreen tree had been planted in one 50 × 50 m plot within each of four blocks in 1988, making it the oldest replicated experiment containing native tropical tree species. Previous land-use history was the same in all plots. The old-growth forest had been cleared in 1955, the slash burned, and pasture of Panicum maximum L. and Melinis minutiflora Pal. established in 1956. Grazing had been continuous until abandonment in 1987. The native species planted in this study included: Hieronyma alchorneoides Allemao; Pentaclethra macroloba (Willd.) Ktze, a nodulated legume; Virola koschnyi Warb.; and Vochysia guatemalensis J.D. Smith. Their corresponding 4-letter acronyms are: HIAL, PEMA, PIPA, VIKO, and VOGU. The single exotic species was Pinus patula ssp. tecunumanii (Eguiluz and J.P. Perry), hereafter referred to as Pinus patula. Fisher [24] provided a complete description of the experiment. To provide a reference for the plantations, a fifth block was established in the remnant of old-growth forest on the same landform, soil consociation, and elevation, in 2003. We quantified approximately 270 tree species traits and ecosystem processes and properties in these experimental plantations [18,20,29,30,41]. Of these measured variables, 59 are parameters in the Century model (Table S2). The five experimental species encompassed a wide range of values for each trait, measured at ages 15-17 years, with all trees growing under the same climate, soil, and land-use history (Table S2).
The soils at the experimental site are thought to be derived from andesitic lava flows estimated to be ~1.5 million years old [42]. These residual soils have been classified as Mixed Haplic Haploperox [43]. Soil C and N data were available for the plantations and old-growth forest [30]. Old-growth forest NPP and biomass C and N data were estimated from the literature [44][45][46].

Model Parameterization and Calibration
We parameterized the model for La Selva's latitude, longitude, and climate, with the latter used in the stochastic weather mode to generate inter-annual variability in weather. To parameterize for old-growth forest in this site, we used trait data for Pentaclethra macroloba (Table S2) because it is the most dominant tree species at La Selva (Clark and Clark 2000). We ran the parameterized model for 1.5 million years, the age of the parent material. We then estimated model parameters that allowed it to predict output variables that corresponded with observed data for four stocks: aboveground biomass C and soil C, N, and organic P. We began with default model parameters used previously for a Puerto Rican rainforest. Simulated ecosystem development was rapid, but accrual of wood, hence forest development at La Selva was not achieved with these default N fixation rates and parent material weathering rates. We adjusted inputs of N and P incrementally until simulated values for these four output variables matched observed values ( Figure 1a, Table S3). Attaining measured wood production in La Selva simulations required an increase in N fixation and annual additions of 0.23 N and 0.003 P g/m 2 . We refer to the calibrated 1.5-million-year model run as the 'baseline' model run. For modeling the experimental plantations, we appended an input file to the baseline run to simulate land-use changes prior to initiation of the experimental site, e.g., deforestation, implementation of pasture and its abandonment after ~35 years, and plantation establishment and management factors (Table S4). We adjusted N and P inputs and traits of the juvenile default tree species (Pentaclethra) in the plantations as in the calibration for the baseline run.  Table  S2), modeled both before deforestation and during the experimental plantation period.

Modeling Experiments
In two separate sets of modeling experiments designed for the wet tropical forest biome, we evaluated the effects on ecosystem processes of: (1) two aspects of climate change that now occur in this region and (2) eight tropical tree traits that have the potential to influence biogeochemical cycles. For each experiment, we conducted a single model run at each of five values for the variable and recorded the output variables for multiple ecosystem properties and processes at the end of each 100yr run, i.e., steady-state.
Without further calibration, we conducted experiments by appending simulations to the baseline run. For a given experiment, we systematically varied a single input variable over the course of five model runs. We ran each experiment for 100 years, but steady-state was reached by ~20 years. The output (or response) variables evaluated were C and N stocks in total (above-and below-ground) biomass and soil and total NPP. In the climate-change experiments, we also assessed streamflow and N leaching. Where appropriate, other response variables were also tested, including potential evapotranspiration (PET), evapotranspiration (ET) that would occur if water were not limiting, and actual evaporation and transpiration, the quantity of water actually removed from the forest via evaporation and transpiration (AET), and net N and P mineralization.
Climate-change experiments were designed to test the effects of two major changes in climatechange scenarios that have been observed in Central America, particularly at this site: (a) Increasing minimum air temperatures (Tmin). From the current mean Tmin of 21.47 °C at La Selva, Tmin was increased incrementally over the series of model runs to 25.77 °C (Table S5). This corresponded to a 2.1 °C change in the daily average T. The difference between Tmax and Tmin across all months of the experimental runs ranged from 7.2-8.9 C, whereas the historical range within the year is 8.4-10.3 °C for this site (http://www.ots.ac.cr/meteoro/default.php?pestacion=2). (b) Decreasing rainfall during the drier months. Annual precipitation was decreased from 420 to 347 cm (Table S6), with all decreases occurring from January to May, the drier period at La Selva.
Tree-species-trait experiments were conducted to evaluate the relative effects of eight traits on ecosystems. The selection of tree species variables to measure and model was based on knowledge of the trait's capacity to influence resource (C and nutrient) capture and cycling, as well as Century's capacity to incorporate these traits and simulate their consequences for ecosystem processes. The five experimental values for each of the eight traits spanned a realistic range based on literature values (where data were available), with the middle value representing 'La Selva old-growth forest' and two values spaced to represent potential variation above and below La Selva's value. 'Production potential' in Century is a conceptual trait, scaled from 0 to 1.0; it represents the potential for gross forest production. For this experiment, values of 0.3-0.7 were estimated based on previous use of this model for tropical forests and observed data on production in this site. For several traits, the values (in parentheses) were based on our field-based experiment: 'Allocation to fine roots' (0.1-0.4); 'Fineroot death rate' (0.02-0.18, with the upper range extended to encompass a broader range of possibilities); 'Leaf C:N' (12)(13)(14)(15)(16)(17)(18)(19)(20)(21)(22)(23)(24)(25)(26)(27)(28); and 'Fine-root C:N' (23-42.5) [18,19]. 'Leaf death rate' was based on values from Reich et al. [47]. 'Leaf lignin' (0.1-0.4) was based on data from Constantinides and Fownes [48] and Aerts [49] and 'Fine-root lignin' (0.15-0.35) on Russell et al. [30] and Cusack et al. [50].

Evaluation
To evaluate whether the simulation of land-use changes had set ecosystem properties on a trajectory that corresponded well with properties at the onset of the experimental plantations, we compared simulated biomass and soil C and N stocks (Figure 1b) with published values for this site [18,24]. Observed data for biomass P stocks were not available. Although true validation was not possible, we evaluated the performance of the model in simulating the effects of individual tree species. Within the input data file for tree parameters, we created a unique subfile for each of the five experimental species, using measured trait values (Table S2). Without any further calibration, the model file that simulated land-use changes, including plantation installment and growth, was appended to the calibrated baseline (old-growth forest) simulation. We developed separate programs for each experimental species simply by identifying the species-specific tree sub-file described above. We compared simulated output parameters and values measured in four blocks over 2003-2005 in the experimental plantations for one process (NPP , Table S7) and four stocks (C and N of aboveground biomass and soil, 0-30 cm [18].
For output parameters for individual species that fell outside of the 95% confidence interval of measured values across the four blocks in the field experiment, we then experimented with adjusting both plant-and site-level input parameters to determine whether it was possible to attain more realistic output values. Plant parameters that we adjusted included

Climate-Change Experiments
With Tmin increasing by 4 °C from the current 21.47 °C, simulated biomass declined by 21%, from 9946 to 7868 g C m −2 and NPP declined by 12%, from 1497 to 1321 g C m −2 year −1 over the series of five experimental model runs, with each run for 100 years (Figure 2a). Potential ET declined monotonically from 212 to 165 cm yr −1 , by 22% across this range of Tmin values (Figure 2b); AET decreased by 20%, from 200 to 159 cm yr −1 . Simulated soil C declined by 16%, from 8420 to 7058 g C m −2 , while soil N decreased by 19%, from 770 to 623 g N m −2 (Figure 2c). Streamflow increased by 18% and N leaching declined by 8% as Tmin increased (Figure 2d). Net N mineralization declined from 34.1 to 30.2 gN m −2 yr −1 and net P mineralization declined from 1.69 to 1.51 gP m −2 yr −1 as Tmin increased (Figure 2e).
Drier period in January-May. By decreasing drier-season rainfall, simulated annual rainfall declined 15% over the range of experimental runs. At the experimentally highest rainfall, i.e., 'normal' MAP of 420 cm, biomass was lowest (Figure 3a). It increased by 6% to peak at the intermediate rainfall level of 376 cm yr −1 , and then biomass declined by 4% as rainfall decreased to the lowest level. Soil C and N both declined monotonically by 16% across the range of increasing rainfall (Figure 3b). The experimental rainfall reduction was during the drier part of the year; nonetheless, streamflow decreased by 16% as rainfall declined by 17% (Figure 3c). N leaching declined by 16%, with leaching increasing relatively more at higher rainfall. Nutrient availability trends were similar to those of biomass in that net N mineralization peaked at intermediate rainfall levels, and then declined by 3% as rainfall increased (Figure 3d). P mineralization declined 4% as rainfall increased, but this process is less influenced by biological controls in Century, and was relatively constant under drier conditions.

Tree-Species-Trait Experiments
An increase in the trait of production potential across the simulated range resulted in an increase in C stocks, 40% in aboveground biomass and 6% in soil; both peaked at the second-highest level and then declined by 5% and 1% respectively (Figure 4a). Aboveground NPP increased by 19% ( Figure  5a). Increasing allocation to fine roots led to a 56% decline in aboveground biomass, an 18% increase in SOC (Figure 4b), and a 34% decline ANPP (Figure 5b). With increasing leaf death rate, aboveground biomass declined monotonically by 75% over the simulated range whereas SOC increased by 7% (Figure 4c); ANPP increased by 34% (Figure 5c). With increasing fine-root death rate, all parameters increased slightly, by 7%, 3%, and 2% in aboveground biomass, SOC, and ANPP respectively (Figures 4d and 5d).
Regarding tissue chemistry traits, with increasing leaf C:N aboveground biomass increased by 72%, SOC by 18% (Figure 4e), and ANPP by 34% (Figure 5e). Both SOC and ANPP peaked at intermediate values. In contrast, with increasing fine-root C:N aboveground biomass and SOC increased by <1% (Figure 4f) and ANPP by 3% (Figure 5f). With increasing leaf lignin, aboveground biomass declined by 12% while SOC increased by 8% (Figure 4g); ANPP decreased by 10% ( Figure  5g). Fine-root lignin had a trend similar to leaf lignin, but the effect was smaller, with all three output variables varying by only 3% across the range of simulated values (Figures 4h and 5h).

Modeled Tree Species
Among the 25 simulated output values for the modeled tree species (five species × five variables), 14 were within the 95% confidence interval (CI) for measured values across the four blocks in the field-based experiment ( Figure 6). For Pentaclethra, four of the five output variables were predicted well, whereas for Hieronyma and Vochysia only two and one values respectively fell within the 95% CI. Aboveground biomass N stocks were well predicted for four of the five species and soil N stocks for three species. Biomass C stocks and ANPP were well predicted for only two of the species. Soil C stocks were the least well-predicted output parameter.
For Pentaclethra and Hieronyma, relatively little adjustment of input parameters was required to attain output values within the 95% CI, and further fine-tuning brought output values closer to the measured means, with the exception of soil C in Pentaclethra ( Figure 6, Table S8). The main adjustments were for biomass-LAI and allocation to wood and fine roots. Similarly, Pinus reached realistic output values with fairly simple adjustments of increasing the wood death rate, and decreasing the leaf death rate and passive SOM decomposition. For Vochysia however, there was no combination of input parameter adjustments that yielded realistic results for all of the simulated ecosystem properties.  [18,30]. An asterisk denotes that uncalibrated simulated values fell outside of the 95% confidence interval of measured values.

Insights from Modeling Activities
Soil dynamics. We gained two insights about soil processes at La Selva as a result of three main model parameter adjustments that were required to achieve realistic soil C stocks during ecosystem development under natural forest and again under the plantation phase. The adjustments were: (1 and 2) additional N and P inputs were required to meet the nutrient demands for accruing C stocks in biomass and soil during the first 2000 years of soil development; and (3) the rate of passive SOM pool decomposition had to be increased during the plantations phase, above the rate during the mature forest phase, i.e., prior to deforestation. Two main questions arise. (1) What is the source of the extra N and P needed to balance the C and nutrient budgets? (2) What unaccounted-for processes are influencing SOM decomposition, such that soil C stocks were not initially well predicted for four of the five species in the field-based experiment?
Regarding additional N inputs required, it is possible that N fixation rates are higher at La Selva, compared with Puerto Rico. For P, one possibility is that the soil age is younger than documented [42], such that the soil was not as P-depleted as would be expected for residual soils derived from the 1-Ma-old parent material. This is consistent with a possibility posed by Porder et al. [51] that nearby volcanoes might have contributed younger, overlying material, lahar or ash deposits, from which soil was derived. Another possibility is that the process of erosion, which would rejuvenate soil with P inputs [51], is not included in the model. In addition, the ability to model P cycling has been hindered by methodological issues in measuring plant-available soil P [52].
Another explanation that addresses both questions is that sequential reduction and oxidation processes, which are not included in the model, may play an important role in C and nutrient dynamics in this site, as found studies in Puerto Rico [53]. This mechanism relates to the fact that tropical soils tend to be high in short-range order Fe and Al [43], which promote formation of organometal complexes and C sorption, and thus contribute to C accrual in tropical soils [54], while O2 availability may limit decomposition of mineral-associated C in soil micro-sites [55]. These processes could be influenced differentially by individual species. Changes in soil water dynamics during landuse change and root traits that differ among these five species [30,41] could affect these oxidationreduction processes, leading to concomitant changes in soil C and nutrient dynamics.
Another mechanism that can regulate soil nutrient availability involves the alteration of soil pH. This 'master' soil variable influences biogeochemical processes and soil properties in the variablecharge soils that occur at this site, including cation exchange capacity (CEC), base saturation, mineral speciation, and particle surface charge [30,56,57]. Previous work suggested that species traits that result in the release or consumption of H + , and thereby alter soil pH, could influence the liberation of occluded cations and P associated with colloids [57]. Nodulation that supports microbial N2 fixation in leguminous plants can acidify soil via ammonification and subsequent nitrification [58] and thus promote occlusion of P and cations. In contrast, Al-accumulating species alter the equilibrium between Al minerals and soluble Al, consuming H + [59], thereby increasing soil pH and liberating cations and P from occluded forms.
Water. Model output indicated that La Selva has high rates of potential evapotranspiration (PET), 46% of the mean annual precipitation (MAP) (Figure 2b) and this agreed with previous measurements in this forest [60]. At first this was puzzling however, because one might expect lower PET in a site that receives >400 cm MAP, given that cloudy and rainy conditions would tend to reduce evapotranspiration. Inspection of the diel rainfall data for La Selva however, revealed that a substantial proportion of rainfall occurs at night, which could explain why ET would not be very strongly diminished by high MAP.

Climate Change
Although the amplitude of climate change in the tropics may be smaller than at higher latitudes, small but fast changes are a cause for concern because tropical species may be adapted to narrower climate ranges [61][62][63]. Air temperature is expected to sustain the greatest changes in ecological systems due to global warming, and precipitation has also been identified as a key variable [2]. Warming in the tropics is of critical concern because tropical trees may already be pushed beyond their thermal thresholds for photosynthesis [64]. It is largely unknown whether tropical trees will acclimate to increased temperatures, i.e., whether tropical forests will remain a sink for CO2 or become a net source, further exacerbating climate change [65], and reducing biodiversity in tropical forests [61].
We focused on specific scenarios of climate change that have been observed in this site; trends of increasing night-time temperature (Tmin) and decreased drier-season rainfall are predicted for Central American forests in general [26,27]. Of these two scenarios, increasing Tmin had the greatest effect on simulated ecosystem processes (Figure 2). The increase in the minimum temperature generally reflects an increase in atmospheric moisture, which usually causes a reduction in PET. The reduction in the difference between the maximum and minimum temperature is also correlated with a reduction in PET [66]. Thus, Century uses the difference between the minimum and maximum temperature to calculate PET, which is high in deserts and low in tropical forests due to atmospheric water levels. Modeled NPP declined with increasing Tmin because the narrower difference between Tmax and Tmin led to a reduction in PET and AET (Figure 3b), which is essentially a surrogate for vapor pressure deficit (VPD) in Century. Consequently, this reduced productivity. Similarly, others [67] found that annual NPP correlated most with AET. With NPP declining as a result of lower AET, relatively less organic C was transferred to soil, thereby reducing N and P mineralization rates as Tmin increased. Also, as AET declined, streamflow increased, driving slight increases in P loss in streamflow. In their analysis of factors limiting productivity in Hawaiian forests, Raich et al. [23] found that modeled leaching losses of P increased over time of soil genesis, driving a concomitant decline in ANPP in P-limited forests. At La Selva, although leaching loss of P increased by only 1% as Tmin increased, that loss may have been enough to depress productivity, given the relatively high P limitation in this soil type.
In other long-term field studies in the old-growth forest at La Selva, woody biomass increments were negatively correlated with Tmin [28,68]. One possible mechanism by which this would occur is that leaf dark respiration (Rdark) increased with Tmin, reducing the storage of non-structural carbohydrates and thus net photosynthesis, as suggested by Asao and Ryan [69]. However, Rdark varies widely within and among species and PFTs [70] and some tropical species can acclimate to higher night-time temperatures [65]. Models that include this acclimation have very different simulated results of carbon fluxes between the vegetation and atmosphere [71]. The plant-level mechanisms by which carbohydrates are allocated, stored, and mobilized are poorly understood [72] however, making it is difficult to evaluate realism in model output at larger scales [71].
Under the simulated 'drier period in January-April' scenario (with lower MAP), productivity was not influenced as much as in the 'increasing Tmin' scenario because AET stayed fairly constant as precipitation was decreased. Nevertheless, biomass declined 19%, presumably due to the decreased rainfall. The impact on streamflow was small because the majority of streamflow occurs during the rainy period, in which rainfall was not altered in this scenario. With less soil water available during the drier period, N leaching also declined slightly with drier-season rainfall. On a global basis, there is concern that we underestimate the vulnerability of trees to extended and hotter droughts [73]. Given the ecological and social importance of trees, there is an urgent need for predicting thresholds for tree responses to more extreme drought.

Tree Species Effects
The tree-species-traits experiments demonstrated that this biome-level model was sensitive to variation in several of the modeled traits of tropical evergreen tree species. The traits of leaf death rate, leaf C:N, allocation to fine roots, and production potential all had a larger impact on biomass C stocks and ANPP in comparison with the climate-change scenario that had a larger impact, higher night-time temperature. The relationships among the traits were consistent with ecological concepts of leaf and stem economics and plant functional types [70,74,75]. The relationships between plant traits and ecosystem processes were also consistent with concepts based on observed data [17,49,[76][77][78][79]. These results have applications in selecting species for plantations. For example, species with high production potential tend to have high nutrient demands. If the demands exceed nutrient availability, as could be the case in highly weathered soils, biomass would be lower than the potential for that species, as in Figure 4a. In a degraded site, managers might select a species with high leaf litter production, i.e., high leaf death rate and leaf production, in order to promote an increase in SOM. Planting this species could result in lower wood production (and lower NPP) however, as a result of less N and P available for allocation to wood, given that leaves have a relatively high demand for N and P (low C:N and C:P) compared to wood. One might also select a species with high leaf and fine-root lignin content as a means to increase soil C stocks. Although the relationships between these traits and soil C stocks are consistent with studies that form the basis for processes modeled in Century [80], more recent work indicates that the residence time of lignin in soils can be much shorter [81], so there is uncertainty about this relationship.
Modeling of the five individual species revealed two main points. First, the species altered sitelevel processes differentially, thus requiring adjustment of site-level parameters in order to obtain realistic output for biomass C and soil C stocks under some of the tree species. In the model structure, these parameters, including soil pH, the decomposition rate of the passive SOM pool, and N fixation, are regarded as constant across a site. While it is known that trees species can influence soil properties and processes differentially and have a greater effect on soil than abiotic factors [17,[76][77][78], these plant-soil feedbacks are not be represented explicitly as plant input parameters (traits) in most models. Representation of these feedbacks from plants would advance the model's realism, especially for variable-charge soils.
This study highlighted that, despite our capacity to include dozens of plant traits in a model, an unusual plant trait not included in the model can thwart the capacity to predict the effects of individual species in the real world. Especially for one species in our study, Vochysia, it was impossible to adjust model input parameters in a combination that yielded realistic output for ecosystem properties and processes. We suspect that the reason that the modeled species could not accrue enough nutrients to attain the observed soil C and N stocks and ANPP was that this species' unusual means for acquiring nutrients, i.e., Al accumulation and or the capacity to alter soil pH and other properties, are not represented in the plant input parameters. Moreover, sequential reduction and oxidation processes that could be common in highly weathered tropical soils are not represented in the model, making it difficult to model plant nutrient availability, and hence NPP.

Trade-Offs in Modeling
Century 4.5 does not incorporate the intricacies of ecophysiology, such as the trade-offs and constraints on the role of woody stem production in energy capture. Thus, perhaps the most remarkable aspect of these simulation experiments is that differences among species in the biogeochemical traits alone provided as much realism as they did, especially considering that elemental cycles other than N and P were not simulated, nor aspects of plant biology other than ones related to biogeochemistry.
Representing plant processes in more detail is an active area of research that has yielded important advances in incorporating ecophysiology and plant demography, and thus greater realism, into Earth Systems Models (ESMs) such as the Community Land Model Version 4.5 (CLM), Functionally Assembled Terrestrial Ecosystem Simulator (FATES), and DayCENT. The added complexity raises interesting questions regarding specific gaps in ecological understanding, methods for interpreting model output, and quantifying uncertainty, which are also areas of collaborative research [10]. Increased complexity in vegetation dynamic models brings challenges in informatics and statistics with regard to input requirements and output dimensionality, however [8,9], although the development of model informatics systems is meeting these challenges [10,11]. These developments are important in advancing the predictive capacity of ESMs. If the trade-off for increased model complexity is reduced accessibility for non-modeling experts however, then simpler models with a lower predictive capacity will continue to provide non-specialists with a reasonable research tool for exploring plant-soil-atmosphere interactions.

Conclusions
Modeling with Century revealed several insights about soil C and nutrient dynamics in this tropical rainforest. Balancing the C, N, and P budgets required additional nutrient inputs, suggesting that processes not currently contained in the model, such as erosion, sequential reduction and oxidation, and changes in soil pH play important roles in soil C and nutrient dynamics in tropical soils with variable-charge clays. Although this forest receives 400 cm MAP, such that one would expect cloudy conditions that could reduce ET, modeled PET was relatively high, and was best explained by the fact that a large proportion of rain falls at night. In our experimental modeling with Century, we found that a biome-based model could reflect underlying plant biology well enough to simulate the effects on ecosystems of changes in plant traits related to biogeochemistry, and thereby provide insights into mechanisms that link plant-soil interactions with ecosystem functioning. Some tree species were modeled more realistically than others that had unusual traits, indicating that novel plant traits and soil processes not usually included in models can frustrate model predictions. These simulations also indicated that relatively subtle changes in a tropical climate, such as increasing night-time temperature, as well as land-use change that involves shifts in species composition, can have feedbacks to the C cycle at the regional level, and potentially the Earth system.

Supplementary Materials:
The following are available online at www.mdpi.com/xxx/s1, Figure S1: Comparison of simulated and measured effects of five tropical tree species on ecosystem processes and properties; Table S1: Site parameters for calibrating Century for the experimental field site at La Selva Biological Station in Costa Rica; Table S2: Century model input parameters regarding tree species and their effects on mature rainforest at La Selva Biological Station, Costa Rica (LSM); Table S3: Estimate of net primary productivity (NPP) in La Selva oldgrowth forest; Table S4: Input file parameters for land-use change prior to and including establishment of plantations at La Selva; Table S5: Monthly air temperature parameters (°C) for modeling experiment regarding increase in night-time temperature (Tmin) using Century, calibrated for the experimental field site at La Selva Biological Station in Costa Rica; Table S6: Monthly precipitation (cm) for modeling experiment regarding climate-change scenario of a drier dry season using Century, calibrated for the experimental field site at La Selva Biological Station in Costa Rica; Table S7: Net primary productivity (NPP) in experimental plantations at La Selva; Table S8: Century input parameter adjustments for individual species in experiment at La Selva biological Station.
Author Contributions: Conceptualization, methodology, software, formal analysis, validation, and funding acquisition: A.E.R. and W.J.P. Data curation, visualization, writing-original draft preparation, project administration: A.E.R. Writing-review and editing: W.J.R. All authors have read and agreed to the published version of the manuscript.