Assessing Impacts of Plant Stoichiometric Traits on Terrestrial Ecosystem Carbon Accumulation Using the E3SM Land Model

Carbon (C) enters into the terrestrial ecosystems via photosynthesis and cycles through the system together with other essential nutrients (i.e., nitrogen [N] and phosphorus [P]). Such a strong coupling of C, N, and P leads to the theoretical prediction that limited nutrient availability will limit photosynthesis rate, plant growth, and future terrestrial C dynamics. However, the lack of reliable information about plant tissue stoichiometric constraints remains a challenge for quantifying nutrient limitations on projected global C cycling. In this study, we harmonized observed plant tissue C:N:P stoichiometry from more than 6,000 plant species with the commonly used plant functional type framework in global land models. Using observed C:N:P stoichiometry and the flexibility of these ratios as emergent plant traits, we show that observationally constrained fixed plant stoichiometry does not improve model estimates of present‐day C dynamics compared with unconstrained stoichiometry. However, adopting stoichiometric flexibility significantly improves model predictions of C fluxes and stocks. The 21st century simulations with RCP8.5 CO2 concentrations show that stoichiometric flexibility, rather than baseline stoichiometric ratios, is the dominant controller of plant productivity and ecosystem C accumulation in modeled responses to CO2 fertilization. The enhanced nutrient limitations and plant P use efficiency mainly explain this result. This study is consistent with the previous consensus that nutrient availability will limit xfuture land carbon sequestration but challenges the idea that imbalances between C and nutrient supplies and fixed stoichiometry limit future land C sinks. We show here that it is necessary to represent nutrient stoichiometric flexibility in models to accurately project future terrestrial ecosystem carbon sequestration.


