Ensemble projections elucidate effects of uncertainty in terrestrial nitrogen limitation on future carbon uptake

The magnitude of the nitrogen (N) limitation of terrestrial carbon (C) storage over the 21st century is highly uncertain because of the complex interactions between the terrestrial C and N cycles. We use an ensemble approach to quantify and attribute process‐level uncertainty in C‐cycle projections by analysing a 30‐member ensemble representing published alternative representations of key N cycle processes (stoichiometry, biological nitrogen fixation (BNF) and ecosystem N losses) within the framework of one terrestrial biosphere model. Despite large differences in the simulated present‐day N cycle, primarily affecting simulated productivity north of 40°N, ensemble members generally conform with global C‐cycle benchmarks for present‐day conditions. Ensemble projections for two representative concentration pathways (RCP 2.6 and RCP 8.5) show that the increase in land C storage due to CO2 fertilization is reduced by 24 ± 15% due to N constraints, whereas terrestrial C losses associated with climate change are attenuated by 19 ± 20%. As a result, N cycling reduces projected land C uptake for the years 2006–2099 by 19% (37% decrease to 3% increase) for RCP 2.6, and by 21% (40% decrease to 9% increase) for RCP 8.5. Most of the ensemble spread results from uncertainty in temperate and boreal forests, and is dominated by uncertainty in BNF (10% decrease to 50% increase for RCP 2.6, 5% decrease to 100% increase for RCP 8.5). However, choices about the flexibility of ecosystem C:N ratios and processes controlling ecosystem N losses regionally also play important roles. The findings of this study demonstrate clearly the need for an ensemble approach to quantify likely future terrestrial C–N cycle trajectories. Present‐day C‐cycle observations only weakly constrain the future ensemble spread, highlighting the need for better observational constraints on large‐scale N cycling, and N cycle process responses to global change.


| INTRODUC TI ON
Over the last decades, the terrestrial biosphere has sequestered roughly a quarter of the CO 2 emitted by anthropogenic activities (Le Quéré et al., 2018). The future evolution of this carbon (C) uptake depends, among others, on the availability of nitrogen (N) to support this increase in terrestrial C (Hungate, Dukes, Shaw, Luo, & Field, 2003). Nitrogen is an essential nutrient for life, and therefore plant growth and soil microbial activity, but its availability is limited in most natural terrestrial ecosystems Hyvönen et al., 2007;LeBauer & Treseder, 2008;Vitousek & Howarth, 1991).
Terrestrial biosphere models (TBMs) that estimate the future C-uptake potential of land vegetation increasingly consider the dynamics of the global N cycle and its linkage to the C cycle (Arora et al., 2019;Zaehle & Dalmonech, 2011). These models generally suggest that N constraints attenuate the land C response to global change (Sokolov et al., 2008;Thornton et al., 2009;Zhang, Wang, Matear, Pitman, & Dai, 2014). This is also the case in the new set of ESMs in the Coupled Model Intercomparison Project Phase 6 (CMIP6; Arora et al., 2019), where about half of the models include a N cycle representation. However, adequate ecosystem-scale characterization of terrestrial N cycle processes remains a challenge because of the complexity of N cycle processes, their spatial heterogeneity and the long-term cumulative effect of comparatively small fluxes such as biological nitrogen fixation (BNF) on ecosystem N availability (Thomas, Brookshire, & Gerber, 2015). As a consequence, the inclusion of N cycling in complex ESMs has been associated with considerable uncertainty, making it difficult to quantify the effect of N limitation in simulations of the biospheric C sink response to future global change (Zaehle & Dalmonech, 2011).
A better understanding of terrestrial N limitation in the Earth system requires the identification of the most important N-process mechanisms (Thomas et al., 2015). In the past, model comparisons that aimed to characterize the contribution of such N-related model uncertainties to model differences in future predictions were always confounded by alternative model structures not associated with the N cycle (Fleischer et al., 2019;Huntzinger et al., 2017;Medlyn et al., 2016;Thomas, Zaehle, Templer, & Goodale, 2013;Zaehle et al., 2014). Du et al. (2018) have recently highlighted the large difference in the simulated steady-state terrestrial N cycle by comparing selected components of the N cycle representation of three biosphere models within a common modelling environment. Here, we use an alternative approach that relies on the implementation and factorial combination of published strategies to represent N cycling in TBMs (Meyerholt & Zaehle, 2015Meyerholt, Zaehle, & Smith, 2016; Table 1) into one common framework to provide a robust assessment of the likely range of future terrestrial C storage and attribution to underlying model assumptions.
We focus on the representation of processes that affect the longterm trajectories of the coupled land C and N cycles in response to changing atmospheric CO 2 , climate and N deposition: processes that affect the ratio of C to N storage in the land biosphere, as well as the balance of ecosystem N in-and outputs (Gruber & Galloway, 2008;Hungate et al., 2003;Walker et al., 2015;Zaehle et al., 2014). Alternative assumptions about the dependence of BNF and N loss on plant N demand can lead to diverging responses of the terrestrial N balance to global changes, and therefore change the amount of N available to store C (Du et al., 2018;Meyerholt & Zaehle, 2018;Meyerholt et al., 2016;Wieder, Cleveland, Lawrence, & Bonan, 2015). Alternative assumptions about the flexibility of the plant and soil organic matter (SOM) C:N ratios affect the response of tissue-specific process rates (such as photosynthesis and respiration), the relative competitive strength of plants and SOM for N, and finally the amount of C that can be stored given a specific N amount (Meyerholt & Zaehle, 2015;Thomas et al., 2013).
We selected two C:N stoichiometry, five BNF and three N-loss algorithms, and used these to generate an ensemble of 30 factorial combinations of these algorithms within the O-CN C-N cycle model (Zaehle, Ciais, Friend, & Prieur, 2011; see Section 2). Here, we first analyse the C-N cycle ensemble with respect to its ability to simulate the global contemporary C and N cycle commensurate with observations. We then analyse the trajectories resulting from historical  and projected  changes of atmospheric CO 2 , climate, and N deposition for two global change scenarios, specifically the representative concentration pathway (RCP) 2.6 (with an average global land air temperature (GLAT, but without Antarctica) increase of 2.4°C between 1850 and 2099) and 8.5 (GLAT increase of 8.0°C; Dufresne et al., 2013;Hempel, Frieler, Warszawski, Schewe, & Piontek, 2013;Meinshausen et al., 2011), and compare these to the C-only reference version of O-CN. We attribute model uncertainty to underlying process formulations and investigate to what extend a range of contemporary C-cycle benchmarks constrain the model spread. Based on this analysis, we estimate the anthropogenic emissions compatible with each RCP scenario, given the RCP-specific atmospheric CO 2 trajectory, simulated oceanic uptake (Dufresne et al., 2013), as well as the projected land C uptake considering N constraints.

