Some ESMs or their land components, the land surface models (LSMs), now include some representation of the N and P cycles, reflecting the relevance of nutrients as an important limiting factor across much of the Earth's land C assimilation (Fleischer et al., 2019;LeBauer & Treseder, 2008;Terrer et al., 2019). In the latest coupled model intercomparison project phase 6 (CMIP6; Eyring et al. (2016)), at least 10 ESMs considered N cycling (Arora et al., 2020;Jones et al., 2016), although only one-ACCESS/CABLE (Kowalczyk et al., 2013;Wang et al., 2011)-included P dynamics (Goll et al., 2012;Thum et al., 2019;Wang et al., 2007;Yang et al., 2019;Zhu et al., 2019). Arora et al. (2020) suggested that model uncertainty is reduced in the CMIP6 simulations when N cycling is included. However, most of these models assume that N and P are unrealistically acquired at no C cost to plants ( Davies-Barnard et al., 2020), or these costs are to an extent implicitly included in autotrophic respiration, but without explicitly representations. This cost is often thought of as C assimilates that would otherwise be allocated as C to plant growth, that is, net primary production (NPP). In fact, omitting the C costs of N and P acquisitions may have significant impacts on the total C uptake estimates, as well as on the size and dynamic of N and P predicted by ESMs (Fisher et al., 2010;Ostle et al., 2009;Shi et al., 2016).

Autotrophic respiration is generally treated in LSMs as maintenance and growth respiration fluxes separately. Maintenance respiration is defined as the C cost to support the metabolic activity of plant living tissues, while growth respiration is defined as the additional C cost for the synthesis of new growth (Oleson et al., 2013). While maintenance respiration is associated with protein and membrane turnover and maintenance of ion concentrations and gradients, growth respiration is associated with growth processes such as synthesis of new structures in growth (Lötscher et al., 2004). Note that neither maintenance, nor growth respiration terms explicitly account for the C costs of nutrient acquisition. The most common pathways in which a plant expends energy or C costs for nutrient acquisition are resorbing (or retranslocating, these terms are used interchangeably) nutrients from senescing leaves, mobilizing and taking up nutrients that are locked in mineral soils or organic matter, maintain symbioses in the rhizosphere (mycorrhizae) through C exudates (Hobbie, 2006;Hobbie & Hobbie, 2008;Högberg & Högberg, 2002;Parniske, 2008), and biological N fixation (Dickinson et al., 2002;Guignard et al., 2017;Rastetter et al., 2001;Shi et al., 2016;Vitousek et al., 2002;Vitousek & Field, 1999;Wang et al., 2007). While it is well-known that plants allocate a part of the available C to acquire N and P from multiple acquisition pathways, only a few studies have attempted to quantify the variable C costs of N and P acquisition at the same time from each uptake pathway separately (Allen et al., 2020;Brzostek et al., 2014;Fisher et al., 2010;Shi et al., 2016). Some evidence suggests that trees can allocate up to 20% of NPP to symbiotic and free-living microbes in the rhizosphere to increase their access to N (Brzostek et al., 2015;Hobbie, 2006;Högberg & Högberg, 2002). By accounting for N and P limitation at the same time, some ESMs project reductions in global NPP by as much as 25% (Allen et al., 2020;Goll et al., 2012;Wang et al., 2010), with some models suggesting that the land surface will become a source of C to the atmosphere by the end of the 21st century (Wieder et al., 2015a), although still highly debatable (Brovkin & Goll, 2015;Sun et al., 2017;Wieder et al., 2015b).

To represent C costs of N acquisition in ESMs, Fisher et al. (2010) introduced the fixation and uptake of nitrogen (FUN1.0) model, which is based on an optimal allocation theory whereby plants optimize the allocation of C used to acquire nutrients from the soil (directly through roots or from mycorrhizal associations), senescing leaves, and symbiotic biological N fixation to maximize growth (Bloom et al., 1985;Rastetter et al., 1997). Brzostek et al. (2014) integrated C-N trade-offs of arbuscular mycorrhizae (AM) and ectomycorrhizae (ECM) fungi into FUN (FUN2.0), which significantly improved dynamic predictions of N retranslocated from leaves and taken up from the soil. Afterward, Shi et al. (2016) coupled FUN2.0 into the Community Land Model version 4.0 (CLM4) (Bonan et al., 2011;Koven et al., 2013;Oleson et al., 2010) and Cai et al. (2016) coupled FUN into Noah-MP (Niu et al., 2011). FUN2.0 was then further developed and carried into CLM5 (Fisher et al., 2019;Lawrence et al., 2019;Wieder et al., 2019).

The incorporation of P dynamics into FUN2.0, mirroring the existing model structure for N, led to a new model version, FUN3.0 (with a name update: FUN, Allen et al., 2020). Here, we couple FUN3.0 into the energy exascale Earth system model (E3SM) land model (ELM; Caldwell et al., 2019;Golaz et al., 2019;Yang et al., 2019;Burrows et al., 2020) to explore the sensitivity of the land surface C uptake to a dynamic C cost of N and P acquisition at the same time at global level. We evaluate variability across plant functional types (PFTs) to understand which ecosystems are the most N and P limited. We investigate how the C costs of N and P acquisition vary spatially, and determine how sensitive the land C sink is to a dynamic prediction of the C costs of N and P acquisition.

Within FUN, plants use C for N or P uptake by preferentially using dynamically varying pathways that are "cheaper" to them. Decreasing nutrient availability in soil and leaves leads to higher C expenditure on nutrient uptake and thus reduced growth (Fisher et al., 2019). In the original FUN1.0 model, plants obtain N through the following alternative and most common pathways: (a) active uptake, (b) retranslocation, and (c) symbiotic biological N fixation, which presents concurrent costs in terms of C, determined by the abundance of N in soils and plant tissues, and the temperature-dependent cost of N fixation (Bloom et al., 1985;Fisher et al., 2019;Jiang et al., 2017;Rastetter et al., 2001;Riahi et al., 2017;Terrer et al., 2018). The full model description and further analysis can be found in previous publications (Allen et al., 2020;Brzostek et al., 2014;Fisher et al., 2010;Shi et al., 2016). The full set of model equations can be found in Text S1 in Supporting Information S1.