Introduction
As an important functional trait, plant tissue-level stoichiometric ratios define relative abundances of carbon (C) and other necessary chemical elements (e.g., N and P) in different plant tissues (Watanabe et al., 2007), such as leaves, fine roots, sapwood (live wood), and heartwood (dead wood). These elemental ratios exert strong controls on how energy and material flow through plants and thus play an important role in plant survivorship, growth, reproduction, and functioning (Elser et al., 2010). According to the stoichiometry homeostasis hypothesis (Sterner & Elser, 2002), plants strive to maintain critical tissue stoichiometric ratios for growth and function, even though external element supplies may dramatically change across space and time. To maintain this homeostasis, plant C assimilation can be reduced, for example, when soil nutrient supply is reduced (Agren & Weih, 2012;Harpole et al., 2011). In this case, the "immediate" reduction of plant biomass production could be due to two reasons: (1) direct functional control of nutrients on biochemical photosynthesis reactions and/or (2) biomass construction limitations. The former case is supported by the unique role of N in maintaining important proteins such as the Rubisco enzyme that drives photosynthesis machinery (Kattge et al., 2009;Niklas et al., 2005;Reich et al., 1995). Phosphorus is another essential element, particularly in P-rich ATP and rRNA, that bind CO 2 molecules to RuBP (ribulose bisphosphate), drive cell growth, and facilitate plant metabolism (Hidaka & Kitayama, 2013;Kiirats et al., 2009;Mate et al., 1993;Reich et al., 2009). In addition, plant tissue construction requires both nonstructural C and nutrients following fundamental stoichiometry rules (Elser et al., 2010;Kerkhoff et al., 2006;Sterner & Elser, 2002). Therefore, the reduction of C productivity could occur when nutrient uptake from soil does not keep pace with C uptake from photosynthesis (Hungate et al., 2003).
Consistent with this theoretical understanding of plant C to nutrient relationships under natural conditions, perturbation experiments (i.e., nutrient fertilization) confirm the hypothesis that, in nutrient-limited ecosystems, C productivity is reduced but can be enhanced by supplying additional N (Ares & Fownes, 2001;Foster & Morrison, 2002;Gundersen, 1998), P (Campo & Vázquez-Yanes, 2004;Herbert & Fownes, 1995;Vitousek & Farrington, 1997), or both (Davidson et al., 2004;Sarmiento et al., 2006;Tanner et al., 1990). Moreover, ecosystems could experience different nutrient limitation conditions, due to distinct histories of long-term soil and ecosystem development. For example, high-latitude ecosystems (i.e., Arctic tundra and boreal forest) tend to be more N limited, because of energetic constraints on reactive N supply through N 2 fixation (Vitousek & Field, 1999) and soil organic N mineralization (Nadelhoffer et al., 1991). However, lowland tropical forests are expected to be more P limited due to slow P supply from weathering of low P concentration parent materials and due to long-term depletion via leaching and mineral occlusion (Vitousek et al., 2010). Therefore, the observed plant C:N ratios increase (Martin et al., 2014;Vitousek et al., 1988) and N:P ratios decrease (McGroddy et al., 2004;Reich & Oleksyn, 2004) from low latitudes (warm) to high-latitude (cold) ecosystems.
The atmospheric CO 2 fertilization effect on ecosystem carbon sequestration is commonly represented as the carbon-concentration feedback (β) (Friedlingstein et al., 2006), where β denotes carbon storage sensitivity to atmospheric CO 2 concentration. Combining the experimental and theoretical understanding of CO 2 fertilization and nutrient limitation effects on photosynthesis and plant growth gives rise to a prediction that future carbon-concentration feedbacks will be significantly dampened (less positive β) at global scales due to limited nutrient availability. Previous studies that investigated global nutrient limitations on future terrestrial C sequestration used C-to-nutrient elemental balance approaches applied to C-only simulations of the global carbon budget (Hungate et al., 2003;Wieder, Cleveland, Smith, & Todd-Brown, 2015) or nutrient-enabled earth system land model simulations (Y-P Wang & Houlton, 2009;Zaehle et al., 2015;Zhang et al., 2014). Both approaches demonstrated that nutrient availability would strongly attenuate the CO 2 fertilization effect on terrestrial ecosystem carbon accumulation (Thomas et al., 2015;Zaehle et al., 2014).
A general consensus of widespread N and P limitation on global land C sequestration has been reached from theory, observation, and modeling. However, large uncertainties in quantifying how much land C accumulation would be inhibited due to nutrient limitation arise because (1) external nutrient supplies may change (Wieder, Cleveland, Lawrence, et al., 2015), (2) plant-soil C nutrient interactions and competition are dynamically evolving (De Graaff et al., 2006;Riley et al., 2018;Zhu et al., 2019), and (3) plant stoichiometry  and the flexibility of stoichiometry (Medlyn et al., 2015) when plants utilize the acquired nutrients to construct new biomass differs among PFTs and changes in response to environmental conditions. This study focuses on plant stoichiometry and flexibility that (1) constrain biome level C:N:P stoichiometry using multiple plant trait observational databases and (2) quantify the effect of C:N:P stoichiometry and flexibility on future land C sequestration due to CO 2 fertilization using an earth system land model. We hypothesize that (1) 21st century C accumulation will be reduced by constraints from plant stoichiometric traits and (2) baseline C:N:P ratios and C:N:P flexibility play distinct roles in controlling future C accumulation under CO 2 enrichment through biomass construction and ecosystem nutrient use efficiency (defined as plant net primary productivity divide by nutrient uptake).

Stoichiometry Traits
We synthesized observed plant C:N:P stoichiometry for plant tissues that are often represented in prevailing earth system land models (i.e., leaves, fine roots, live coarse roots, dead coarse roots, live stems, and dead stems). Leaf C:N:P stoichiometry data primarily came from the TRY plant trait database (Kattge et al., 2011), from which we used 40,374 records covering 6,438 species. We either directly used the reported C:N, C:P, and N:P ratios or calculated them when both C and N (or C and P) were available. Fine-root C:N:P stoichiometry data were from the Fine-Root Ecology Database (Iversen et al., 2017), TRY, and additional syntheses from published literatures (supporting information, Data S1), which provided 1,126 observations. Compared with widely observed leaf and fine-root C:N:P stoichiometry, woody tissues, including coarse root and stem wood, were less common in these databases. We, therefore, synthesized additional wood C:N:P ratios from other literature (Data S1) and combined those data with TRY. In total, we acquired 338 observations of wood stoichiometry. Classification of species-level observations into biome level categories was based on site information in the Köppen-Geiger climate code (Peel et al., 2007) (tropical, temperate, boreal, and tundra), leaf type (broadleaf and needleleaf), life span (deciduous and evergreen), photosynthesis pathway (C3 and C4), and plant type (grass, shrub, and forest) ( Figure 1).
In this study, we focused on two particular traits of plant stoichiometry: the C:N and N:P ratios emergent at biome level and overall C:N and N:P flexibility through time and space for all biomes. Temporal changes mainly reflect genotype control on plant C:N:P stoichiometry for different biomes, while spatial heterogeneity reflects species-level variation, soil nutrient availability, and climate conditions (Elser et al., 2011;Markert, 1989).