| Model description
Model simulations were performed with the O-CN model , an extension of the ORCHIDEE model (Krinner et al., 2005), the land surface component of the Institute Pierre Simon Laplace (IPSL) climate model (Dufresne et al., 2013). O-CN describes a fully prognostic N cycle that accounts for ecosystem N input from atmospheric deposition and biological fixation, and N losses due to leaching of dissolved organic and inorganic N, as well as gaseous losses associated with volatilization, nitrification and denitrification. Once soil mineral N is taken up by plants, it cycles as organic N between different plant, litter and SOM pools, where its tissue-specific concentration affects among others plant photosynthesis and respiration, plant allocation to fine roots and vegetation growth. Soil mineral N availability also affects litter and SOM decomposition rates, as well as the rate of net N mineralization, which depends further on the difference in C:N stoichiometry of decomposing litter and SOM pools. Overall, the O-CN model has shown good performance against a range of global biosphere benchmarks (Le Q uéré et al., 2018;. To isolate the effect of N-cycle parameterization from other model structural assumptions of O-CN, we generated a C-only version of O-CN. This version assumes for each plant tissue-type or soil biogeochemical pool, vegetation-type specific, time-invariant N concentrations constrained by observed estimates (Kattge et al., 2011). In this model version, N is added to the plant labile N pool (soil mineral NH 4 pool) whenever limiting N availability would reduce plant growth (SOM decomposition). As a result, this model version responds to changes in atmospheric CO 2 , temperature, precipitation and so on as the model without a full N cycle would. In particular, the response of terrestrial C gain to CO 2 is driven by the biochemical response embedded in the leaf-photosynthesis model (Kull & Kruijt, 1998), and only modulated due to the scaling from leaf to canopy, and the response of vegetation dynamics and turnover to changes in production . In a similar manner, other ecosystem processes, such as autotrophic and heterotrophic respiration respond to temperature and soil moisture, and are unaffected by N availability .
For the purpose of this study, the N cycle representation in O-CN was expanded to include alternative algorithms to represent the key N mechanisms of C:N ratio flexibility, BNF and N loss (Table 1; Sections 2.1.1-2.1.3). Each of the added algorithms was previously published and applied in N-enabled land C cycle models. The details of the equations, parameters and the implementation into the O-CN code are described by Meyerholt and Zaehle (2015; C:N flexibility), Meyerholt et al. (2016;BNF) and Meyerholt and Zaehle (2018;N losses). We combined these algorithms in a factorial manner to generate an ensemble of 30 C-N models (O-CN revision 295, see Meyerholt, Sickel, & Zaehle, 2019). TA B L E 1 Overview on the alternative nitrogen cycle algorithms employed in this study to represent carbon:nitrogen ratio (C:N) flexibility, biological nitrogen fixation (BNF) and ecosystem-level nitrogen losses. All algorithms, including their parameterizations are described in full in Meyerholt and Zaehle (2015;M1;C:N flexibility), Meyerholt et al. (2016;M2;BNF) and Meyerholt and Zaehle (2018;M3; nitrogen losses)

Model
Description Reference Documentation C:N flexibility, see Section 2.1.1

R2
C:N in organic plant and soil pools flexible within prescribed, PFT-specific bounds  FLX in M1 Biological nitrogen fixation (BNF), see Section 2.1.2

F1
Linear relationship with time-invariant climatology of actual evapotranspiration, set to zero above an soil inorganic N pool of 2 g N/m 2 Cleveland et al. (1999),  FOR in M2

F2
Monotonically increasing, saturating function of the simulated time-variant annual net primary production Cleveland et al. (1999), Goll et al. (2012), Thornton et al. (2009)  Carbon:nitrogen ratios influence, for instance, leaf-specific photosynthetic efficiency through its dependence on foliar N content (Field & Mooney, 1986). They also affect the net N mineralization rate through its dependence on the difference between the C:N ratios of plant litter and SOM (e.g. Manzoni, Trofymow, Jackson, & Porporato, 2010). Finally, they affect the overall amount of N required for C storage (Hungate et al., 2003). Some models represent C:N ratios as time-invariant, tissue type-dependent constants, assuming that the observed difference of, for example, needleleaved and broad-leaved foliage C:N ratio result from a tissue typespecific stoichiometric homeostasis of leaf growth and N supply (e.g. Sokolov et al., 2008;Thornton et al., 2009; C:N algorithm R1).
The implications of flexible stoichiometry with increased N availability are on the one hand increased biochemical rates (photosynthesis, respiration, gross N mineralization and immobilization), and on the other hand decreased the ability to sequester C per unit plant or soil N.

| Biological nitrogen fixation
Biological nitrogen fixation through asymbiotic or symbiotic pathways is the most important form of N input in many natural ecosystems (Vitousek et al., 2002). Reflecting the lack of consensus about the mechanisms that control ecosystem-level BNF, there is a large range of published approaches ranging from a global, timeinvariant map of BNF to process-oriented formulations that allow BNF to adjust to ecosystem N requirements. This choice does not only affect the geographic patterns of BNF, but also the response to ecosystem productivity, and therefore both the N capital of ecosystems in equilibrium and their response to global change. The simplest strategy, BNF algorithm F1, assumes that the large-scale pattern of BNF corresponds to patterns in long-term average climate conditions (Cleveland et al., 1999), but does not respond to environmental change at timescales considered in land surface models (as used e.g. in . BNF algorithm F2 relates BNF to net primary productivity as a saturating function (e.g. Thornton, Lamarque, Rosenbloom, & Mahowald, 2007), which, while derived from the study of Cleveland et al. (1999) has limited support from the data presented in that paper. BNF algorithm F3 (e.g. Goll et al., 2012;Smith et al., 2014), based on correlating BNF with simulated evapotranspiration (as in Cleveland et al., 1999), suggests that BNF increases with foliar coverage and productivity, but is reduced by water limitation or increasing atmospheric CO 2 . In contrast to these phenomenological strategies, BNF algorithm F4 formulates BNF as a time-lagged response to vegetation N demand (or deficit to sustain growth), which furthermore considers a light limitation factor for non-tropical plants, reducing BNF in closed-canopy ecosystems (Gerber, Hedin, Oppenheimer, Pacala, & Shevliakova, 2010).
BNF algorithm F5 is based on the principle of optimality of resource investment into nutrient acquisition, where BNF is regulated by the relative C cost of N uptake through BNF or root uptake, indirectly responding to plant N demand and supply (Meyerholt et al., 2016;Rastetter et al., 2001).

| Ecosystem nitrogen losses
Ecosystem nitrogen losses occur due to leaching of inorganic N or gaseous losses during volatilization, nitrification and denitrification (Gruber & Galloway, 2008). Some models simulate these pro- loss algorithm L2) is to hypothesize that the majority of N loss is directly related to soil N turnover, that is gaseous losses are calculated as a fixed fraction of the mineralization flux. In these models, plant N demand has limited control on N-cycle openness and therefore plant N availability. An intermediate strategy (Thornton et al., 2007), described by loss algorithm L3, applies a hierarchical structure, removing a fixed fraction of excess N left by soil immobilization, plant uptake and turnover-based denitrification losses to represent volatilization of reactive N species, further denitrification and leaching loss. respectively. From 1850, we performed a number of simulations for the 1850-2099 period that varied in the applied transient forcing of atmospheric CO 2 , N deposition  and climate (see also Figure S2a-c). Since the climate data were only available from 1901, for the period 1850-1900 climate of the years were randomly chosen as during the spin-up period (see Table 2).

| Ensemble evaluation
We evaluated the simulated terrestrial productivity from the s3 experiment against two alternative estimates: a direct comparison against independent data-driven estimates based on the upscaling of site-based eddy-covariance measurements (Jung et al., 2011). As an alternative metric, we evaluated the long-term mean  polewards increase in the seasonal cycle amplitude of atmospheric CO 2 . Although this increase reflects the balance of terrestrial productivity and respiration, it is strongly correlated with the largescale latitudinal distribution of terrestrial productivity (Dalmonech & Zaehle, 2013;Heimann et al., 1998). As a simplified metric, we diagnosed the polewards increase in the seasonal cycle by regressing the detrended seasonal amplitude against the sine of latitude at the 13 long-term monitoring stations (Ω lat , ppm/sin(lat)).
For comparison to atmospheric CO 2 measurements, the simulated monthly net land-atmosphere CO 2 fluxes for the years 1982 to 2011 from our 30 + 1 model versions were transported to 13 longterm atmospheric monitoring stations (Dalmonech & Zaehle, 2013; see Table S3 for the list of stations used) using the Jacobian representation of the TM3 atmospheric transport model, version 3.7.22 (Kaminski, Heimann, & Giering, 1999;Rödenbeck, Houweling, Gloor, & Heimann, 2003), with interannually varying wind fields (Kalnay et al., 1996, updated), together with estimates of the net ocean-atmosphere C flux Mikaloff Fletcher et al., 2006;, as well as estimated global fossil fuel emission fluxes (Boden, Marland, & Andres, 2013). We further evaluated the sensitivity of the atmospheric CO 2 growth rate to interannual climate variability (γ IAV , ppm/K) by correlating the observed or simulated interannual CO 2 growth rate

| Attribution of model spread to a particular algorithm
In order to characterize the effect of a particular model algorithm on the ensemble mean, we grouped ensemble members according to the algorithm used, for example model group R1 includes all ensemble members that employ the C:N algorithm R1. We then calculated the mean for this group (M group ) and compared it to the mean of the entire ensemble (M ensemble ). For comparability across different model outputs and regions, we calculated a normalized z-score as where SD ensemble is the standard deviation of the entire ensemble. Note that the number of ensemble members (n) for the calculation of M group differs across the process groups of C:N (n = 15), BNF (n = 6) and N loss algorithms (n = 10; see Table 1).

| Carbon-cycle sensitivities
For a better understanding of the N-cycling effect on projected land C storage (vegetation, litter and soil C) in the year 2099, we decomposed the change in land C into C-cycle sensitivities to changes in atmospheric CO 2 (β L , Pg C/ppm) and climate change (γ L , Pg C/K), as in Arora et al. (2013). Because our simulations consider transient changes in N deposition, an additional term to describe the carbonnitrogen sensitivity (η L , Pg C/Pg N), that is the change in land C in response to the cumulated amount of atmospheric N deposition was added to the framework.
where ΔC L is the simulated change in land C storage, C A is atmospheric CO 2 (ppm), T is the global 30-year mean land air temperature (excluding Antarctica) and ∑ N is the cumulative N deposition on land, all between 1850 and 2099. These C-cycle sensitivities were derived from the factorial simulations s0-s3 as follows: The C-concentration sensitivity (β L = ΔC L /ΔC A ) was derived from the s1 to s0 difference in land C storage in the year 2099 and solely reflects the effects of atmospheric | 3983 CO 2 . The C-N sensitivity ( L = ΔC L ∕ ∑ N) was derived from the s2 to s1

TA B L E 2 Overview of the simulations performed
difference in land C storage in the year 2099, and therefore accounts for potential synergistic interactions between the C-concentration and C-N sensitivities. The C-climate sensitivity ( L = ΔC L ∕ΔT) was derived as the difference between s3 and s2 land C storage in 2099 and therefore accounts for potential synergistic interactions between the C-concentration, C-N and C-climate sensitivities. The comparison of these sensitivities across studies needs to consider that these properties are known to be scenario dependent .
Estimates of the C-cycle response to C and N perturbations at global scale are scarce. In an attempt to provide an indication of plausibility of the simulations, we evaluate the ensemble by comparison of meta-analyses on plant growth responses of elevated CO 2 (Baig, Medlyn, Mercado, & Zaehle, 2015) and N addition experiments (Schulte-Uebbing & de Vries, 2018).

| Implied RCP-compatible emissions
In order to diagnose the N-related imprint on the anthropogenic where k (=2.12 Pg C/ppm) is the conversion factors of atmospheric CO 2 abundance to mass , and i is one particular ensemble member.

| Present-day model performance
First we evaluate the ensemble with respect to its capacity to predict the present-day global C and N cycles in terms of overall stocks, gross fluxes and the net land-atmosphere C flux over the recent period (1996-2005 mean). The magnitude of the simulated net land-atmosphere C flux during the 1990s and 2000s, as well as the increment from the 1990s to the 2020s is within the range 138 Pg C/year for the C:N ensemble and C-only model respectively) and net primary production (NPP: 61 ± 4 Pg C/year vs. 62 Pg C/ year). These numbers are broadly consistent with independent estimates (Beer et al., 2010;Saugier & Roy, 2001). Simulated global BNF ranges from 37 to 117 Tg N/year and thus accounts for most of the large literature range of plausible BNF estimates Vitousek, Menge, Reed, & Cleveland, 2013), although one other study has suggested significantly larger rates (Xu-Ri & Prentice, 2017). The range of global N-use efficiency of vegetation production, that is the ratio of NPP to plant N uptake, is 46-66 g C/g N. A detailed overview of the model estimates is given in the Supporting Information (Table S1).
The limiting effect of N availability on C cycling in the ensemble is most prevalent polewards of 40° latitude, where the C-N ensemble shows lower productivity. For the Northern Hemisphere, this is in general terms in better agreement with independent data-driven estimates of gross primary productivity (Jung et al., 2011) than the C-only reference (Figure 2a,b). This finding is consistent with the comparison to atmospheric CO 2 observations (Figure 2c). Most ensemble members simulate the observed polewards increase of the seasonal cycle amplitude of atmospheric CO 2 concentrations (Ω lat ) better than the C-only reference. All ensemble members simulate the seasonal phasing and amplitude, as well as the sensitivity of atmospheric CO 2 to interannual climate variations (γ IAV ) largely in agreement with observations (Supporting Information, Figure S1 and Table S2).  increases under both RCP scenarios during the first half of the 21st century (Figure 4a,b). The spread then declines again under RCP 2.6 as atmospheric CO 2 levels stabilize, whereas it continues to increase in the RCP 8.5 scenario with continuously growing atmospheric CO 2 levels (compare also Figure S2a). Integrated over the scenario pe-   , as well as the scenario periods (2006-2099) for representative concentration pathway (RCP) 2.6 (blue) and RCP 8.5 (red), with different shading indicating model groups. Note that N deposition increases in RCPs 2.6 and 8.5 contribute to the increased N losses independent of changes in BNF or internal nitrogen cycling. Model groups (n = 15, 6, 10 for the groups of carbon:nitrogen, BNF and N loss models) comprise all models in Table 1. See Figure S2 for time series of these changes CO 2 after the year 2030 in RCP 2.6 causes the model spread in net land-atmosphere CO 2 exchange to slowly decrease, as ensemble members with stronger N limitation gradually acquire the 'missing' N through increased BNF that outpaces N losses.

