Carbon-nitrogen interactions in idealized simulations with JSBACH (version 3.10)

. Recent advances in the representation of soil carbon decomposition and carbon-nitrogen interactions implemented previously into separate versions of the land surface scheme JSBACH are here combined in a single version which is set to be used in the upcoming 6 th phase of coupled model intercomparison project (CMIP6). Here we demonstrate that the new version of JSBACH is able to reproduce the spatial variability in the reactive nitrogen loss pathways as derived from a compilation of δ 15 N data (R=.76, RMSE=.2, Taylor score=.83). The inclusion of carbon-nitrogen 5 interactions leads to a moderate reduction (-10%) of the carbon-concentration feedback ( β L ) and has a negligible effect on the sensitivity of the land carbon cycle to warming ( γ L ) compared to the same version of the model without carbon-nitrogen interactions in idealized simulations (1% increase in atmospheric carbon dioxide per yr ). In line with evidence from elevated carbon dioxide manipulation experiments, pronounced nitrogen scarcity is alleviated by (1) the accumulation of nitrogen due to enhanced nitrogen inputs by biological nitrogen ﬁxation and reduced losses by leaching


Introduction
The version of the Max-Planck-Institute Earth System Model (MPI-ESM) used in the 5 th phase of the coupled model intercomparison project (CMIP5) experienced pronounced biases in simulated soil carbon (Todd-Brown et al., 2013), soil hydrology (Hagemann and Stacke, 2014), and the lack of carbon-soil nutrient interactions Wieder et al., 2015), hampering the reliability of the simulated response of land system to increasing carbon dioxide (CO 2 ), climate and land use and land cover changes. Recent model developments addressed these issues (Goll et al., 2012Hagemann and Stacke, 2014) in separate versions of the land surface scheme of the MPI-ESM, JSBACH, but have not been yet combined in a single model version. 5 The projected carbon balance in JSBACH was substantially affected by recent model developments: The implementation of carbon-, nitrogen-and phosphorus interactions reduced accumulated land carbon uptake by 25% between 1860-2100 under a business as usual scenario (Goll et al., 2012), while the implementation of a new decomposition model (YASSO) reduced the accumulated land carbon uptake by about 60% in the same period . The exchange of the former CENTURY type soil decomposition model (Parton et al., 1993) with the YASSO decomposition model (Tuomi et al., 2008(Tuomi et al., , 2009(Tuomi et al., , 2011 10 improved the present-day state of the carbon cycle compared to observations as well as the response of decomposition to soil warming, and substantially reduced the uncertainties of land use change emissions for a given land use change scenario . The strong impact on the carbon balance of each of both developments underline the importance of combining them in a single version. The capacity of land ecosystems to increase their nitrogen storage as well as to enhance recycling of nitrogen in organic 15 matter are major constraints on their ability to increase carbon storage under elevated CO 2 concentrations (Hungate et al., 2003;Thomas et al., 2015;Liang et al., 2016). The respective response patterns of nitrogen processes governing the balance and turnover of organic nitrogen are crucial (Niu et al., 2016) to asses the likelihodd of the occurence of (progressive) nitrogen limitation (Luo et al., 2004). Recent advances in the interpretation of soil δ 15 N global data sets provide a promising tool (Houlton et al., 2015;Zhu and Riley, 2015) by allowing a more detailed evaluation of the nitrogen loss pathways in land 20 carbon-nitrogen models than previously done (e.g. Parida (2011); Goll et al. (2012)).
Since future scenarios of CO 2 concentrations differ among CMIP phases, an idealized setup of an annual increase in CO 2 concentration by 1% is used to foster the analysis of the carbon cycle feedbacks among models, and to compare emerging properties of different model versions in various CMIPs (Eyring et al., 2016). We adopt this approach taking advantage of existing simulations of climatic changes in this idealized setup of the CMIP5 intercomparison (Arora et al., 2013) to drive the 25 land surface model JSBACH uncoupled from the atmosphere and ocean components of the earth system model. This article documents the modifications to the soil carbon decomposition  and nitrogen cycle (Parida, 2011;Goll et al., 2012) submodels and the combination of both developments in a recent version of JSBACH including an advanced soil hydrological scheme (Hagemann and Stacke, 2014), scheduled to be used in CMIP6. We further analyzed the state of the nitrogen cycle using soil δ 15 N data and quantified the carbon cycle feedbacks to increasing CO 2 concentrations 30 and climate change. The analysis aims at facilitating the interpretation of the models dynamics in the upcoming round of CMIP experiments (Eyring et al., 2016), and allows a straightforward comparison to the result from the previous round of CMIP (Taylor et al., 2012).

Model description
The implementation of the nitrogen cycle and the soil carbon and litter decomposition model YASSO is described in detail in Parida (2011); Goll et al. (2012) and Goll et al. (2015), respectively. In the following, a brief summary of the major concepts is given and afterwards the modifications to the original developments needed to combine them are documented in detail. The 5 notation applied here follows Goll et al. (2012Goll et al. ( , 2015 and a scheme of the cycling of carbon and nitrogen as well as their interactions are given in Figure 1. The decomposition model (YASSO) is based on a compilation of litter decomposition and soil carbon data and distinguish organic matter fractions according to litter size and solubility (Tuomi et al., 2008(Tuomi et al., , 2009(Tuomi et al., , 2011. In JSBACH we use two litter size classes, which correspond to litter from non-lignified and lignified plant material . Each of the two litter 10 classes is further refined into four solubility classes (acid-soluble (C A ), water-soluble (C W ), ethanol-soluble (C E ), nonsoluble (C N )) (Eq. 1). One additional pool (C H ) represents humic, slowly-decomposing substances.
The interactions between nitrogen availability and carbon fluxes, namely primary productivity and decomposition, are based on the concept of CO 2 -induced nutrient limitation (CNL) (Goll et al., 2012). In this framework, we distinguish between CNL and background nutrient limitation. The latter is assumed to be indirectly considered in the original parametrization of 15 carbon cycle processes as they are based on measurements in present day ecosystems and therefore reflect present day nutrient conditions. CNL is an additional nutrient limitation caused by the increase in atmospheric CO 2 and is computed dynamically according to nutrient supply and demand. In case microbial and vegetation nitrogen demand cannot be met by the supply, all carbon fluxes of which the donor compartment has a higher C:N ratio than the receiving pool (i. e. the fluxes of carbon from the solubility classes pools to the humus pool) are down-regulated. The concept of CNL allows to introduce carbon-nitrogen 20 interactions to YASSO, as the needed conditions are met, e.g. the parametrization of YASSO indirectly reflects present day nutrient effects on decomposition as it is based on leaf litter experiments.
Following Goll et al. (2012), CNL affects the decomposition of all pools except the slowly-decomposing nutrient-rich pool . The litter decomposition data on which YASSO is based is not suited to link the fate of nitrogen in litter to the respective solubility pools. Therefore, we assume one single nitrogen pool representing all nitrogen linked to the four carbon 25 solubility pools per litter class (Eq. 7). This can be refined in the future if appropriate data becomes available.
The nitrogen cycling is primarily driven by carbon fluxes using constant N-to-C ratios of organic pools (Eqs. 6-7), with the exceptions of the non-lignified litter pool (Eq. 8) (Parida, 2011). Further exceptions are the processes linking the terrestrial carbon cycle with the atmosphere (biological nitrogen fixation and denitrification) and the aquatic systems (leaching) which are computed either as substrate-limited (Eqs. [15][16] or, for the case of biological nitrogen fixation (Eq. 14), as driven by demand 30 due to the ample supply of N 2 in the atmosphere (Parida, 2011). The nitrogen cycle and its interactions with the carbon cycle are not modified. The only exception is that the turnover times of the nitrogen litter and soil organic matter pools are derived from the YASSO decomposition model (Eq. 10) instead of the former decomposition model.
All parameters and variables are given in Table 1&2 Figure 1. Schematic representation of carbon (top) and nitrogen (below) cycling in JSBACH. Vegetation is represented by four pools: "active" (leaves and non-lignified tissue) and "wood" (stem and branches), "reserve" (sugar and starches) and "mobile" (labile nitrogen) (Goll et al., 2012). Dead organic matter is represented by "non-lignified litter", "lignified litter" (lignified litter and fast-decomposing soil organic matter), and "humus" (slow-decomposing organic matter) (Raddatz et al., 2007). All organic matter pools have fixed C:N ratios, except the pools "reserve", "labile" and "non-lignified litter". While the first two pools have no corresponding pool, the C:N ratio of the latter pool varies according the balance between immobilization demand and supply. The carbon in the litter compartment is further refined into the acid-soluble (A), water-soluble (W), ethanol-soluble (E), and non-soluble (N) compounds  which have no C:N ratio assigned. Soil mineral nitrogen is represented by a single pool (soil mineral pool). The opposing triangle marks carbon fluxes which are downregulated in case the nitrogen demand exceeds the nitrogen supply.

Nitrogen effect on decomposition
Matrix C describes the soil carbon pools (A, W, E, N ) of the two litter classes (i) in JSBACH, excluding recalcitrant humic substances (C H ): The dynamics of the soil carbon pools are described as as a function of climatic conditions (F); and matrix I i is the carbon input of type i to the soil. The dynamics of the humus pool (C H ) are described as where p H is the relative mass flow parameter and k H the decomposition rate of the humus pool. A detail description of decomposition can be found in the supplementary information of Goll et al. (2015), here we only focus on the modification of the original implementation.
The climate dependence of the decomposition rate factor k j,i of the carbon pools was originally implemented by Goll et al. (2015) based on Tuomi et al. (2008) as: where T is air temperature and P is precipitation, β 1 , β 2 ,γ are parameters, and α j,i are decomposition rates at references conditions (T = 0 and P → ∞) of pool i of litter class j. α j,i is the product of reference decomposition rate r j of solubility class j and the litter diameter d i of litter class i. YASSO uses precipitation instead of the more direct driver of soil moisture due the lack of adequate soil moisture observation to relate the decomposition date the model is based on (Tuomi et al., 2008). 20 We introduced a scaling factor, namely the nitrogen limitation factor (f N limit ), to account for the down-regulation of decomposition when nitrogen is in short supply f N limit is calculated based on a supply and demand approach (Parida, 2011;Goll et al., 2012). In a first step, potential carbon fluxes are computed from which the gross mineralisation, immobilization and plant uptake of mineral nitrogen is diagnosed. In 25 a second step, all fluxes consuming nitrogen (donor compartment has a higher C:N ratio than the receiving pool as well as plant  uptake) are down-regulated in case nitrogen demand cannot be met by the nitrogen supply. Hereby, a common scalar (f N limit ) (see appendix) is used thereby no assumption about the relative competitive strengths of microbial and plant consumption has to be made. In case nitrogen demand is met by the supply, the fluxes computed in the first step are taken as actual ones without any modification.
2.1.2 Dynamics of nitrogen in litter and soil organic matter 5 Nitrogen in litter and soil organic matter is separated into three pools, namely slowly-decomposing organic matter (humus) C H , lignified litter and fast decomposing organic matter C lw , as well as non-lignified litter and fast decomposing organic matter C la (Goll et al., 2012). We assigned each of the three nitrogen pools to one or more corresponding YASSO pools (Table 3). A refinement of the representation of nitrogen in decomposing material following strictly the carbon classification is not straightforward as the carbon pools (A, W, E, N ) defined by their respective solubility characteristics do not correspond to 10 substance classes with distinguished stoichiometries.
In JSBACH, nitrogen in compartments with a fixed N-to-C ratio, namely nitrogen in lignified litter (N lw ) as well as nitrogen in slowly decomposing organic matter (N H ), are derived from the corresponding YASSO carbon pools (C j,i ) by:  Goll et al. (2015) β2 -1.4 × 10 −3 • C −2 Temperature dependence of decomposition Goll et al. (2015) γ1 -1.21 m −1 Precipitation dependence of decomposition Goll et al. (2015) αj,i where r H is the N-to-C ratio of former slow carbon pool (C H ) now applied to the humus pool of YASSO, and r lw of former lignified (woody) litter pool (C lw ) now applied to the sum of the solubility class pools for lignified litter of YASSO.
The dynamics of nitrogen in non-lignified litter & fast decomposing organic matter (N la ) were not modified from the original nitrogen-enabled version of JSBACH (Parida, 2011) and are given by: where the first term describes the nitrogen influx from active, non-lignified plant tissue (N a ), the second term describes the flux of nitrogen from herbivores excrements which is not directly available to biota, and the third term arises from the nitrogen released by biological mineralization of litter and fast-decomposing soil organic matter. We assume that active plant material (N a ) is consumed by herbivores at a constant rate ( ) and immediately excreted (Parida, 2011). We separate the excrement into 10 labile (β 3 ) and fast decomposing (1 − β 3 ) nitrogen compounds, the latter enters the non-lignified litter pool (N la ).
The decomposition rate d la of nitrogen in litter and fast decomposing soil organic matter equals the decomposition rate of the sum of the YASSO carbon pools C A,a + C W,a + C E,a + C N,a and is given by so that d la can be derived from 15 d la = C A,a (t + 1) + C W,a (t + 1) + C E,a (t + 1) + C N,a (t + 1) C A,a (t) + C W,a (t) + C E,a (t) + C N,a (t) − 1 As we calculate potential decomposition fluxes in a first step to derive nitrogen demand (see Goll et al. 2012) we know the state of pools for time t and t + 1.
The dynamics of the soil mineral nitrogen (N smin ) were not modified and are given -as originally formulated by Parida where F extr is the net of fluxes connecting the compartments considered in the model and outside (here: biological dinitrogen (N 2 ) fixation, leaching, N 2 and nitrous oxide (N 2 O) emissions), D micr and D veg are the nitrogen demands of vegetation and microbes, respectively. Due to the lower nitrogen content of litter compared to humus, the decomposition of lignified and non-lignified litter corresponds to a net immobilization of nitrogen, which is part of the D micr . The term (r w − r lw )F C w lw 25 represents nitrogen leaching from freshly shedded wood given by the decomposition and the stoichiometries assigned to wood (r w ) and lignified litter (r lw ). The decomposition rate of nitrogen in the humus pool (N H ), d H , equals the decomposition rate of the corresponding YASSO carbon pool C H , k H . This rate according to Eq.(4) is given by: Note that there is no nutrient effect on the decomposition of N H and k H is calculated exactly like described in Goll et al. (2015).
For the calculation of the microbial (soil) nutrient demand (D micr ) we substituted the pools C la and C H with the corresponding YASSO pools in Eq.(15) of Goll et al 2012: The fluxes F C la s and F C lw s are the net fluxes of carbon to the humus from the solubility pools (AWEN) of non-lignified and lignified litter, respectively. F C la a and F C lw a are the respective sums of respiration fluxes of the AWEN pools.

The processes governing the terrestrial nitrogen balance in JSBACH
Nitrogen enters terrestrial ecosystems by biological nitrogen fixation (BNF), as well as atmospheric deposition, while nitrogen is lost via leaching, erosion and denitrification. 10 BNF in global models is commonly represented by an empirical relationship based on a compilation of site measurements (Cleveland et al., 1999). Due to the lack of a process-based alternatives, we use this approach as described in Parida (2011) despite its shortcomings (Thomas et al., 2013;Sullivan et al., 2014;Wieder et al., 2015). In this approach BNF (BN F ) is derived from the annual average of daily net primary productivity (N P P ) using the empirical relationship between BNF and evapotranspiration (Thornton et al., 2007): where f emp = −0.003 dayg −1 (C) is an empirical relationship from Cleveland et al. (1999) , f bnf = 0.7 g(N)m −2 day −1 is a calibrated constant to achieve a global sum of BNF of 100 Mtyr −1 for a simulated NPP of 65 Gtyr −1 based on estimates for present day (Galloway et al., 2013;Ciais et al., 2013), and w N and w C the standard atomic weights of nitrogen and carbon, respectively.

20
The losses of nitrogen are given priority over immobilization and plant uptake each time step. Following Meixner and Bales (2002); Thornton et al. (2007); Parida (2011), daily losses by leaching are derived from dissolved nitrogen in soil water and the fraction of soil water lost to rivers per day (f h2o ) assuming a homogeneous distribution of mineral nitrogen (N smin ) in the soil volume : where f s is the fraction of mineral nitrogen (N smin ) in soil solution. f h2o is computed dynamically accounting for evapotranspiration, precipitation, and changes in the soil water storage using a 5 layer soil hydrological scheme (Hagemann and Stacke, 2014) Following Parida (2011); Goll et al. (2012), daily losses by denitrification are assumed to be at most 0.02% (k denit = 0.002 day −1 ) of the soil mineral (N smin ): where α is a JSBACH internal indicator of soil moisture stress [0-1] which is dynamically computed from soil moisture and used to scale biological activity (Raddatz et al., 2007). 5

Calibration & parametrization of the model
The parametrization of YASSO (version 3.20) and of the nitrogen cycle in JSBACH was not changed and is described in Goll et al. (2012Goll et al. ( , 2015. The only exception is the re-calibration of losses of nitrogen by leaching to the new hydrological model in JSBACH (Hagemann and Stacke, 2014). This is achieved, following Goll et al. (2012), by tuning the fraction of mineral nitrogen in soil solution (f s ) so that the assumption regarding the absence of CNL in the pre-industrial state is met which equals 10 to a negligible (<2%)) effect of nitrogen on global net primary productivity and carbon storage. We tuned the the fraction of soil mineral nitrogen in soil solution to f s = 0.1, which is comparable to fractions used in other models (Wang et al., 2010) as well as in observations (Hedin et al., 1995).

Simulation setup
We force the land surface model JSBACH with half hourly climatic data simulated by the MPI-ESM instead of running 15 JSBACH coupled with the atmospheric and ocean components of the MPI-ESM. Therefore, our simulations, in contrast to simulations of the MPI-ESM, do not account for the feedback between land and atmosphere with respect to the water and energy cycle. However, the resulting inconsistencies between climate and land surface should not change the results of the present study and are anyway partly implicit to the underlying CMIP5 simulations because of the prescribed atmospheric CO 2 levels in case of biogeochemical feedbacks (Taylor et al., 2012). For the sake of simplicity, we will refer to the JSBACH 20 simulations driven by the climate from respective ESM simulations, with the respective label of the ESM simulations.
The climatic forcing is derived from MPI-ESM simulations performed for the CMIP5 project (Table 4) (Taylor et al., 2012).

Spinup
The concept of CO 2 induced nutrient limitation (CNL) assumes that nitrogen effects on the carbon cycle are marginal under pre-industrial conditions. Therefore the cycles of carbon and nitrogen can be equilibrated in a two-step procedure in which the 25 carbon cycle is first brought into equilibrium (less than 1% change in global stocks per decade) using the climatic forcing from the pre-industrial control run (Goll et al., 2012). In a second step, we then initialize the nitrogen pools using the prescribed C:N ratios and the equilibrated carbon stocks as well as extremely high mineral nitrogen pools. The model is run again with the climatic forcing from the pre-industrial control run to equilibrate mineral nitrogen dynamics using the same criterion as for the first step. The resulting length of the simulation is 5.5 kyr and 2.6 kyr for step one and step two, respectively. Atmospheric Table 4. Simulations performed with JSBACH with and without carbon-nitrogen interactions using climatic forcing from MPI-ESM simulations performed for the CMIP5 project (Taylor et al., 2012 To analyze the effect of nitrogen limitation on the response of the land carbon cycle to increasing CO 2 concentration and climate change, we perform simulations with JSBACH with and without activated nitrogen cycle (Table 4). The simulations are forced with the climatic conditions from a set of 140 yr long CMIP5 simulations with the MPI-ESM in which atmospheric CO 2 concentration increases at a rate of 1% yr −1 from pre-industrial values until concentration quadruples (Arora et al., 2013). 5 The set of MPI-ESM simulations consist of a simulation where increasing CO 2 affects the climate but not the terrestrial biogeochemistry (radiatively coupled), a second simulation where increasing CO 2 affects the terrestrial biogeochemistry but not the climate (biogeochemically coupled), and a third simulation where increasing CO 2 affects both, climate and biogeochemistry (fully coupled). The biogeochemically-coupled and the radiatively coupled simulations allow us to disentangle the carbon-concentration feedback β L and carbon-climate feedback γ L , respectively (Friedlingstein et al., 2006;Arora et al., 2013). 10 β L is derived from the biogeochemically coupled simulations by dividing the difference in the total land carbon between the first and the last decade by the difference in the atmospheric CO 2 concentration of the same periods. γ L is derived from the radiatively coupled simulations by dividing the difference in the total land carbon between the first and the last decade by the difference in global land temperature of the same periods. The MPI-ESM simulations do not include the confounding effects of changes in land use, non-CO 2 greenhouse gases, aerosols, etc., and so provide a controlled experiment with which to compare for MAP -i.e. using indices more closely related to the governing processes; and (c) by using non-linear regression to fit a statistical model that is explicitly based on the isotopic mass balance equation of (Houlton and Bai, 2009).
The fraction of N loss in gaseous form (f denit ) was estimated using the principle described by e.g. Houlton and Bai (2009); Bai et al. (2012), but using a process-based statistical model for the relationship between soil δ15N data and environmental 20 predictors fitted to publicly available data on soil δ15N (Patino et al., 2009;Cheng et al., 2009;McCarthy and Pataki, 2010;Fang et al., 2011Fang et al., , 2013Hilton et al., 2013;Peri et al., 2012;Viani et al., 2011;Sommer et al., 2012;Yi and Yang, 2007). It was assumed that soil δ 15 N reflects the source (atmospheric) δ 15 N modified by isotopic discrimination that occurs during leaching (slight) and gaseous losses (much larger). For simplicity, the source δ 15 N was assumed to be zero and discrimination during leaching was neglected. Mean annual runoff (q, in [mm]) was estimated from precipitation and potential evapotranspiration 25 following (Zhang et al., 2004), with ω = 3. Following Xu-Ri et al. (2008) we assumed that leaching losses increase to a maximum dependent on soil water capacity, yielding an annual runoff factor f(q): with W max = 150 mm. Mean monthly soil temperatures (T m , in K) were estimated for 0.25 m depth following Campbell and Norman. We assigned a generic activation energy of E a = 55kJmol −1 K −1 (Canion et al., 2014) and summed the monthly index values over the 12 months yielding the annual soil temperature factor f (T ), where T ref = 293K.

5
The data were then fitted via the gaseous discrimination factor and a constant k by non-linear least-squares regression to the relationship where δ is soil δ 15 N, δ 0 is the δ 15 N of the N inputs. Assuming the leaching discrimination factor is 0, f denit can be expressed from the first principle (Houlton and Bai, 2009). Re-arranging Equation 20 and 21 we get A spatial map of f denit was derived from the empirical relationship between temperature, runoff and f denit using simulated values of f (q) and f m (T m ) from JSBACH. Thereby, model biases in climate are accounted for in the data derived f denit which 15 allows a straightforward comparison with simulated f denit . In addition, we derived maps of f denit based on monthly grids of observed mean climate from 1961-1990 covering the global land surface at a 10 minute spatial resolution (CRU CL2.0) (New et al., 2002) which are shown in the appendix.

Results & discussion
3.1 Model evaluation: pre-industrial state 20 The model simulates nitrogen stocks and fluxes under pre-industrial conditions well within the wide range of the few available observation based estimates ( Table 5). Most of the estimates are for present day conditions and thus are not directly comparable due to the human influence on the nitrogen cycle (Galloway et al., 2013).
The organic nitrogen stocks and fluxes are given by the prescribed C:N stoichiometry and the state variables of the carbon cycle and thus are not affected by the changes we introduced here, except for non-lignified litter and fast decomposing soil 25 organic matter which shows in general good agreement with observed C:N ratios for most biomes (see appendix). Mineralisation of organic nitrogen is the major source of nitrogen for vegetation and the simulated flux is less than existing model based  Zaehle et al., 2010). In models, nitrogen mineralisation is solely a by-product of decomposition of soil organic carbon and we thus attribute the differences between simulated mineralisation to the use of YASSO decomposition model compared to the use of the CENTRUY decomposition model Zaehle et al., 2010) as the soil C:N stoichiometries are comparable among models. We refer to the evaluation of the carbon cycle in JSBACH elsewhere (Anav et al., 2013;Goll et al., 2015), as the concept of CO 2 induced nutrient limitation 5 prevents an effect of nitrogen on the carbon cycle under pre-industrial CO 2 concentrations.
Estimates of global fluxes and stocks of nitrogen are often lacking or associated with large uncertainties, thus a detailed analysis of the simulated nitrogen cycle is hampered (Zaehle, 2013). However, recent advances in the use of δ 15 N data (Houlton et al., 2015), which are one of the few sources of spatially extensive data relevant to the nitrogen cycle, allow the evaluation of the respective importance of nitrogen loss pathways in space. Due to the different environmental controls of the loss pathways, 3.2 Changes in the land carbon cycle in the 1% CO 2 increase simulations 5 JSBACH simulates a strong increase in net plant productivity (NPP) due to increasing CO 2 from pre-industrial level to 4 × pre-industrial level (Figure 3). The simulated increase in NPP of 16.0% for a rise in atmospheric CO 2 from 370 ppm to 550 ppm is somewhat lower than the estimated increase of 25% from 4 free air carbon dioxide enrichment (FACE) experiments  (Table 6). A lower increase than in the FACE experiment can be expected as the long-term effect of elevated CO 2 is likely to be less than the one derived from short duration FACE experiments (Norby et al., 2010) on early successional forests . The simulated increase in GPP of 23.1% for an increase in atmospheric CO 2 concentration from the level of the year 1900 to 2013 is close to the 25% increase for the same increase in CO 2 estimated from intramolecular isotope distributions (isotopomers), a methodology for detecting shifts in plant carbon metabolism over long times (Ehlers

Changes in the land nitrogen cycle in the 1% CO 2 increase simulations
Increasing atmospheric CO 2 concentration leads to the accumulation of nitrogen in the terrestrial system (Figure 4a An increase in the nitrogen stock by 8.5% was found in a 15yr ecosystem scale CO 2 enrichment experiment (Shi et al., 2016), which is more pronounced than the simulated increase in the nitrogen stock of 3.2% for a comparable increase in CO 2 (year 29-69; 369-551ppm). A strong stimulation of BNF by 44%, with a strong decline in leaching by 42% and no significant changes in denitrification, mineralisation and soil organic nitrogen were found in a compilation of CO 2 enrichment experiments (Liang et al., 2016). However the representativeness of these findings was questioned recently (Rütting, 2016). In addition to that, CO 2 was increased abruptly in CO 2 enrichment experiments while it increased gradually in our simulations.
As the different nitrogen processes have different response patterns, t hey are likely to react differently to an abrupt than to a gradual increase in CO 2 . Although the relative contributions of reduced losses and increased inputs to an accumulation remain 10 somewhat elusive due to methodological biases (Rütting, 2016) and limited data, an accumulation of nitrogen under elevated CO 2 is a plausible scenario.
We find that the processes governing the nitrogen balance operate on different time scales ( Figure 5). The mineral nitrogen stocks decline from 1.33 Gt to 1.06 Gt during the first 35 yr and thereby reduce the substrate driven nitrogen losses ( Figure 5) from 98 Mtyr −1 to 83 Mtyr −1 . However, losses by leaching start to increase afterwards and are higher by the end of the 15 simulation than at the start. While the reduction in losses and gains in inputs contribute to equally parts to the accumulation in the first decades of the simulation, the stimulated BNF dominates the accumulation in later ( Figure 5). This highlights the importance of long term manipulation experiments for improving our understanding about the long term effects of increasing CO 2 on the terrestrial biosphere.
We find that the effect of increasing CO 2 and climate change on the nitrogen balance differ. Elevated CO 2 alone leads to a 20 shift from inorganic nitrogen to organic nitrogen (Figure 4a&b), whereas climate change is dampening this shift as warming stimulates the decomposition of organic nitrogen ( Figure 4d) and thereby slow down the progressive immobilization of mineral nitrogen into biomass and soil organic matter. Climate change alone leads to a loss of nitrogen from the system (Figure4a): The enhanced net mineralisation of organic nitrogen due to warming leads to increased losses of nitrogen via leaching and denitrification. We further found that changes in the water cycle due to climate change are increasing losses by leaching, while 25 denitrification follows primarily the changes in substrate availability despite the influence of soil moisture in Eq. 16 ( Figure 5).
Overall, we find that total nitrogen losses are intensified relative to the substrate availability at the end compared to the start of the simulations ( Figure 5).
The simulated increase in tightness of the nitrogen cycle as mineral nitrogen stocks deplete is in line with the substrate-based mechanisms proposed based on recent compilation of ecosystem nitrogen addition experiments (Niu et al., 2016), in which 30 the mineral nitrogen exerts a major control on the mineral nitrogen consuming processes. However, the respective observed response patterns of ecosystem nitrogen processes remain to a large degree unknown and are represented in a strongly simplified way in JSBACH . In general, we find that the effect of nitrogen availability on carbon storage is rather moderate in all simulation ( Figure 3) due to the adjustments of the nitrogen balance to changes in the carbon cycle ( Figure 5).

The effect of nitrogen on the carbon feedbacks
We quantify the strengths of the climate carbon-feedback (γ L ) and the carbon-concentration feedback (β L ) from the radiatively "coupled" and biogeochemically "coupled" simulations, respectively (Figure 6a,c). Both land feedbacks, with a global β L of 0.61 Pg(C)ppm −1 and a global γ L of -27.5 Pg(C) • C −1 , are 34% and 53% smaller than the multi-model averages of the CMIP5 models, despite that the CMIP5 version of JSBACH simulated global β L and γ L which are 59% and 42% times larger 5 than the average of CMIP5 models (Table 7).
In CMIP5 the two ESMs with nitrogen limitation (which shared the same terrestrial biosphere component) had feedback strengths 70-75% lower than averaged across models (Arora et al., 2013), suggesting a prominent role of nitrogen in dampening both carbon cycle feedbacks. Here, we find that the dampened response is primarily related to the modifications of the soil and litter carbon decomposition module, rather than to the inclusion of the nitrogen cycle. The global β L in the simulation without 10 nitrogen cycle is only 10% larger than in the simulation with nitrogen, while there is hardly any difference in global γ L between simulation, as the small positive and negative differences cancel out ( Figure 6). The contrasting findings regarding the effect of nitrogen on the land carbon feedbacks illustrates the need of a multitude of carbon-nitrogen models to draw general conclusions.
The large difference between the CMIP5 version of JSBACH and the recent version described here can primarily be attributed to the new decomposition model. The drastically reduced γ L is primarily caused by the smaller initial soil carbon stock, as well as by the long-term acclimation of decomposition to warming due to substrate depletion in YASSO . 5 The importance of the initial soil carbon stock for the carbon losses due to warming was illustrated by Todd-Brown et al.
(2014). The lower β L can be attributed to a much lower fraction of biomass which is converted into stable soil organic matter.
Therefore, the increasing productivity translates to a much smaller increase in carbon storage. JSBACH does not account for the stimulation of decomposition of recalcitrant carbon under elevated CO 2 due to increases in labile organic matter (from of priming effect) observed in lab incubation experiments (Kuzyakov et al., 2000), which could potentially alter the response of 10 soil carbon to increasing CO 2 , but for which there is no evidence of its relevance on multi-decadal time scales (Cardinael et al., 2015).
A dampening effect of the nitrogen cycle on the response of the terrestrial carbon cycle to climate change and increasing CO 2 is in line with the majority of carbon-nitrogen model studies (Ciais et al., 2013;Zaehle, 2013), but not all (Esser et al., 2011;Warlind et al., 2014). The mechanism by which nitrogen can dampen β L is via an increasing decoupling of gross pri-15 mary productivity from biomass accumulation under increasing CO 2 concentration: the incorporation of nitrogen into biomass reduces the mineral nitrogen availability (Luo et al., 2004;Liang et al., 2016) which negatively affects growth (Norby et al., 2010) and increases root respiration (Vicca et al., 2012;Mccormack et al., 2015). The dampening of γ L is mainly via an enhanced nitrogen mineralisation in cold regions due to warming (Zaehle, 2013;Warlind et al., 2014;Koven et al., 2015), which cannot be fully captured by JSBACH due to assumption of CO 2 -induced nutrient limitation, and therefore the model is prone 20 to underestimate the effect of nitrogen on γ L .

Model limitations and future development directions
The current understanding of processes governing the terrestrial nitrogen balance is still rather limited (Zaehle, 2013), and several processes which might be of importance, in particular stand dynamics  which can potentially alter biomass turnover (Brienen et al., 2015), interactions between plants and microbes which can stimulate nitrogen scarcity 25 by non-altruistic symbioses (Franklin et al., 2014), plasticity in stoichiometry and leaf nutrient recycling Meyerholt and Zaehle, 2015), and the availability of other nutrients (Goll et al., 2012) are not represented in JSBACH. In addition, the loss of organic matter in general due to erosion, although potentially of importance (Lal , 2002), is not yet represented in global land surface models, but development are underway (Naipal et al., 2015(Naipal et al., , 2016. Due to the concept of CO 2 -induced nutrient limitation in JSBACH the nitrogen cycle serves primarily as an additional 30 constraint on the carbon uptake. The advantage of the approach is its low complexity and avoidance of assumptions about the initial state of nutrient limitation thereby taking into account (1) the lack of data regarding the nitrogen cycle (Zaehle, 2013) as well as (2) the large uncertainty about the nutrient constraint on plant productivity (Letters et al., 2007;Zaehle, 2013).
The shortcomings of this approach are that it limits the applicability of the model to carbon cycle projections for scenarios of increasing atmospheric CO 2 and that it cannot capture any stimulation of the plant productivity due to changes in nitrogen availability itself: In addition to direct increase in nitrogen availability by nitrogen deposition and fertilization, a stimulation of plant productivity can occur due to reduced losses of nitrogen by pathways which are not under control of biota, like fire, leaching, or erosion . As a result, the model might underestimate the importance of nitrogen cycling for carbon uptake under elevated CO 2 .

5
Regarding the processes resolved in JSBACH, the extent of BNF increase has to be regarded as highly uncertain, despite its agreement with short-term experiments (Liang et al., 2016): the formulation of BNF used here is based on an empirical correlation between evapotranspiration and BNF and therefore the rate at which BNF rates increase strictly follows the increase in productivity whereas in reality the different processes leading to changes in BNF on ecosystem scales operate on a different time scale: The control of plants on their symbiotic partners via glucose export, and in case of nodules via oxygen regulation, 10 result in changes in BNF from hours to months. On longer time scales, the composition of the ecosystem, namely the fraction of BNF associated species, affects nitrogen inputs to the system. While for tropical ecosystems there is evidence that any governing mechanism(s) ought to operate at a synoptic scale (Hedin et al., 2009), higher latitude system might experience longer lag times. Additionally, other nutrients, like phosphorus or molybdenum, might slow down or reduce the potential of BNF to increase (Vitousek et al., 2013). BNF models which better resolve the governing mechanisms, for example (Gerber 15 et al., 2010;Fisher et al., 2012), should be incorporated into ESMs to increase the reliability of the simulated pace of changes in BNF (Meyerholt et al., 2016).
Models which simulate simultaneous competition for soil nitrogen substrates by multiple processes match the observed patterns of nitrogen losses better than models like JSBACH which are based on sequential competition (Niu et al., 2016). Here we find that despite the sequential competition, the simulated behavior is in general agreement with the dynamics of substrate-20 based mechanisms derived from manipulation experiments (Niu et al., 2016) and the spatial variability in the respective loss pathways is to a large degree in line with δ 15 N-derived patterns, despite the low performance of another sequential competition model (Houlton et al., 2015;Zhu and Riley, 2015) The high-latitude permafrost processes are not represented here but were shown to be of importance for the effect on warming on carbon and nitrogen losses. Permafrost regions store about 1,000 Gt C within the upper few meters of soil (Hugelius et al., 25 2014). The thawing of permafrost and deepening of the active layer in response to global warming can potentially lead to a much stronger climate carbon feedback (Schneider Von Deimling et al., 2012;Schuur et al., 2015). The recent study of Koven et al. (2015) with the CLM model showed these carbon losses in high latitudes can be partly offset by increased nitrogen mineralisation, and in turn productivity and input to the soils. manipulation experiments is a research priority to increase the model reliability (Luo et al., 2012). In this study, the use of the ESM climatology is justified by a focus on the feedback analysis in the framework of idealized simulations as suggested in the 35 climate-carbon cycle model intercomparison project C4MIP (Anav et al., 2013;Jones et al., 2016). For further evaluation of the nitrogen limitation, however the preferable setup includes site-level simulations driven by observation-derived climatology .

Conclusions
The simulated response of primary productivity to increasing CO 2 , simulated litter stoichiometries, as well as the simulated 5 spatial variability in nitrogen loss pathways are in good agreement with observation based estimates. Here we show that a simple representation of mineral nitrogen dynamics can achieve a high agreement with observation in respect to nitrogen loss pathways. Further refinements of denitrification should address the relationship between denitrification and low soil moisture availability and as well as introduce a temperature scaling function.
The effect of nitrogen cycling on the land carbon uptake in idealized simulations with JSBACH is globally minor, but 10 not negligible. In particular, the carbon-concentration feedback is affected by mineral nitrogen availability, but the extent is moderate compared to earlier studies (Arora et al., 2013;Zaehle, 2013). During the first decades of the simulations, nitrogen limitation is circumpassed by a strong initial decline in loss terms in combination with increases in biological nitrogen fixation.
Afterwards progressive increases in biological nitrogen fixation drive the accumulation of nitrogen in ecosystems. On top of that, warming enhances mineralisation and counteracts the immobilization of nitrogen in biomass. Our study is in line with the 15 majority of carbon dioxide enrichment studies (Liang et al., 2016), showing that progressive nitrogen limitation under elevated carbon dioxide concentrations is less likely to occur than originally suggested (Luo et al., 2004).
The timescale and the extent to which the nitrogen cycle adjusts to increasing carbon dioxide and changing climate depend on the response of the input and loss processes of nitrogen as well as turnover of organic nitrogen. Here, we illustrate that the processes counteracting nitrogen scarcity operate on different time scales and have different trajectories due to differences 20 in the respective environmental drivers, which indicates a picture more complicated than drawn from environmental response functions (Niu et al., 2016;Liang et al., 2016). It is difficult to assess to what extent the timescales in our experiments are realistic, as timescale on which these processes operate is well beyond the typical duration of manipulation experiments.
Our results suggest that nitrogen limitation of land carbon uptake of natural ecosystems could be temporally restricted, being the result of the inertia of the balancing processes (Altabet et al., 1995;Hedin et al., 2009). Ultimately, other nutrients 25 like phosphorus which sources are depleted over time, are likely to dominate the long-term capacity of carbon storage.

Code availability
The JSBACH model version 3.10 used here includes the soil module YASSO and nitrogen components. The model version corresponds to the revision 8691 from the 19th July 2016 in the Apache version control system (SVN) of the Max Planck Institute for Meteorology (https://svn.zmaw.de/svn/cosmos/branches/mpiesm-landveg). This version will be used in 30 the CMIP6 simulations, where other components (landuse, dynamic vegetation, fire) will be included as well. The source   The nitrogen limitation factor, f N limit , is calculated based on a supply and demand approach (Parida, 2011;Goll et al., 2012). In a first step, potential carbon fluxes are computed from which the gross mineralisation, immobilization (D micr ) and plant uptake of mineral nitrogen (D veg ) is diagnosed. In a second step, all fluxes consuming nitrogen (donor compartment has a higher C:N ratio than the receiving pool as well as plant uptake) are down-regulated in case nitrogen demand cannot be met by 5 the nitrogen supply. Hereby, a common scalar (f N limit ) is used thereby no assumption about the relative competitive strengths of microbial and plant consumption has to be made.
where the term in square bracket is the maximum rate at which the soil mineral nitrogen pool can supply nitrogen. Note that in the discretized formulation the mineral nitrogen pool can at most be depleted during a single model time step (∆t). We thus set this maximum rate to dNsmin ∆t .

Appendix B: Evaluation of dynamically computed C:N ratios
The only ecosystem compartment in JSBACH which has a flexible stoichiometry is non-lignified litter and fast-decomposing 5 organic matter. The simulated carbon to nitrogen ratios of this compartment for the six plant functional types in JSBACH are in rather good agreement with observations of foliage litter from the ART-DECO database (Table A1)  The data for tropical species is very scarce and the variability among species is large, which hamper the interpretation of the mismatch between model and observation for the tropical broadleaved deciduous trees.

Appendix C: Consistency of nitrogen loss pathways with earlier estimates
The reconstructed f denit map from observed climatology ( Figure (A1) is generally similar to one presented by Houlton et al. (2015), with high fractions (ca 80%) in the tropics and mid-latitude deserts, a strong gradient of decreasing fractions with 15 decreasing temperature towards high altitudes and latitudes, and values in the range 0-20% reached in cold, wet climates in the north. However, some differences are apparent, most obviously connected with the use of mean annual temperature (MAT) by Houlton et al. (2015) to index microbial activity. MAT becomes extremely low in Eurasia towards the northeast, for example, and accordingly, Houlton et al. (2015) estimates of the denitrification fraction become very low there. Craine et al. (2015) noted that climates with very low MAT (including sites in NE Siberia) showed anomalous values of soil δ 15 N, more similar to 20 those of warmer climates. Our approach takes account of this by the use of an index that is much more responsive to the warm summers than to the extreme cold winters found in hypercontinental climates.
When simulated climatology is used to upscale the empirical relationship between temperature, runoff and soil δ 15 N, the influence of biases in simulated climatology on f denit become apparent. The overestimation of precipitation and subsequently runoff of about 20% in MPI-ESM (Weedon et al., 2011;Hagemann et al., 2013) leads to a pronounced peak in the histogram 25 of f denit at about 0.1-0.2 ( Figure A1), which is mostly in the mid and high latitudes regions in northern hemisphere.  Cardinael, R., Eglin, T., Guenet, B., Neill, C., Houot, S., and Chenu, C.: Is priming effect a significant process for long-term SOC dynam-