Earth System Land Model
To assess the impacts of plant stoichiometry on the global C cycle, we used the Earth, Energy, Exascale System Model (E3SM) developed by the U.S. Department of Energy (Golaz et al., 2018). The land model of E3SM used in this study is ELMv1-ECA (Riley et al., 2018;Zhu et al., 2019). ELMv1-ECA considers C, N, and P elemental cycles in plants and soil, including soil biogeochemical reactions, plant allocation and physiology, competition between various consumers, and abiotic processes (e.g., leaching). Primary C processes represented in ELMv1-ECA are as follows: (1) C enters into the terrestrial ecosystem through photosynthesis via gross primary productivity (GPP), which is constrained by leaf-level nutrient content (Kattge et al., 2009;Reich et al., 2009;Walker et al., 2014); (2) a fraction of the GPP is respired by plant maintenance and growth respiration depending on the N concentration of plant tissues and growth costs; (3) plant litter flux transfers C from living biomass to soil through leaf and root phenological cycles, background turnover, tree mortality, and disturbance (i.e., fire); (4) wood debris fragmentation and litter decomposition result in litter C accumulation in soil and soil organic matter formation; and (5) soil C is returned to the atmosphere through heterotrophic (microbial decomposition of litter and soil organic carbon). Represented nutrient cycles closely coupled with C include N mineralized by microbial activity in the soil and N supplied from symbiotic N 2 fixation facilitated by nitrogenase enzyme activity (Cleveland et al., 1999). Because of the large energy cost of nitrogenase synthesis, plant N 2 fixation is inhibited when the fine-root system is able to take up nitrogen more efficiently than N 2 fixation (quantified with marginal N gain per carbon investment on N 2 fixation versus fine-root growth) (Y-P Wang et al., 2010). ELM first calculates the potential cost of carbon for N 2 fixation based on soil temperature (Houlton et al., 2008) and compares it with the carbon cost of root nitrogen uptake (Rastetter et al., 2001). The plant either grows new fine roots to take up nitrogen or directly fixes N 2 , according to the carbon costs. This approach is different from the FUN module implemented in CLM5 (Brzostek et al., 2014;Lawrence et al., 2019), which explicitly calculates the carbon cost of N 2 fixation and subtracts the cost from net primary productivity. Reactive P is generated from parent material weathering, plant phosphatase activity, and P deposition. The partitioning of available N among plants, microbial immobilizers, nitrifiers, and denitrifiers is resolved using the equilibrium chemistry approximation (ECA) competition theory (Huang et al., 2018;Medvigy et al., 2019;B Wang & Allison, 2019;Zhu, Iversen, Riley, Slette, & Vander Stel, 2016;Zhu & Riley, 2015;Zhu et al., 2017;Zhu, Riley, Tang, & Koven, 2016). Similarly, we use ECA to resolve P competition and partitioning among plants, microbial immobilizers, and mineral surfaces. The ELM used in this study (ELMv1-ECA) differs from another ELM version (ELMv1-CTC) in three major aspects: ELMv1-ECA uses flexible C:N:P stoichiometry, dynamic C:N:P resource allocation (Friedlingstein et al., 1999), and ECA nutrient competition. In contrast, ELMv1-CTC assumes fixed C:N:P stoichiometry, constant C:N:P allocation, and relative demand-based nutrient competition.