| Carbon-cycle sensitivities
To further understand the cause of the model spread in projected land C change from 1850 to 2099 under RCP 8.5, we decompose the projections in terms of their sensitivity to atmospheric CO 2 concentration increase, climate change and increased N deposition. The change in land C due to increasing atmospheric CO 2 concentrations, that is the C-concentration sensitivity β L (Figure 6a,b) is the primary cause for the spread in projected land C storage change (Figure 4a). We find that N constraints reduce β L in 29 out There is a tight correlation between the response of NPP to elevated levels of CO 2 and β L (r 2 = .83, p < .01; Figure S3), because in the models, NPP is the primary driver of land C stock changes. An evaluation of the simulated, transient global NPP response to elevated CO 2 is challenging, because most experiments were done in artificial setting, for a short time and following a step-increase. It is nevertheless interesting to note that most  The change of land C due to climate change, that is the C-climate sensitivity (γ L ) averages −20 ± 5 Pg C/K (mean ± SD; Figure 6c,d), corresponding to a 19 ± 20% attenuation of the C losses simulated by the C-only reference (γ L = −25 Pg C/K). This attenuation is primarily the result of N fertilization of vegetation due to enhanced soil N mineralization, which partially compensates increased for soil C losses as a result of increased soil respiration in a warmer world. Such an effect has been observed by one large-scale forest soil warming experiment (Magill et al., 2004), but the number of studies on full ecosystem warming effect on N cycling is too small to corroborate this effects at large scales. The attenuation is strongest with fixed stoichiometry (C:N algorithm R1) and N loss scheme L1. In the other two loss schemes, increased net mineralization directly results in an increase in ecosystem N losses through the assumption that a significant fraction of these losses are coupled to the net N mineralization rate, which prevents a larger fertilization effect.
A correlation between the C-climate sensitivity at centennial (γ L ) timescales and interannual (γ IAV ), apparent in the CMIP5 ESM ensemble (Cox et al., 2013), does not exist in our ensemble (Figure 7).
The γ IAV is similar between the C-N ensemble and the C-only reference, suggesting that the direct effects of interannual temperature variations on C cycle processes drive γ IAV in the ensemble. Indirect effects resulting from temperature-related anomalies in N mineralization and plant uptake appear to play a lesser role at this timescale.
The contribution of C:N and loss algorithms to the uncertainty in γ L clearly indicates that these processes contribute notably to the long-term climate response of the biosphere (Figure 6). This finding suggests that the correlation in the CMIP5 ensemble may be a result of the simplified model structure considered, and may be altered when the long-term effects of nutrient-related soil-vegetation feedbacks on the response of ecosystems to perturbations are taken into account.
In the ensemble, β L is highly anti-correlated with the sensitivity of terrestrial C storage to N deposition η L (r 2 = .89, p < .01; Figure 8).
This correlation reflects the differing degree of N limitation in the ensemble: ensemble members with low present-day N availability and therefore low productivity (Figure 2) are more sensitive to N added from deposition, and at the same time also more limited in their ability to respond to CO 2 fertilization with increased production and C storage. Ensemble members with strong simulated N limitation (low response to elevated CO 2 ) also retain more ecosystem C in a warmer climate through a larger fertilization effect from warming-induced increases in net N mineralization. Consequently, γ L in the C-N ensemble is weakly anti-correlated with β L (r 2 = .56, p < .01) and positively correlated with the sensitivity of land C storage to N deposition η L (r 2 = .65, p < .01).
The correlation between the C cycle sensitivities suggest that there is a compensation of uncertainty in terms of future land C storage in N-enabled models. To assess the magnitude of this effect, we evaluated Equation (3) for the range of β L , γ L and η L in the ensemble, but assuming no correlation between them. This simple propagation of errors of the mean yields a change in total land C between 1850 and 2099 of 382 ± 109 Pg C, compared to the spread in the ensemble (382 ± 53 Pg C). This reduction in the model spread by more than half through the interactions of N cycle F I G U R E 7 Correlation between the land climate sensitivity γ L and the response of the atmospheric CO 2 growth rate to interannual climate variability γ IAV  β L (Pg C/ppm) processes highlights the need to account explicitly for these processes and their effect on C cycling in assessments of future biosphere dynamics.