Following the same approach as the original FUN1.0 model, FUN2.0 includes differential costs of active N uptake between ECM, AM, and non-mycorrhizal plants building upon the simultaneous N uptake from different pathways such as resistors in parallel using the Ohm's law (Brzostek et al., 2014). Allen et al. (2020) integrated P dynamics into FUN2.0 by incorporating the direct C cost of P uptake on top of the C cost of N uptake into a new model formulation, that is, FUN3.0. FUN3.0 was validated over 45 temperate forest plots in Indiana, USA, and 18 tropical dry forest plots in Guanacaste, Costa Rica. The validation sites varied in P availability and mycorrhizal associated tree types. FUN3.0 provides an important framework for predicting coupled N and P limitation within ESMs to enhance predictions of ecosystem response to global change.

The optimal allocation framework of FUN2.0 has been coupled into CLM4.5 and CLM5 (Brzostek et al., 2014;Fisher et al., 2019;Lawrence et al., 2019;Shi et al., 2016Shi et al., , 2019)). Here, we take a similar approach to integrating the P dynamics of FUN3.0 into the E3SM land model (ELM: Ricciuto et al., 2018;Riley et al., 2018;Holm et al., 2020). ELM-FUN is based on the Converging Trophic Cascade (CTC) biogeochemistry framework, which includes P dynamics (Yang et al., 2016) and plant storage pools as new capabilities.

Root biomass impacts nutrient root uptake as described in Equation S3 in Supporting Information S1 and nutrient allocation distributed among different plant organs follows Metcalfe et al. (2017). For a complete description of the CNP cycles in ELM, including plant allocation and other fluxes please refer to Burrows et al. (2020).

Features in ELM include: (a) prognostic P cycle and C-N-P interactions; (b) mortality rates for tropical forests based on soil physical factors (Galbraith et al., 2013); and (c) C, N, and P non-structural vegetation storage pools. The P pools and fluxes are based on a previous modeling study that implemented the P cycle into CLM4 (Yang, Thornton, et al., 2014). The P cycle is vertically resolved, capturing the dynamics of inorganic P pools, the conversion between organic P and inorganic P in each soil layer, as well as the vertical transport of P pools between layers (Yang et al., 2019). Allocation to structural plant pools is calculated using available C and available non-structural nutrients after the remaining C (assimilation minus maintenance respiration) is allocated to a non-structural C storage pool following Metcalfe et al. (2017). Flexible plant tissue stoichiometry is not considered in the model.

Vegetation is divided into multiple PFTs with independent parameterizations controlling photosynthesis, leaf gas exchange, biomass allocation, and other processes (Sulman et al., 2021). The current PFT configuration of ELM follows CLM5; Lawrence et al. (2019), which are based on broad plant groups relevant for global-scale configurations (Wullschleger et al., 2014).

To generate the trade-offs between AM, ECM, and non-mycorrhizal root uptake, FUN3.0 needs an estimate of the percentage of aboveground biomass that associates with each mycorrhizal type for each model grid cell (Braghiere, Fisher, et al., 2021;Brzostek et al., 2014;Shi et al., 2016). Each PFT was classified based on commonly known associations between the specific plant species and either AM-and ECM-fungi. As a result, some PFTs are largely AM-dominated (e.g., grasslands, crops) and some are ECM-dominated (e.g., boreal forest) (Allen et al., 1995;Phillips et al., 2013;Read, 1991).

These PFT fractions are coarse and do not extensively capture the spatial heterogeneity in mycorrhizal association, particularly in mixed-mycorrhizal PFTs, such as tropical (Waring et al., 2016) and temperate forests (Phillips et al., 2013). However, recent comparisons with a global map of mycorrhizal association (Steidinger et al., 2019) indicate the robustness of this assumption (Braghiere, Fisher, et al., 2021). Regarding N fixation, following the implementation in CLM5 (Fisher et al., 2019) each PFT has an associated maximum fraction of the C acquisition that can be used for fixing N, defined as a constant fraction for natural PFTs (25%). The fractions of the AMand ECM-associated plants in association with ELM PFTs are shown in Table S3 in Supporting Information S1. The spatial patterns in nutrient uptake pathways are driven by a combination of the PFT distribution plus any other non-PFT spatial information that impacts the FUN model. In LSMs, a vegetated unit with multiple PFTs is associated with a single soil column. FUN3.0 was coupled to ELM using thirteen input variables (Table S1 in Supporting Information S1) with most of them being directly calculated by the C-N-P processes of ELM including: available soil N (N soil ) and available soil P (P soil ), fine root biomass (C root ), leaf N before senescence (N leaf ), leaf P before senescence (P leaf ), NPP (C NPP ), soil temperature (T soil ), and evapotranspiration (ET). The derivation and direct usage of variables from ELM to generate the necessary inputs for FUN3.0, as well as the key modifications needed to enable the coupling between FUN3.0 and ELM can be found in Supporting Information S1. A diagram of the coupled model is shown in Figure 1.