Modeling Protocol and Experiments
We conducted ELMv1-ECA simulations (1.9°by 2.5°resolution) with a 400-year accelerated decomposition spinup, in which soil organic C turnover rates are accelerated to more rapidly achieve equilibrium . After accelerated decomposition spinup, another 400-year regular spinup was conducted. Both spinup simulations used a constant atmospheric CO 2 concentration (285 ppm) and repeated 1901-1920 Global Soil Wetness Project Phase 3 (GSWP3) 6-hourly climate forcing including temperature, precipitation, relative humidity, wind speed, and surface pressure (Danger et al., 2008). Transient simulations were performed from 1901 to 2005 with transient atmospheric CO 2 concentrations, GSWP reanalysis climate forcing (Danger et al., 2008), transitional N deposition (Lamarque et al., 2005), and phosphorus deposition (Mahowald et al., 2008). To isolate the role of enhanced atmosphere CO 2 on ecosystem C dynamics, simulations between 2006 and 2100 were driven with a 20-year (1986 to 2005) repeated GSWP3 climate forcing and future atmospheric CO 2 concentration from the Representative Concentration Pathway (RCP) 8.5 scenario (Friedlingstein et al., 2006). We selected the worst-case scenario (RCP8.5) increase in atmospheric CO 2 as a test case for the purpose of investigating strong nutrient limitation effects on the C cycle.
Three core simulations were conducted, all of which follow the abovementioned simulation protocol but which differed in stoichiometric traits. The baseline simulation used the default unconstrained plant stoichiometry and did not allow flexibility of plant stoichiometry (BASE). The second simulation used the plant stoichiometry from aboveground and belowground observations and did not allow flexibility of plant stoichiometry (FIXED). The third simulation used both observed plant stoichiometry and observed stoichiometric flexibility (FLEX). To focus on land C sequestration, we used four observationally constrained benchmarks to evaluate model performance at global scale: (1) FLUXNET-MTE (multi-tree ensemble) global gross primary productivity upscaled from FLUXNET observations (Beer et al., 2010), (2) global mean carbon use efficiency (net primary productivity divide by gross primary productivity) estimated from MODIS (Moderate Resolution Imaging Spectroradiometer) products (He et al., 2018), (3) living vegetation biomass including aboveground and belowground components (Saatchi et al., 2011), and (4) Soil Harmonized Database global top 1-m soil C stocks (Hiederer & Köchy, 2011).

Plant C:N:P Stoichiometry
Plant stoichiometry is a complex emergent property that varies across space, time, and species (Reich & Oleksyn, 2004). The variation of plant stoichiometry is controlled by multiple factors including genotype (Elser et al., 2011;Markert, 1989), climate conditions, soil biogeochemistry, substrate nutrient availability, and plant physiology (Agren & Weih, 2012). To use the highly variable plant stoichiometry data collected from observations to inform ELMv1-ECA parameterization, we synthesized the observations into 14 natural plant functional types (PFTs) according to structural, phenological, physiological, and climatic features (Table S1). We did not aim to resolve species-level stoichiometric trait differences which were subject to small-scale variability, for example, in soil nutrient availability. Rather, we grouped different species into 10.1029/2019MS001841