| RCP-compatible anthropogenic CO 2 emissions
We lastly assess the impact of the projected land C uptake on future anthropogenic CO 2 emission pathways that are compatible with the representative CO 2 concentration pathways of the RCP scenarios.
When the N effect is taken into account, the predicted compatible emissions for the 2006-2099 period are reduced to 342 ± 18 Pg C (mean ± SD) for RCP 2.6 and 1,805 ± 42 Pg C for RCP 8.5, compared to the C-only reference of 371 Pg C for RCP 2.6 and 1,863 Pg C for RCP 8.5 (Figure 9).
The estimates of RCP-compatible emissions are strongly correlated with the present-day latitudinal distribution of productivity (Figure 9), as measured by Ω lat (as defined in Figure 2; RCP 2.6: r 2 = .92; RCP 8.5: r 2 = .91). Using the Ω lat constraint to weight ensemble member predictions (see Supporting Information) slightly decreases the mean estimate and range of RCP-compatible emissions (Figure 9), but the effect is much lower than in a previous study based on C-cycle-only models owing to the larger spread in their ensemble (Booth et al., 2017). For RCP 2.6, our ensemble estimate of RCP-compatible emissions under N constraints is 339 ± 11 Pg C (mean ± SD), a reduction from the C-only estimate by 32 Pg C or 9%. For RCP 8.5, we predict RCP-compatible emissions of 1,796 ± 28 Pg C, a reduction relative to the C-only reference by 67 Pg C or 4%.