Global simulations were performed at hourly 1.9° × 2.5° spatiotemporal resolution. All three model simulations were spun up using ELM for 200 years of accelerated decomposition mode, followed by 600 years of regular spin up (Koven et al., 2013;Thornton & Rosenbloom, 2005), while soil P pools were initialized from observations (Yang et al., 2013) at the beginning of the regular spin up. Some processes governing the global distribution of soil P stocks operate in geologic time scales, which makes it impossible at this time to perform a fully prognostic model spin up from "bare ground" conditions to initialize P pools with slow turnover times. Instead, ELM relies on global maps of soil phosphorus pools synthesized from observations (Yang & Post, 2011;Yang, Post, et al., 2014), which have previously been used successfully in land model simulations using prescribed atmospheric forcings (Burrows et al., 2020;Yang et al., 2016). Spin up was followed by historical transient simulations (1850-2010) with historical atmospheric CO 2 , N deposition (Lamarque et al., 2005), P deposition (Mahowald et al., 2005), land use change (Hurtt, 2018), and using the Global Soil Wetness Project 3 (GSWP3) climate forcing (Dirmeyer et al., 2006). We analyzed a 12 year (1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005) period and investigated the role of N and P in affecting the responses of NPP to three configurations of the model: ELM only (CNP cycles without the FUN model), ELM-FUN2.0 (CNP cycles with the FUN model for N only), and ELM-FUN3.0 (CNP cycles with the FUN model for N and P). These factorial experiments allow us to disentangle the impacts of N and P in isolation and in combination.