Journal of Advances in Modeling Earth Systems
one category if they belong to the same PFT. By doing this, we acknowledge that the simulated "plant" represented in the global model might differ from locally observed species. However, the PFT framework reduces global model computational complexity, benefits model interpretability (Poulter et al., 2015), and generates a tractable testbed for how stoichiometric traits affect the global carbon C. Therefore, we focus on two statistical properties of the observed stoichiometry data: median (50th percentile) and variability (defined as the range between the 25th and 75th percentiles). Here, we refer to the median and variability as baseline and flexibility of stoichiometry, respectively. We note that flexibility includes observed stoichiometry fluctuations over time, intraspecies and interspecies differences, spatial differences, and potential differences in the age, size, or functional classes of the sampled tissue (Iversen et al., 2017). We use this flexibility range to characterize implications of changing PFT stoichiometry ratios on present-day and future carbon dynamics.
At the PFT level, observationally constrained plant baseline leaf C:N ratio is dramatically different from the ELMv1-ECA default leaf C:N ratio (Figure 2), which was inherited from the Community Land Model (CLM4.5) (Oleson et al., 2013). Unlike the default leaf C:N ratios in the model, the observations show significant differences among grasses, forests, and shrubs. C 3 Arctic grass, C 3 non-Arctic grass, and C 4 grass have the lowest observed baseline leaf C:N ratios, and forest and shrub leaf C:N ratios exhibit large differences. Within the non-grass PFTs, tropical broadleaf evergreen forest has the lowest baseline leaf C: N ratio, indicating a relative abundant N supply from soil. This pattern is consistent with the "temperature-biogeochemistry" hypothesis that high temperatures enhance soil N mineralization and thereby N uptake by roots, and lower leaf C:N ratios are supported by faster cycling of litter with low C:N ratios (Reich & Oleksyn, 2004). As a result, leaf C:N ratios increase as mean annual temperature declines from~30°C in, for example, the tropical forest biome, to~5°C in, for example, the boreal forest biome ( Figure S1a). However, leaf C:N ratios start to decline again when mean annual temperature is below 5°C ( Figure S1a), mainly because high-latitude ecosystems are dominated by low C:N ratio C 3 Arctic grass in the model (Figure 2) rather than a weak N limitation over the Arctic system ( Figure S1b).
Compared with leaf C:N ratios, leaf N:P ratios were much more consistent between default and observationally constrained baseline values (Figure 3). In both scenarios, the tropical broadleaf evergreen and tropical broadleaf deciduous forests have the highest leaf N:P ratios. This pattern implies that, given the same nitrogen supply, tropical forests have the lowest phosphorus supply compared with other ecosystems. This prominent feature is explained by the "soil substrate age hypothesis." Tropical soils are older, infertile, and phosphorus-depleted compared with temperate and Arctic soils (Reich & Oleksyn, 2004). As a result, across the temperature gradient from tropical to Arctic ecosystems, observationally constrained leaf N:P ratios monotonically decrease ( Figure S1b).
Fine-root and woody C:N:P stoichiometry is less-often observed than leaf stoichiometry and therefore largely unconstrained in the default model, which assumes constant stoichiometry ratios for (1) fine roots (C:N = 42 and N:P = 24), (2) livewood (sapwood) (C:N = 50 and N:P = 60), and (3) deadwood (heartwood) (C:N = 500 and N:P = 6) (Oleson et al., 2013;Yang et al., 2014). The observationally constrained fine-root C:N baseline ratios generally agree with the default fine-root C:N ratio (Iversen et al., 2017). However, the observationally constrained fine-root N:P baseline ratios are very different and have larger inter-PFT variation than the model's default values. Furthermore, wood stoichiometry differences between observations and the default model are very large. First, stoichiometry data indicate that forest PFTs have relatively higher wood C:N and N:P ratios than shrub PFTs (Figures 2 and 3). Second, default model values consistently Only median values are showed here as used in the model analyses. ArcC3Gr: Arctic C3 grass; BorBDS: boreal broadleaf deciduous shrub; BorBDT: boreal broadleaf deciduous tree; BorNDT: boreal needleleaf deciduous tree; BorNET: boreal needleleaf evergreen tree; C4Gr: C4 grass; non-ArcC3Gr: non-Arctic C3 grass; TempBDS: temperate broadleaf deciduous shrub; TempBDT: temperate broadleaf deciduous shrub; TempBES: temperate broadleaf evergreen shrub; TempBET: temperate broadleaf evergreen tree; TempNET: temperate needleleaf evergreen tree; TroBDT: temperate broadleaf deciduous tree; TroBET: tropical broadleaf evergreen tree. underestimate sapwood C:N and overestimate sapwood N:P ratios, even though the default C:P ratios can be close to observations. Also, the default model underestimates by more than an order of magnitude the heartwood N:P ratios for forest PFTs. C:N:P stoichiometric flexibility is calculated from the 25th and 75th percentile of the data for each PFT and plant tissues. Due to insufficient data to infer stoichiometry distributions for woody tissues and for some undersampled PFTs (see discussion in section 4), we harmonized the stoichiometric flexibility for all PFTs and modeled inter-PFT and intertissue differences and their impacts on C cycle with ensemble simulations. The percentage flexibility of C:N stoichiometry follows a Gaussian-like distribution with a range of 5% to 45%, while percentage flexibility of N:P stoichiometry is more evenly distributed (Figure 4). Although the probability density distributions are different, the median values of percentage flexibility in C:N and N:P ratios are the same (~25% flexibility). According to the 25th, 50th, and 75th percentiles of the distributions, we conduct nine ensemble simulations with C:N flexibility of 17%, 25%, and 28% and N:P flexibility of 15%, 25%, and 43% (see discussion in section 3.4).