| D ISCUSS I ON
In this paper, we analysed future projections of the terrestrial C cycle using an ensemble of 30 C-N cycle models with alternative representations of C-N stoichiometry, BNF and N losses. We find that the ensemble is generally within or at least close to the range of observations with respect to contemporary global C-cycle benchmarks. Our consistent framework combining alternative model process representation elucidates the otherwise hidden consequence of model assumptions and thereby provides a more robust estimate of the N constraints on future C budget compared to previous studies (Wieder, Cleveland, Smith, et al., 2015;Zaehle et al., 2015;Zhang et al., 2014). Based on this analysis, we find that N dynamics reduce the global CO 2 fertilization effect under the RCP 8.5 scenario by 24 ± 15% and reduce the C loss associated with global warming by 19 ± 20%. In combination, the result is a reduction of projected land C storage under the RCP 8.5 scenario of 21% (−9% to +41%). The relative reduction in future land C uptake is slightly lower in the climate change stabilization scenario RCP 2.6, where it averages 19% (−3% to 40%). Integrated assessment models that are used to calculate the compatible emissions based on ocean and land C uptake do currently not account for terrestrial N dynamics.
The results of our study suggest that an additional mitigation effort is required to maintain atmospheric CO 2 at levels prescribed in the RCP scenarios of 9 ± 3% (RCP 2.6) and 4 ± 2% (RCP 8.5).
Previous studies (Sokolov et al., 2008;Thornton et al., 2009; have suggested ranges for the C-concentration sensitivity (β L ) in N-enabled models between 0.6 and 1.2 Pg C/ppm, which includes our best estimate of 0.75 ± 0.13 Pg C/ppm (mean ± SD). A direct comparison of these estimates is challenging as they are known to be scenario dependent .
In terms of the relative reduction of β L compared to a C-only model version, which ranged 50% and 80% (Sokolov et al., 2008;Thornton et al., 2009;Wårlind et al., 2014), we find a smaller reduction due to the inclusion of N cycling (24 ± 15%). Our quantification of the C-climate sensitivity (γ L ) is also consistent with previous studies (Sokolov et al., 2008;Thornton et al., 2009;Wårlind et al., 2014;. In all these models, N interactions partly mitigated the positive C-climate feedback because N mineralization from SOM decomposition fertilizes plant growth in N-limited ecosystems and thus leads to enhanced vegetation C storage. However, none of the ensemble members showed a shift to a negative feedback, as suggested by the CLM4 model (Thornton et al., 2009).

| Process attribution
One advantage of our ensemble approach is that it enables direct attribution of ensemble uncertainty to process formulations (Figures 2-6). The comparison demonstrates that alternative model choices of each of the processes considered have characteristic and important implications for the simulated level of present-day and future N limitation. Our analysis demonstrates further that in combination, these three process formulations capture a wide range of C-N cycle interactions: for instance, assuming a combination of a vegetation N demand-driven BNF (BNF algorithm F4), a flexible stoichiometric constraint on growth (C:N algorithm R2), as well as a hierarchical N loss scheme (L3) leads to a model with prognostic N cycle that has little constraining effect on the decadal to centennial C-cycle projections.
Conversely, assuming tight stoichiometric constraints (C:N algorithm R1), direct dependence of BNF on N-limited plant production or CO 2sensitive evapotranspiration (BNF algorithms F2 and F3 respectively), as well as a concentration-based N loss scheme (L1) results in a comparatively strong N limitation effect, but still does not reproduce the very strong N limitation suggested initially by CLM4 (Thornton et al., 2009). Our regional analysis further reveals that different processes are important in different locations/climates, illustrating the regional nature of simulated N limitation.
In the ensemble, the BNF response to elevated CO 2 is the largest source of uncertainty. This reflects the relatively poor understanding of the processes, and their representation in models, affecting BNF at large scale. Increases in response to CO 2 have been observed in controlled experiments (e.g. between 30% and 62% in a recent meta-analysis by Liang, Q i, Souza, & Luo, 2016), but the responses were low or declining over time due to other limiting factors (Hofmockel & Schlesinger, 2007;Hungate et al., 2003Hungate et al., , 2004. Two BNF schemes stand out in their CO 2 response: coupling BNF to evapotranspiration (BNF algorithm F3) leads to the geographically wide-spread prediction that BNF decreases in response to CO 2 by about 7 Tg N/year or 7% globally, despite increasing N limitation of plant growth. Such a behaviour does not appear to be supported by the available evidence. Conversely, assuming a BNF dependence on instantaneous vegetation N demand and light availability (BNF algorithm F4) leads to a large projected increase of BNF (on average by 126 Tg N/year or 208% between 1850 and 2099) in response to increasing CO 2 , particularly in temperate and boreal regions. While such a demand-driven approach appears to be appealing from a process understanding perspective, its current implementation lacks sufficient representation of other limiting factors for the BNF response. Removing these two BNF algorithms F3 and F4 from the ensemble changes projected 1860-2099 land C storage only slightly, although it reduces the uncertainty in the ensemble: the projected land C storage for RCP 8.5 (RCP2.6) of 384 ± 52 (265 ± 28) Pg C is reduced to 379 ± 39 (263 ± 21). Understanding the global patterns of BNF, as well as its environmental controls in response to changing C and nutrient availability should therefore be a priority in reducing N-related uncertainty (see e.g. Zheng, Zhou, Luo, Zhao, & Mo, 2019).
Furthermore, all algorithms considered in our study have in common that they do not consider community assembly and dynamics as part of the processes that constrain BNF responses and only relate BNF to instantaneous changes in ecosystem N properties. It will likely be necessary to explicitly represent the effects of community dynamics and assembly to represent the timescale of BNF shifts at ecosystem level (Menge, Hedin, & Pacala, 2012;Menge, Levin, & Hedin, 2008;Thomas et al., 2015).
The choice of fixed versus flexible C:N ratios contribute to variations in our estimates of both β L and γ L . Assuming fixed stoichiometry leads to a tighter N constraint and thus a more profound impact on the CO 2 and climate response. With flexible stoichiometry (algorithm R2) increasing N availability increases tissue N concentrations, and vice versa. Increased tissue N concentrations lead on the one hand to higher photosynthesis and respiration rates per unit biomass, and on the other hand a higher N requirement for growth (see Meyerholt & Zaehle, 2015, for a more detailed discussion), which in combination reduces the C storage response to added N. While some flexibility has been observed (Ainsworth & Long, 2005;Magill et al., 2004;McNulty et al., 2005), recent model-data intercomparison studies have suggested that it is overestimated by current model formulations (Meyerholt & Zaehle, 2015;Zaehle et al., 2014), clearly demonstrating the need for a better representation of this process.
However, the impact of stoichiometric flexibility on overall land C sequestration considering all forcings is attenuated because of the compensating effects of climate warming and N deposition on plant N availability, and the large share of low C:N SOM pools with limited stoichiometric flexibility in the total land C storage. Therefore, although important from a process representation point of view, the overall contribution to land C storage trends is small compared to that of BNF in our study.

| 3991
The N loss formulation does not have a clear effect on the biospheric response to elevated CO 2 , but shows a strong impact on the climate sensitivity of the C cycle (γ L , Figure 6c). The assumption that ecosystem N loss is primarily controlled by the rate of soil N turnover (loss algorithms L2 and L3) leads to higher predicted soil N losses in response to warming. In these cases, vegetation has limited potential to act as a N sink, leading to a lower fertilization effect compared to loss algorithm L1 (Figure 6c). In particular, we found that a combination of flexible stoichiometry and the hier- When assessing different loss algorithms, feedbacks in the simulated C and N cycles can lead to equifinality with respect to C cycle observations (Meyerholt & Zaehle, 2018). In order to better constrain N loss process parameterization, there is a need for ecosystem-scale experiments or monitoring studies that focus on quantification of the relevant N-loss pathways and their systematic use for model evaluation. An alternative source of information may be linked to the use of 15 N abundance data, which could be used both in a process-based manner by comparison to 15 N enrichment studies providing insight into the fate of N in terrestrial ecosystems (Goodale, 2017), or by attempting to interpret spatial and or temporal trends of 15 N abundance in terms of changes in the terrestrial N cycle (Craine et al., 2018(Craine et al., , 2019Hiltbrunner, Körner, Meier, Braun, & Kahmen, 2019;Houlton, Marklein, & Bai, 2015).

| Limitations of the study
Further alternative process representations in current N-enabled TBMs exist, for example the competition of plants and soil organisms for N or the downregulation of photosynthesis (Zaehle & Dalmonech, 2011). These additional processes would likely increase the effect of the internal (C:N) versus external N cycle (N in-and output) uncertainty, but are unlikely to affect our conclusions on the importance of BNF and N loss process representations, given the magnitude of the projected change particularly in BNF and the long timescale at which these changes occur. One process only implicitly considered in the ensemble of this study is nutrient mining in response to increased below-ground C allocation under elevated CO 2 .
In the ensemble, models with flexible C:N stoichiometry allow for a redistribution of N from soil to vegetation under N stress, because of increasing SOM C:N ratios. However, this process has been insufficient to explain the observed redistribution of soil N to the vegetation in free-air CO 2 enrichment experiments at Duke and ORNL forest . Increased below-ground C allocation has been found to increase SOM turnover (Drake et al., 2013;Hungate et al., 2009), and has been hypothesized to contribute to the sustained response of N-limited forests to elevated CO 2 in ecotomycorrhizal-dominated forests (Norby et al., 2017;Terrer, Vicca, Hungate, Phillips, & Prentice, 2016). One recent global model accounting for this effect suggested a substantial redistribution effect of N, essentially increasing land C storage, as increased soil C losses are more than compensated by increased above-ground C storage following this indirect fertilization effect (Sulman et al., 2019). This could potentially alter our conclusion about the importance of BNF in determining future N limitation in favour of processes determining ecosystem internal N cycling.
The absolute values of land C uptake and its sensitivities to environmental drivers depend on the C cycle assumptions built into O-CN. Nevertheless, we believe that the insights into sign and magnitude of the N cycle process contributions ( Figure 6) will be transferrable to other TBMs. Each of the N process algorithms applied in this study (Table 1) relies on few, insufficiently constrained parameters. We chose the parameter values in accordance with previous work (Meyerholt & Zaehle, 2015Meyerholt et al., 2016), which resulted in flux and stock magnitudes after model spin-up commensurate with current understanding (Supporting Information, Table S1). The simulated process responses to perturbations may be sensitive to the precise parameterization of the algorithms. Therefore our ensemble range may be a conservative estimate of the full N cycle uncertainty range. Considering that most algorithms respond linearly to small parameter perturbations (e.g. Meyerholt et al., 2016), it is unlikely that this uncertainty will affect the general conclusions drawn here about the relative importance and sign of the individual process contributions or the overall sign and range of the N effect.
As many other studies (Arora et al., 2019), this study does not consider the limitation of the terrestrial C cycle by the phosphorus (P) cycle. A first field-scale free-air CO 2 enrichment experiment in an Australian Eucalypt forest with low P availability has shown that P limitation can severely constrain the response (Ellsworth et al., 2017). Model simulations made for this and a planned experiment in the P-limited Amazon rain forests have demonstrated that the effect of P availability can be much stronger than N effects in P-limited ecosystems, but large uncertainties persist in these simulations (Fleischer et al., 2019;Medlyn et al., 2016). Considering P in addition to N dynamics in our ensemble would likely result in attenuated CO 2 responses in predominantly P-limited ecosystems, such as the Amazon rain forest or large parts of Australia.

| Model evaluation
We evaluated the model ensemble against a range of contemporary C-cycle observations, and demonstrated good, but varying performance of the ensemble members. In our ensemble, the difference in high-latitude productivity results directly from N cycle processes, as the C-cycle only reference considered all other climate constraints but the N cycle. The larger than observed high-latitude productivity of the C-cycle control may be O-CN specific, and could likely be corrected without necessarily including a fully prognostic N cycle. The variable performance of current C-cycle models and C-N cycle models (Le Quéré et al., 2018) suggest that currently, these benchmarks cannot be conclusively used to demonstrate that N cycle representations are necessary to match contemporary observations. However, they do allow us to state that all ensemble members in our study represent terrestrial C-N dynamics to a degree comparable to other state-of-the-art in biosphere models and that they are therefore valid candidates to assess potential future biosphere dynamics.
Despite the good present-day performance, substantial spread remains in the projected N effect on future land C sequestration.
These results support the concept of an ensemble approach to make a more robust assessment about the likely effect of N on projections of the global C cycle. The results further imply that advanced global C-cycle benchmarking systems (Collier et al., 2018;Righi et al., 2019) may insufficiently constrain long-term dynamics of the terrestrial biosphere, as they do not sufficiently constrain N availability and its effect on terrestrial C cycling. However, we note that we have not fully exploited the available observations. In our ensemble, some of the uncertainty resides in high-latitude ecosystems. Applying regionally explicit benchmarks, such as satellite-derived estimates of vegetation C (Liu, Dijk, McCabe, Evans, & Jeu, 2013) would likely help to constrain the model spread better. However, these data would not directly constrain simulated N dynamics, but other model features more directly affecting vegetation turnover.
Reducing uncertainty in the N effect requires the development of adequate and robust constraints on the terrestrial N cycle and its effect on vegetation growth and SOM decomposition, which necessarily requires a multifaceted approach. There is a need to improve the general understanding of ecosystem N dynamics and their covariation with C and water fluxes through a systematic observation and reporting of relevant N cycle characteristics in ecosystem monitoring to fully make use of the available C cycle observations to constrain coupled C-N cycle models (e.g. Vicca et al., 2018). A complementary source of information could be obtained by testing the simulated covariation of plant stoichiometry with climate and soils within plant functional types using existing plant trait databases (e.g. Kattge et al., 2011). The increasing availability of hyperspectral remote sensing will provide a new tool to better evaluate regional trends in foliar N content or related leaf properties (Croft et al., 2020;Knyazikhin et al., 2013;Ollinger et al., 2008;Townsend, Serbin, Kruger, & Gamon, 2013).
However, the evaluation of biosphere models with observations from monitoring networks will always and unavoidably be confounded by co-occurring effects of global change which obscure the effects of C-N coupling. It will therefore remain important to perform hypothesis-driven model evaluation against manipulation experiments to test model formulations with respect to those processes that determine the mid-to long-term response of ecosystems to perturbations , potentially involving new field experimentation dedicated to understanding the longterm response of ecosystems to changes in C or N availability (Norby et al., 2016). Even with improved constraints on current terrestrial N dynamics, the divergence of simulation results over the 21st century suggests that the long-term effect of N dynamics will remain to some extend unconstrained by current observations. This highlights the need for an ensemble approach to adequately reflect uncertainties in our process understanding in projections of the future C cycle.

DATA AVA I L A B I L I T Y S TAT E M E N T
Globally aggregated output of the ensemble is available at doi: https://doi.org/10.5281/zenodo.3620826. The O-CN model code is available from the authors upon request.

ACK N OWLED G EM ENTS
The authors acknowledge funding from the European Union's Horizon 2020 research and innovation programme (grant agreement no. 641816; CRESCENDO), as well as the European Research Council (ERC; grant agreement no. 647204; QUINCY).

AUTH O R CO NTR I B UTI O N
J.M. and S.Z. designed the study and performed the analysis; J.M.
implemented, tested and calibrated the models; S.Z. and K.S. performed the simulations. All the authors contributed to the writing.