We use the International Land Model Benchmarking (ILAMB) v2.6 package for model validation (Collier et al., 2018) focusing on global patterns of the C and water cycle variables, including datasets of: (a) aboveground living biomass based on inventory plots upscaled to the globe using remote sensing imagery (Saatchi et al., 2011); (b) emulated CO 2 concentrations from the NOAA/ESRL flask network (http://www.esrl.noaa.gov/gmd/ccgg/); (c) FLUXCOM GPP and ecosystem respiration (Jung et al., 2019(Jung et al., , 2020)) -Brown et al., 2013); and (viii) evapotranspiration (ET) from the Global Land Evaporation Amsterdam Model (GLEAM) v3.3a (Martens et al., 2017).

The relationships of the evaluated variables with precipitation, temperature, and surface incident shortwave radiation were performed using data products from the global precipitation climatology project (GPCP) Monthly Analysis version 2.3 (Adler et al., 2018), the climatic research unit (CRU) monthly temperature version 4.02 (Harris et al., 2014), and the clouds and the earth's radiant energy system (CERES) surface irradiances edition 4.1 (Kato et al., 2018;Loeb et al., 2018).

We show the results for the ILAMB overall score for the absolute values (S overall ), as well as their relationships with precipitation, incident shortwave radiation, and temperature ( 𝐴𝐴 𝐴𝐴 relationship overall

) defined as:

where S bias is the spatially integrated bias score, S RMSE is the root-mean-squared error (RMSE) score doubly weighted to emphasize its importance, S phase is the phase shift score, S iav is the interannual variability score, and S dist is the spatial distribution score in Equation 1.1. 𝐴𝐴 𝐴𝐴 x RMSE is the spatial RMSE score for each x meteorological variable in Equation 1.2. For the whole set of equations of each term in Equations 1.1 and 1.2, and further details refer to Collier et al. (2018).

NPP results from our simulations are also compared with the NPP from 11 CMIP6 models (the same 11 models evaluated in Arora et al. (2020)), MODIS (Running & Zhao, 2015;Running et al., 2004;Zhao et al., 2005), and the International Satellite Land-Surface Climatology Project, Initiative II (ISLSCP II) collection from the International Geosphere Biosphere Program (IGBP) (Cramer et al., 1999). Projects such as the CMIP6 (Eyring et al., 2016) have allowed for consistent comparison of the response of the C cycle under climate change from existing state-of-the-art ESMs. Using CMIP6 outputs together with satellite observations allow us to estimate the range uncertainty in C cycle variables due to model structural variation. We note that our simulations are done offline with atmospheric forcing data, such as in LUMIP (Lawrence et al., 2016) and TRENDY (Sitch et al., 2015), while CMIP models also simulate atmospheric states in coupled runs through the historical period.

In addition, we compared estimates of biological N fixation with a recent compilation of 12 studies by Davies-Barnard and Friedlingstein (2020), plus estimates by Sullivan et al. (2014) and Peng et al. (2020). We also added an explicit comparison with predicted N and P limitation values presented in Du et al. (2020) and generated an agreement map between both data products.

Coupling both recent versions of FUN into ELM improved the overall normalized standard deviation between model and observation for all the evaluated variables in ILAMB (Figure 2). Moreover, FUN3.0 improves model performance even further when compared to FUN2.0. Although most of the evaluated variables were improved to some degree, the largest improvements were found for biomass (Figure S1a in Supporting Information S1), LAI (Figure S1d in Supporting Information S1), and soil C (Figure S1f in Supporting Information S1).

The relationships of the evaluated variables with precipitation, incident shortwave radiation, and temperature also improved with ELM-FUN3.0 (Figure 2). All the ILAMB plots comparing ELM, ELM-FUN2.0, and ELM-FUN3.0 are available from https://braghiere.github.io/ILAMB_ELM_FUN_reduced_v1/. [1994][1995][1996][1997][1998][1999][2000][2001][2002][2003][2004][2005], which is 13.6 Pg C yr -1 less than that simulated by ELM for the same period. This implies a downregulation of global annual NPP by 21.0%. Zonally, NPP was reduced across all latitudes in ELM-FUN2.0 relative to ELM (Figure 3a). In ELM-FUN3.0, NPP was 32.2 PgCyr -1 on average from 1994 to 2005, which is 32.6 PgCyr -1 less than that simulated by ELM for the same period, a reduction of 50.3% in global annual NPP. Zonally, NPP was reduced across all latitudes in ELM-FUN3.0 relative to ELM (Figure 3). The drops in global NPP are larger than the direct costs of nutrient uptake due to feedbacks on LAI and GPP in the model. Additionally, the global mean C Use Efficiency (CUE; defined as NPP/ GPP) is negatively impacted when considering the C costs of N and P acquisition (Figure S11b in Supporting Information S1), going from ∼30% in ELM to ∼20% in ELM-FUN3.0. Although CUE is a highly uncertain metric, previous studies indicate that natural ecosystems present a CUE of 32%-47% (Campioli et al., 2015;Malhi et al., 2009).

The NPP zonal mean (along the latitudinal component) from all three simulations with ELM are shown in scatter plots compared to the zonal mean of MODIS NPP (Figure 3b), IGBP NPP (Figure 3c), and the zonal mean NPP from the ensemble member of 11 CMIP6 models (Figure 3d). The implementation of FUN into ELM improved the relationships between ELM and all three data products. For MODIS NPP, considering the C costs of N acquisition in ELM increased the r 2 of the linear relationship from 0.32 to 0.84 and reduced the root mean squared error (RMSE) by 52% from 223.2 gCm -2 yr -1 to 107.7 gCm -2 yr -1 , while the slope also changed from 1.26 to 1.03. Adding the C costs of P acquisition on top of the costs of N acquisition increased r 2 to 0.86 and reduced the RMSE by 54%-102.0 gCm -2 yr -1 , while the slope changed to 0.73.

For the IGBP NPP linear relationship, the RMSE decreased 56% from 152.4 gCm -2 yr -1 with ELM to 67.6 gCm -2 yr -1 with ELM-FUN2.0%, and 2% to 148.61 gCm -2 yr -1 with ELM-FUN3.0. For the CMIP6 ensemble member NPP, the RMSE decreased 42% from 142.9 gCm -2 yr -1 with ELM to 82.6 gCm -2 yr -1 with ELM-FUN2.0, but increased 9% (155.9 gCm -2 yr -1 ) with ELM-FUN3.0. The r 2 , RMSE, and slope of all model runs and products are shown in Figure 3. In general, the relationships between modeled zonal NPP with ELM and these three data products seem to be improved when the C costs of nutrient acquisition are considered (Figure 3b), although the addition of both nutrients seems to underestimate NPP, potentially because of implicit compensating mechanisms already present in the ELM NPP without FUN. Although FUN3.0 directly represents the C costs of N + P acquisition, that does not mean the remaining parts of ELM are improved as well. This study does not intend to improve all the processes represented in the model. FUN is dynamically coupled to ELM in the historical runs, that is, 1850-2010. However, during spin up (200 years in accelerated decomposition mode + 600 years in normal mode), the CNP pools are brought to equilibrium without the FUN model, that is, ELM only for all three model configurations as described in Section 2.2. The standard spin-up protocol is used to achieve C, water, and energy equilibrium at the start of the simulation (Lawrence et al., 2019). The specifics of a spin up method often vary according to model type and objectives, but generally it involves either running a model repeatedly using a single period of forcing data, or running a model for many years with historical forcing data and evaluating model outputs after a time window long enough to achievement of steady state of the variable of interest.

The year 1850 equilibrium conditions were calculated by integrating over a repeating 20-year period of an atmospheric reanalysis data set (i.e., years 1901-1920 from the forcing data sets) along with fixed atmospheric CO 2 , N deposition, and P deposition. As with earlier versions of ELM, it is prohibitively expensive to run the full model for the period of time required to achieve a quasi-steady state. Thus, the spin-up procedure involves an "accelerated decomposition" methodology introduced in Koven et al. (2013) and Thornton and Rosenbloom (2005). During the accelerated decomposition phase, the decomposition of the slow C pools (e.g., the long turnover time soil C and coarse woody debris pools) are artificially increased to allow faster convergence on the equilibrium state (Lawrence et al., 2019). Due to the processes in FUN accessing soil P, which operates in geologic time scales, ELM-FUN coupling starts after model spin up. Therefore, it is expected that in the initial years of the transient simulation, simulations with ELM-FUN versions show a pulse in autotrophic respiration during the first years because of a new C source to autotrophic respiration from FUN, but over time respiration declines relative to ELM (Figure S10 in Supporting Information S1). This behavior is also described in Shi et al. (2016) when evaluating the impacts of FUN in CLM. Nevertheless, the evaluations performed in this study focus on the period 1994-2005, or 144 years after model equilibrium was reached in order to minimize the possible impacts of the model initialization from cold start (Meier et al., 2018;Shi et al., 2016).

Among the few global nutrient cycling products available for model benchmarking, Fisher et al. (2012) coupled remotely sensed estimates of NPP and ET to identify areas where light and water are not limiting factors on productivity. Overall, the NPP downregulation fraction simulated by ELM-FUN3.0 presents a similar global spatial distribution to that of the nutrient limitation product from Fisher et al. (2012) (see Figure S2 in Supporting Information S1). Figure 4b shows a scatterplot between the zonal mean of nutrient limitation product from Fisher et al. (2012) and the zonal mean of the normalized C use ratio for N acquisition in FUN2.0, and N + P acquisition in FUN3.0. The C use ratio is defined as the ratio of C spent on N acquisition (FUN2.0) or N + P acquisition (FUN3.0) to the total C available to the plant. The normalized C use ratio is the C use ratio divided by the maximum C use ratio. The agreement between the modeled nutrient limitation and the observation-based estimate from Fisher et al. (2012) increases when N and P are considered together. The r 2 increases from 0.73 to 0.83 and the RMSE decreases from 4.3% to 3.7%.

Globally, the amount of N uptake and the pathway of N acquisition varied as a function of ecosystem productivity, dominant mycorrhizal type, soil nutrient availability, and environmental conditions. ELM-FUN3.0 simulated the mean annual global total N uptake as 841.8 Tg N yr -1 .

Mycorrhizal uptake and direct root uptake follow the spatial distribution of the most productive areas of the globe, but with significantly different magnitudes (Figure 5). While the global total mycorrhizal root uptake is 659.9 Tg N yr -1 , mostly due to AM fungi associations (73.0% of the total mycorrhizal root uptake or 482.1 Tg N yr -1 ), with the remaining 27.0% due to ECM fungi association (177.8 Tg N yr -1 ) (Figure 5); the direct root uptake is 84.3 Tg N yr -1 . These results are due to the differences in the C cost of uptake by mycorrhizal roots compared to non-mycorrhizal roots as defined in Table S2 in Supporting Information S1 and Allen et al. (2020), the spatial distribution of available soil N (Figure S3 in Supporting Information S1), and the spatial distribution of mycorrhizal types. Spatial differences among distinct mycorrhizal types are given by the PFT distribution following Table S3 in Supporting Information S1. The extra-tropics have lower rates of N uptake in comparison to the tropics, with higher values associated with mycorrhizal uptake and retranslocation in boreal forests, and direct root uptake over the tundra ecosystems. Arid and semi-arid areas of the world present the lowest N uptake (Figure 5).

Total biological N fixation comprises 35.3 Tg N yr -1 , nearly an order of magnitude lower than the total of mycorrhizal and non-mycorrhizal root uptake. In the northern hemisphere, temperate forests have higher symbiotic N fixation values than other ecosystems (Figure 5). The highest rates of symbiotic N fixation are in the tropics. This estimate is near the lowest boundary of a global compilation study suggesting that biological N fixation is between 52 and 130 Tg N yr -1 (Davies-Barnard & Friedlingstein, 2020) (Figure S8 in Supporting Information S1). Other studies highlight the need for more field measurements of N fixation (Peng et al., 2020), and while our low N fixation estimate is likely due to the fact that all fixed N goes into the plant directly, rather than going into soils first such as N from other active uptake pathways, a study presents evidence that N fixation may be lower than previously assumed (Sullivan et al., 2014).

The total amount of retranslocated N and P from leaves globally are 97.6 Tg N yr -1 and 7.3 Tg P yr -1 , respectively. Tropical and temperate forests, as well as grasslands and crop areas in the southeastern USA, Asia, and southern South America have relatively higher retranslocation rates (Figure 5). These high N and P retranslocation values are driven by high LAI and the proportional high leaf N and leaf P concentrations simulated by ELM-FUN3.0. Areas with low soil mineral N and low soil mineral P have high retranslocation rates and areas with low leaf N and low leaf P, such as the tundra and arid/semi-arid regions, have the lowest retranslocation rates (Figures 5c and6b). Moreover, the retranslocation efficiency values for N and P, that is, the ratio of the retranslocated nutrient to the amount of nutrient in dead leaves prior to senescence, varies from less than 15% in grassland areas in China and southern South America, to up to 50% in the boreal regions for N and P, and savannas for N, which is in agreement with observations (Aerts, 1996). The tropics have a highly variable N retranslocation efficiency ranging from around 30% in the broadleaf evergreen forests to 50% over grasslands (Figure 7a), and a P retranslocation efficiency of around 40% (Figure 7b), which seems slightly underestimated in comparison to observational studies that indicate 50% (Aerts, 1996;Reed et al., 2012).

Similarly, both the rate of P uptake and the pathway of P acquisition varied as a function of ecosystem productivity and dominant mycorrhizal type globally. ELM-FUN3.0 simulated the mean annual global total P uptake of 48.1 Tg P yr -1 , compared to that simulated by E3SMv1.1-CTC of 42.0 Tg P yr -1 and by E3SMv1.1-ECA (Equilibrium Chemistry Approximation) of 63.0 Tg P yr -1 (Burrows et al., 2020). The spatial patterns in the global distribution of the three major P uptake pathways is shown in Figure 6. The mycorrhizal root uptake and the direct root uptake (i.e., the non-mycorrhizal root uptake) follow the spatial distribution of the most productive areas of the globe, but with significantly different magnitudes (Figure 6). The global total mycorrhizal root uptake is 20.0 Tg P yr -1 , due largely to AM, which is 55.4% of the total mycorrhizal root uptake (11.1 Tg P yr -1 ) (Figure 6), with 44.6% of the remaining from ECM (8.9 Tg P yr -1 ) (Figure 6); direct root uptake is 20.9 Tg P yr -1 . This is due to the differences in the C cost of uptake by mycorrhizal roots compared to non-mycorrhizal roots as defined in Allen et al. (2020) and the spatial distribution of soil available P (Yang, Thornton, et al., 2014, Figure 5. The global mean annual (1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005) (a) symbiotic biological N fixation, (b) direct root uptake (i.e., non-mycorrhizal root uptake; direct root uptake is called in the figures for short), (c) retranslocation, (d) arbuscular mycorrhizal (AM)-associated uptake, and (e) ectomycorrhizal (ECM)-associated N uptake (g N m -2 yr -1 ) simulated by ELM-FUN3.0. Yang, Post, et al., 2014) (Figure S3 in Supporting Information S1). Most of the tropics, southeast South America, eastern USA, and southeast China have the largest values of soil N and P uptake globally. The boreal forests and tundra ecosystems have lower rates of P uptake in comparison to the tropics, with higher values associated with mycorrhizal uptake and retranslocation. Arid and semi-arid areas have the lowest P uptake (Figure 6).

ELM-FUN3.0 estimated that 8% of NPP or 2.5 Pg C yr -1 was spent by plants to take up N. The two largest amounts of C were spent on mycorrhizal root uptake (i.e., 1891.7 Tg C yr -1 or 74.9%) and on direct root uptake (i.e., 623.0 Tg C yr -1 or 24.7%), respectively. The C spent on retranslocation was 20.6 Tg C yr -1 , which corresponds to 0.42% of the total used C. The C spent on N fixation was 0.8 Tg C yr -1 , which corresponds to 0.03% of the total used C. ELM-FUN3.0 estimated that 5% of NPP or 1.6 Pg C yr -1 was spent by plants to acquire P. The two largest C expenditure fluxes were direct root uptake (i.e., 859.6 Tg C yr -1 or 52%) and mycorrhizal root uptake (i.e., 781.6 Tg C yr -1 or 47%), respectively. C spent on retranslocation was 4.2 Tg C yr -1 , which corresponds to 0.26% of the total used C.

Direct root uptake and mycorrhizal uptake were the two pathways where the majority of C expenditure was predicted. However, the expenditure of C is heavily dependent on PFT (Figure 8). Grasslands had the highest C use rate for P acquisition (22.8 ± 6.2 g C m 2 yr -1 ) with 32% of the total C spent on P acquisition used to support retranslocation (Figure 8). Shrublands showed similar total C expenditure to deciduous broadleaf forests, with 38% and 40% of the total C spent on P acquisition used to support mycorrhizal root uptake, respectively. Evergreen broadleaf and evergreen needleleaf forests had a total C use rate of 3.3 ± 0.5 g C m 2 yr -1 and 2.1 ± 0.3 g C m 2 yr -1 , respectively. In terms of the main uptake pathway, mycorrhizae represent 49% and 55% of the total C spent on P acquisition in evergreen broadleaf and evergreen needleleaf forests, respectively. In deciduous needleleaf forests, retranslocation represents 89% of the total 0.4 ± 0.4 g C m 2 yr -1 spent on P acquisition.

One metric of global N and P limitation is the ratio of leaf N and P to leaf N and P retranslocation-the stoichiometric homeostasis theory (Sterner & Elser, 2002), in which the demand ratio of N versus P for a mature leaf should be relatively constant to maintain leaf functioning (Du et al., 2020). ELM-FUN3.0 indicates a stronger P limitation in the tropics (Figure 9), except over grasslands, which are more N limited (see Figure S4 in Supporting Information S1).

In FUN3.0, leaf N and P retranslocation are related to the environmental cost of N and P by the optimization of C expenditure on retranslocation as the nutrients become more expensive with extraction. Estimates here indicate that 6.1% of the natural terrestrial land area is limited primarily by N, whereas 13.9% is primarily P limited. The remaining 80.0% of the land area is co-limited by N and P or weakly limited by either nutrient alone. Although there is a difference in absolute terms, these results are in alignment with those proposed in Elser et al. (2007) and Du et al. (2020) (Figure S9 in Supporting Information S1), in which a larger part of the global terrestrial land area is more relatively P limited and less significantly limited by N.

Much of the uncertainty in ESMs projections of the future land C uptake is driven by how these models represent nutrient constraints on NPP (Ciais et al., 2013;Davies-Barnard et al., 2020;Fernández-Martínez et al., 2014;Fisher et al., 2012;Terrer et al., 2019;Wieder et al., 2015a). Using a process model of the dynamics and C economics of N and P acquisition from the environment coupled into ELM addresses part of this uncertainty by dynamically predicting the C costs of N and P acquisition at the same time, noting that book-keeping approaches may provide unreliable results (Sun et al., 2017). Specifically, we evaluated the total amounts and the spatial distributions of N and P acquisition (Figures 5 and6); the amounts and spatial distributions of C used for N Figure 9. Nearly 80% of the land surface is estimated to be co-limited (gray) by nitrogen (N) and phosphorus (P), with the remaining 20% predominantly limited by either N (red) or P (blue). Global mapping of ELM-FUN3.0 predicted NP limitation for the period 1994-2005 following the stoichiometric homeostasis theory. Global N and P limitations were calculated using the ratio of leaf N and P over the ratio of leaf N and P retranslocation following the stoichiometric homeostasis theory (Sterner & Elser, 2002), in which the demand ratio of N versus P for a mature leaf should be relatively constant to maintain leaf functioning. Estimates indicate that 6.1% of the natural terrestrial land area is limited primarily by N, whereas 13.9% is primarily P limited. The remaining 80% of the land area is co-limited by N and P or weakly limited by either nutrient alone. and P acquisition (Figure 8 and Figure S5 in Supporting Information S1); and changes in NPP as a result of N and P acquisition and allocation at the global scale (Figure 4). ELM-FUN3.0 is able to simulate the distribution of limitation comparably with existing (but limited) empirical estimates (Figure 9 and Figure S9 in Supporting Information S1), as well as the spatial distribution of the C used for acquiring these nutrients (Figure 4b), which provide an important dynamic constraint on predictions of the future land C uptake. In parallel to the addition of new processes, ESMs benefit from systematic benchmarking against observed datasets throughout model development, within initiatives like ILAMB (Collier et al., 2018). Nonetheless, global nutrient cycling benchmarks are limited (e.g., Sun et al. (2020)) and more upscaling exercises from in situ observations are needed (e.g., Du et al. (2020)), as well as the implementation of new technologies to remotely detect information in the rhizosphere (Fisher et al., 2016;Sousa et al., 2021). In this study, we benchmarked all the ELM-FUN model versions (Figures 2,3 and4) indicating that considering the C costs of N and P acquisition are important for improving the representation of the C cycle.

Accounting for N and N-P limitation lowered historical estimates of NPP by 20% and 50%, respectively (Figure 4a). Our result aligns with a previous study using CMIP5 models for N (Wieder et al., 2015a), but shows a much stronger effect of P on NPP (twice as much). Although previous authors have reported the unlikelihood of the land actually becoming a net source of CO 2 (Brovkin & Goll, 2015), mainly due to changes in mineralization and additional processes governing the recycling of P, our simulations explicitly represent P dynamics and indicate a significant decrease in land C uptake. Similar results were reported for the Amazon, showing that P availability reduces the projected CO 2 -induced biomass growth by 50% compared to estimates from C models (Fleischer et al., 2019), although inconclusive based on limited observations (Turner et al., 2018). ELM-FUN3.0 is on the low end of the CMIP6 NPP range, given that part of the available C is spent on P acquisition. ELM presents a positive GPP bias in comparison to FLUXCOM (0.65 g m -2 day -1 ), also reflected in LAI compared to MODIS (0.95 m 2 . m -2 ), while the model versions with FUN work to reduce the bias, especially in the tropics. Nevertheless, ELM-FUN3.0 reduces tropical and arctic-boreal LAI turning the positive bias into a negative one (-0.17 m 2 . m -2 ). The cumulative Net Biome Production (NBP, Figure S11a in Supporting Information S1) highlights that although none of the model versions are able to capture the C source strengthening up until 1940, and weakening afterwards, all versions are contained within the envelope of observational uncertainty after the 60's.

It is likely that mycorrhizal associations can explain a large fraction of the variance in plant response to elevated CO 2 (Braghiere, Fisher, et al., 2021;Drake et al., 2011;Kivlin et al., 2013;Orwin et al., 2011;Sulman et al., 2017;Terrer et al., 2016Terrer et al., , 2018)), because it has been estimated that fungal root symbionts may be responsible for as much as 80% of plant N and P uptake (van der Heijden et al., 2015). Nevertheless, the large-scale spatiotemporal distribution of these plant-fungi associations are highly uncertain (Norby et al., 2017;Sulman et al., 2019). However, in order to enable global simulations with ELM-FUN3.0, an explicit representation of the functional differences between different types of plant symbiotic associations was assumed following the PFT distribution in the original ELM (Table S3 in Supporting Information S1).