Impacts of Stoichiometry Traits on Present-Day Carbon Cycle
Most Earth system land models participating in the Coupled Model Intercomparison Project Phase 6 (CMIP6) consider the N cycle (e.g., CLM5 in CESM [Lawrence et al., 2019], LM3 in GFDL [Gerber et al., 2010]) and some include P dynamics (e.g., ELMv1-ECA in E3SM , JSBACH-CNP in MPI-ESM [Goll et al., 2012]). Although CMIP protocol standardizes external forcings, including N and P deposition rates (Jones et al., 2016), there are several factors relevant to nutrient cycles that distinguish Earth system model (ESM) simulated nutrient constraints on the C cycle. Here, we focus on a dominant factor: how plant C:N:P stoichiometry coupling affects whole ecosystem C dynamics. Existing land models use PFT-based C:N:P stoichiometry to drive the plant C and nutrient coupling. However, the PFT-level C:N:P ratios are commonly derived from individual studies, empirical knowledge, or small datasets (Goll et al., 2012;Sitch et al., 2003). To inform C cycle uncertainty stemming from assumed C:N:P ratios, we drive ELMv1-ECA with two sets of plant C:N:P stoichiometry. The BASE simulation uses default C:N:P stoichiometry based on empirical knowledge and a small dataset (White et al., 2000;Yang et al., 2014). The FIXED simulation uses fixed C:N:P stoichiometry derived from the TRY plant trait database (Kattge et al., 2011), the Fine-Root Ecology Database (Iversen et al., 2017), and our new synthesis (Data S1).
The role of flexible plant stoichiometry has been shown to be important in understanding and modeling ecosystem properties and functions, for example, nutrient fertilization responses (Meyerholt & Zaehle, 2015).
Here we further demonstrate that considering plant stoichiometric flexibility leads to much better model performance in term of major global C fluxes and pools (Zhu & Zhuang, 2015). The FLEX simulated GPP

Journal of Advances in Modeling Earth Systems
(119 PgC year −1 ) is much larger than the other two scenarios and is close to the GPP benchmark (118 Pg C year −1 ) upscaled from FLUXNET in situ observations. CUE in the FLEX and BASE simulations are comparable (~42%), although it is still lower than CUE derived from MODIS products (~48%). Taking GPP and CUE together, allowing flexibility in plant stoichiometry (FLEX) and using observationally constrained C:N:P baseline values leads to more productive (higher GPP) and more efficient (higher CUE) land ecosystems. Flexibility in plant stoichiometry also benefits C accumulation in living biomass and soil due to faster biomass growth and higher litter and coarse woody debris inputs. Therefore, the top 1-m SOC stock in FLEX (970 Pg C) is higher than the other two scenarios and comparable to the 1,050 Pg C benchmark (Nachtergaele et al., 2010).