The rationale behind this assumption is based upon previously documented geographic distributions of plant symbiosis (Allen et al., 1995;Read, 1991;Shi et al., 2016), where a transition from AM to ECM dominance was reported at single site level observations following a common shift from P to N limitation along increasing latitudes (McGroddy et al., 2004;Reich & Oleksyn, 2004). Although until very recently only a few regional maps of present and past plant symbiotic status were available (Brundrett, 2017;Fisher et al., 2016;Jo et al., 2019;Menzel et al., 2016;Swaty et al., 2016), in the last couple of years a number of different methods were used for extrapolating spatially sparse observations into explicit global maps suitable for applications within ESMs (Braghiere, Fisher, et al., 2021;Soudzilovskaia et al., 2019;Steidinger et al., 2019;Sulman et al., 2019). Braghiere, Fisher, et al. (2021) show that the map presented in Steidinger et al. (2019) follows the spatial distribution proposed in this study (tropics-AM dominated vs. extra-tropics-ECM dominated), indicating that climate variables are the main drivers of global biogeography of plant-fungi symbioses and that fixed values of mycorrhizal associations can be prescribed following spatial distributions of PFTs in ELM.

Although the C cost of mycorrhizal root uptake of N and P from soil is calculated as a function of availability of soil N and P, and the amount of fine root biomass that can access the nutrient, as described in Equation S3 in Supporting Information S1, currently mycorrhizae do not directly impact availability of N or P, or their interaction with soil organic matter. Further model developments should also include some of these important mechanisms.

All nutrient-vegetation models are predicated on the correct functioning of numerous interacting nutrient fluxes in both the soil and vegetation. Here we focused on vegetation processes. The ELM model originated from CLM4.5 with a vertically resolved soil decomposition cascade, which is known for underestimating the production of mineral N resulting in a highly N-limited terrestrial C sink (Koven et al., 2013). In fact, the modeled soil biogeochemistry in ELM-CTC is more sensitive to soil N dynamics than to vegetation biogeochemistry (Tanga & Riley, 2018). Likewise, although ELM considers major P cycle processes, including phosphatase activity, which can relieve P limitation with increasing P demand under elevated CO 2 , there are large uncertainties in the availability and dynamics of soil P (Yang et al., 2019). Thus, the soil mineral N and P uptake simulated by ELM-FUN3.0 is likely underestimated. As such, the apparent nutrient limitation down-regulation on NPP we attribute to FUN is in part attributable to the underestimated soil pools already in ELM. However, coupling FUN3.0 with ELM in E3SM opens up new possibilities to explore different acquisition pathways that could alleviate nutrient limitation.

For example, Richter et al. (2006) suggested that occluded P can become available to plants through mycorrhizal associations. The new global implementation of mycorrhizal associations in ELM-FUN3.0 would allow the development and exploration of this hypothesis. Also, Helfenstein et al. (2020) suggests current LSMs largely overestimate mean residence time of inorganic P in soils, as well as the temporal dynamics of P pools vary largely between different soil types, which could explain part of the uncertainty related to soil P. The active role of plants in changing allocation and flexible stoichiometry are not included in these simulations, but could help alleviate P limitation under elevated CO 2 (Buendía et al., 2014). Likewise, new model developments should include environmental mechanisms controlling phosphatase activity (Alewell et al., 2020), alternative methods for plants to access occluded P pools (Schubert et al., 2020), and organic N (Näsholm et al., 2009), such as rhizosphere priming, defined as the increment or suppression of soil organic matter decomposition by live roots and associated rhizosphere microorganisms (Cheng et al., 2014). Further studies are needed to estimate and clarify the relative role of P limitation and acquisition strategies, which can also be further explored with ELM-FUN3.0.