Impacts of Stoichiometry Traits on Future Carbon Cycle
The observationally constrained plant stoichiometry traits improved model performance in representing critical C processes at present day and could potentially provide more reliable implications under future climate change scenarios. Here, we focus on CO 2 fertilization effects on ecosystem C dynamics and their dependency on plant stoichiometry traits. Therefore, from 2006 to 2100, ELMv1-ECA is driven by repeated historical climate forcing and RCP8.5 CO 2 concentrations. We hypothesize that C:N:P baseline ratio versus flexibility traits play distinct roles in controlling future carbon accumulation under RCP8.5 CO 2 enrichment conditions through plant nutrient use efficiency (defined as plant NPP divide by nutrient uptake). Stoichiometric flexibility in terrestrial vegetation may be able to partly overcome nutrient deficits and maintain larger land C sinks in the CO 2 -enriched future.
Of the three scenarios, FLEX had the largest net ecosystem C gain (331 Pg C) from 2006 to 2100, followed by FIXED (256 Pg C), which maintained constant plant stoichiometry. For these scenarios, the gained C is mostly in SOC (34-37%), then in biomass and coarse woody debris (27-31%), and then in litter (9%) (Figure 6). The net C gain mainly resulted from the CO 2 fertilization effect on plant photosynthesis. Although CUE declines, net primary productivity is still enhanced by higher CO 2 concentrations ( Figure S2).
Further, C accumulation in living biomass mainly occurred over tropical land surfaces for all three scenarios (Figure 7a). In turn, the C accumulation in other tissues (litter, coarse wood debris, and soil) was widely distributed across different ecosystems (Figures 7b-7d). Due to high temperatures and abundant rainfall, tropical soil biogeochemistry cycles nutrients faster than cold or dry regions, thus tropical ecosystems are most efficient in terms of recycling nutrients to accumulate C in plants and soils. Although tropical soil P supply is limited, the flexibility of plant N: P ratio partly offsets P limitation under elevated CO 2 and results in larger vegetation growth (Figure 7) (Fleischer et al., 2019). Consistently, tropical ecosystem P uptake enhancement in FLEX is small (<10%), while the plant phosphorus use efficiency enhancement is larger (35%) (Figure 8d). Unlike P, N uptake over the tropics is enhanced by 30%, indicating imbalanced N and P supplies in tropical ecosystems. Stoichiometric flexibility in this case served as an important strategy for plants to adapt to these source elemental imbalances (Sistla & Schimel, 2012).
By comparing the FIXED and FLEX stoichiometry cases, we find that N uptake, nitrogen use efficiency, and P uptake positively respond to CO 2 concentration change (Figure 8), implying a stronger modeled C-N and C-P coupling. Furthermore, the strength of seasonal plant C-N coupling does not significantly change over the 21st century (Figure 9a), while the C-P coupling strengthens slightly over temperate ecosystems and dramatically over tropical ecosystems (Figure 9b). Due to strong seasonal cycles in temperate ecosystems, C-N and C-P coupling also significantly strengthened during summer when plant growth and nutrient demands are high (Figure 9).
Stoichiometric flexibility is an important mechanism for plants to maintain productivity and function when external elemental supplies (e.g., C, N, P) change rapidly and in an imbalanced way (Sistla & Schimel, 2012). Here, we demonstrated large (29% more land C accumulation; Figure 9) impacts of plant C:N:P flexibility on 21st century terrestrial ecosystem C accumulation. However, we acknowledge that the degree of such stoichiometric flexibility is uncertain and poorly constrained by observations. Therefore, to estimate uncertainty  in 21st century C accumulation associated with uncertain stoichiometric flexibility, we conducted sensitivity analyses using 25th, 50th, and 75th percentiles of the probability distributions for C:N and N:P percentage flexibility (Figure 4). Nine ensemble simulations with combinations of C:N flexibility (17%, 25%, 28%) and N:P flexibility (15%, 25%, 43%) show that the sensitivity of simulated future GPP and NPP to stoichiometric flexibility ( Figure S2, gray lines) is 85% less than the sensitivity from whether or not stoichiometric flexibility is allowed ( Figure S2, red versus magenta lines). Consistently, simulated future GPP continues to increase (39%) due to CO 2 fertilization effects, while the CUE declines (7%) throughout the 21st century. Furthermore, the decline of CUE is dependent on the flexibility range of plant stoichiometry, particularly during the second half of the 21st century.

Reduced Land Model Uncertainty
Most ESM land models that participates in the Coupled Model Intercomparison Project Phase 6 (CMIP6) project consider nutrient cycles (Jones et al., 2016). Predicted terrestrial C budgets and climate feedbacks will