Many CMIP6 ESMs now consider nutrient cycling (Arora et al., 2020;Jones et al., 2016). Terrestrial C budgets, as well as climate feedbacks, are different from C-only versions because plant C production and allocation are strongly limited by N or P (or both) (Zaehle et al., 2015). Although the newly developed ELM-FUN3.0 expands the abilities of determining interactions in the rhizosphere, including mycorrhizal and direct root uptake, as well as retranslocation of N and P from leaves at the same time, there is room for model development and improvement.

Additional complexity in ESMs, while often justifiable from basic process understanding, has typically been accompanied by an ever-increasing number of parameters (Fisher & Koven, 2020), many of which are poorly constrained by observations (Braghiere, Wang, et al., 2021). New processes can also hinder transparency, robustness, and predictive power due to accumulation of uncertainty (Bonan & Doney, 2018;Braghiere et al., 2019;Franklin et al., 2020;Prentice et al., 2015). In regards to this, there are a few steps to be taken in future ESM development, particularly when referring to C-N-P dynamics in these models, as well as uncertainty of the parameters used in FUN3.0 itself.

While ELM-FUN adds on multiple new model capabilities, inherent model uncertainties encompass parametric and structural uncertainty, which dominate the uncertainty in the land C cycle, while uncertainty in the ocean C cycle is more dominated by data forcing (Lovenduski & Bonan, 2017). Parametric uncertainty in LSMs has been traditionally explored through experimentation with different parameter values in order to test different hypothesis on how these parameters impact predictions (Dagon et al., 2020;Hawkins et al., 2019;Huo et al., 2019;Ricciuto et al., 2018). Originally in FUN, the cost functions (see Equation 3.0 in Supporting Information S1) are a set of parameters that determine the environmental cost of nutrient uptake from soils and retranslocation from leaves. For instance, root biomass does not necessarily need to be high if nutrient concentration is high, given non-limitation of other factors (Aerts et al., 1991;Fisher et al., 2010). Therefore, the cost parameters are defined to control the balance to which available soil nutrient and fine root biomass impact the C cost uptake (Allen et al., 2020).