Journal of Advances in Modeling Earth Systems
be dramatically different from C-only versions because plant C production and allocation are strongly limited by N or P (or both) (Wieder, Cleveland, Smith, et al., 2015;Zaehle et al., 2015). However, unlike phytoplankton stoichiometry that is tightly constrained by the "Redfield" C:N:P ratio (106:16:1) (Takahashi et al., 1985), the stoichiometry of terrestrial plants varies across biomes and has large variation across time and space (Meyerholt & Zaehle, 2015;Norby & Iversen, 2006). Compared with observationally constrained fixed stoichiometry, we show that using unconstrained prior stoichiometry could significantly bias terrestrial C accumulation rates at present day and in the future (Figures 5 and S2). To avoid likely large biases in climate predictions induced by model-specific stoichiometry, we suggest using a harmonized and data-constrained PFT-specific C:N:P stoichiometry, which we have provided here (Tables S1 and S2).
Reducing ESM land model uncertainty with observationally constrained stoichiometry traits also requires continuous effort to collect and process plant tissue samples, including standardizing trait definitions and measurement protocols. For example, observed fine-root C:N ratios may reflect "absorptive" parts of the fine-root system in one study (i.e., root order 1, finest root) versus "transport" parts of the fine-root system in another study (i.e., root orders 4 and higher) (Iversen et al., 2017), thus generating ambiguity and uncertainty. Also, available stoichiometry data used in this study are mostly for leaves, and much less information is available for fine-root and woody tissues. Future field campaigns focused on fine-root or wood stoichiometry can improve confidence when applying observationally constrained stoichiometry in ESM land models. In addition, future stoichiometry measurements should prioritize undersampled regions and PFT groups, for example, tropical rain forests, Eurasian boreal forests, and Arctic tundra ecosystems, and less-measured tissue P concentrations that also require more plant tissue for analysis ( Figure 1). Furthermore, in order to effectively make use of those observational stoichiometry data, effective scaling the data from local to larger scales is also important. In this study we scaled the site level data to PFT level, which is compatible with ESM PFT-based plant physiology parameterizations (Poulter et al., 2015). Another possible scaling could be from site level to grid cell level, so that the land model could resolve spatial variation of vegetation dynamics within the same PFT group but across different soil environments (e.g., fertility) (Zhu & Zhuang, 2013). However, the latter approach requires more available site level data and a higher spatial coverage. Further investigating relationships among aboveground and belowground chemistry may allow the prediction of root chemistry from more easily measured leaf chemistry (Liu et al., 2010).

Limitation and Future Work
Although plant C:N:P stoichiometric traits, including baseline ratios and flexibility, are both demonstrated to affect land C dynamics, these conclusions have been made with the assumptions that (1) emergent stoichiometric traits at the PFT level well represent the diverse plant species within that group across space and time and (2) stoichiometry flexibility is the same in different plant tissues. The first assumption might be valid under stable environmental conditions. However, in future simulations when environmental drivers change, plant species composition might significantly shift and the PFT level C:N:P stoichiometric traits could be inconsistent with those observed at present day (van Bodegom et al., 2014). Future work should supplement the stoichiometry traits dataset with perturbation experiment data (i.e., CO 2 enrichment, warming, and disturbance experiments) that measure the continuous changes of plant species composition, stoichiometry traits, and C dynamics. The second assumption is partly due to data unavailability (particularly for woody components) and partly due to the fact that ELM has no preferential allocation (resources are proportionally allocated to different plant tissues based on allocation fraction). Although our sensitivity analysis ( Figure S2) shows the limited impacts of this assumption, future work is needed to better understand how leaf, fine roots, and stem stoichiometry traits are distinct from each other (Medlyn et al., 2015).

Conclusions
Previous efforts demonstrated the importance of nutrient constraints on future carbon accumulation, from a C, N, and P mass balance point of view with prescribed and fixed plant C:N:P stoichiometry (Cleveland et al., 2013;Hungate et al., 2003;Wieder, Cleveland, Smith, et al., 2015). In this study, we argue that ecosystem C accumulation constraints due to nutrient supply will be partly alleviated, because N and P biogeochemical cycles are accelerated and ecosystem nutrient use efficiency is enhanced when flexible stoichiometry is considered. Particularly over tropical ecosystems, elevated P use efficiency led to large C sinks in the future under CO 2 enrichment scenarios. We also show that ELMv1-ECA prescribed with 10.1029/2019MS001841

Journal of Advances in Modeling Earth Systems
different fixed plant C:N:P stoichiometry (default versus observationally constrained) led to large differences in simulated C dynamics. This result suggests that uncertainty stemming from plant stoichiometry ratios could result in large uncertainty in CMIP6 simulations. We therefore encourage focused, interdisciplinary efforts to develop detailed CNP earth system land models, measure undersampled plant species and tissues, and perform syntheses to assemble newly available and unpublished data to improve the quantification of plant stoichiometry traits in earth system land models.