For instance, FUN simulates the ability of AM fungi to act as scavengers by reducing the AM C cost parameter relative to that of ECM fungi for N and increase for P to reflect the greater efficiency of different C-intensive strategies (Raven et al., 2018). Likewise, the AM and ECM cost parameters were chosen so that the thresholds in nutrient availability and root biomass mirror empirical shifts in the abundance of these two types of fungi across fertility and latitudinal gradients (Allen et al., 1995;Brzostek et al., 2014;Phillips et al., 2013). This abstraction, however, means that uncertainties associated with the cost functions are large (Fisher et al., 2019). Fisher et al. (2019) assessed the sensitivity of CLM5 C and N cycles to the FUN cost parameters by conducting one-at-a-time perturbations of the parameters around the default state at four flux tower sites, varying the cost parameters an order of magnitude in each direction from the default value, since reasonable ranges of these parameters are unknown and their physiological controls are poorly understood. For CO 2 and N enrichment experiments, the authors found moderate decline in LAI, NPP, and GPP when cost parameters increased, but highlighted that they do not appear to be a first-order control parameter when compared to allocation parameters in CLM5. Similarly, low values of C costs parameters appeared to reduce the responsiveness of the system to CO 2 fertilization. The cost parameters seem to have a direct impact on fertilization responses illustrating the sensitivity of the model to relatively unconstrained processes. Measuring and modeling these parameters remains an outstanding challenge (Wieder et al., 2021).

Similarly, Braghiere, Fisher, et al. (2021) assessed the sensitivity of CLM5 C and N cycles to different spatial representations of mycorrhizal distributions globally and found differences of up to 8 TgN.yr -1 between different maps. It is possible that cost parameters vary following different functions depending on multiple variables other than only latitude and nutrient availability, such as temperature and moisture, species composition, and soils characteristics. Further research is needed in this regard in order to characterize uncertainty and establish a more process-based understanding of parametric controls in FUN.

Although the addition of FUN3.0 into ELM allows new capabilities for simulating global P uptake, including mycorrhizal dynamics, there are a number of processes that are still not represented in the model, such as P retranslocation from all senescing tissues, not only from leaves, soil biological traits and trade-offs to optimize C allocation for P acquisition, such as capturing microbial species composition or changes in soil pH in response to P availability, as well as fixed allocation and resorption of P under soils with limited nutrient availability (Jiang et al., 2019). As incorporation of P dynamics into ESMs develops (Goll et al., 2012;Thum et al., 2019;Wang et al., 2007;Yang, Thornton, et al., 2014, Yang, Post, et al., 2014;Zhu et al., 2019), these models are paving the way for increased inclusion of different P cycle processes (Reed et al., 2015). While large uncertainties in soil P dynamics are due to organic P, inorganic P is a key player in P dynamics (Helfenstein et al., 2020), and a separation of leaf P into organic and inorganic pools would be desirable to improve predictions of photosynthesis (Jiang et al., 2019). Integrating plant parametric adjustments to nutrient use efficiency, such as updating the maximum rate of Rubisco carboxylase activity or the maximum rate of photosynthetic electron transport (Braghiere et al., 2020;Walker et al., 2014), C-N-P ratios of newly formed soil organic matter during decomposition, priming effects on C, N, and P mineralization, and consideration of flexible plant tissue stoichiometry are among the high priority processes needed to be incorporated in further developments of ELM-FUN3.0. The integration of new processes involving C-N-P dynamics can now feasibly be added into the cost structure of FUN3.0, given its highly modular nature (e.g., the cost functions) based on optimality approaches making use of economic concepts, recently advocated as one of the organizing principles to be used in next generation ESMs (Fisher & Koven, 2020;Franklin et al., 2020).

Lastly, the response of N cycling to C fertilization is to some degree a function of the abundance of N fixing plants (Fisher et al., 2019). In ELM-FUN3.0 this is a fixed entity (25% for natural ecosystems). Increasing use of models that resolve vegetation demographics (Koven et al., 2020;Longo et al., 2019) should allow prognostic estimates of how N fixing strategies vary in high CO 2 conditions (Meyerholt et al., 2016).

Using a process model of the dynamics and C economics of N and P acquisition from the environment coupled into the E3SM land model, we evaluated the total amounts and the spatial distributions of N and P acquisition; the amounts and spatial distributions of C used for N and P acquisition; and changes in NPP as a result of N and P acquisition and allocation at the global scale. We benchmarked all the ELM-FUN model versions indicating that considering the C costs of N and P acquisition are important for improving the representation of the C and water cycles, and that accounting for N and N-P limitation lowered historical estimates of NPP by 20% and 50%, respectively. Modeled and observed nutrient limitation agreement increases when N and P are considered together (r 2 from 0.73 to 0.83), but the likely nutrient limitation down-regulation on NPP we attribute to FUN is in part attributable to the underestimated soil nutrient pools in ELM.