Controls of terrestrial ecosystem nitrogen loss on simulated productivity responses to elevated CO 2

The availability of nitrogen is one of the primary controls on plant growth. Terrestrial ecosystem nitrogen availability is not only determined by inputs from fixation or deposition, but also regulated by the rates with which nitrogen is lost through various pathways. Estimates of large-scale nitrogen loss rates have been 10 associated with considerable uncertainty, as process rates and controlling factors of the different loss pathways have been difficult to characterize in the field. Therefore, the nitrogen loss representations in terrestrial biosphere models vary substantially, adding to nitrogen cycle-related uncertainty and resulting in varying predictions of how the biospheric carbon sink will evolve under future scenarios of elevated atmospheric CO2. Here, we test three commonly applied approaches to represent ecosystem level nitrogen loss in a common carbon-nitrogen 15 terrestrial biosphere model with respect to their impact on projections of the effect of elevated CO2. We find that despite differences in predicted responses of nitrogen loss rates to elevated CO2 and climate forcing, the variety of nitrogen loss representation between models only leads to small variety in carbon sink predictions. The nitrogen loss responses are particularly uncertain in the boreal and tropical regions, where plant growth is strongly nitrogen limited or nitrogen turnover rates are usually high, respectively. This highlights the need for 20 better representation of nitrogen loss fluxes through global measurements to inform models.


Introduction
Given potential detrimental implications of increasing global atmospheric carbon (C) dioxide (CO 2 ) 25 concentrations, research on the terrestrial compartment of the Earth system has focussed on the potential for the biosphere to sequester atmospheric CO 2 in the future (Bonan, 2008). One important constraint to the terrestrial vegetation C sequestration potential may be the availability of plant nutrients, primarily nitrogen (N) (Hungate et al., 2003;Thornton et al., 2007;Gruber and Galloway, 2008;Zaehle et al., 2010a;Fernández-Martinez et al., 2014). 30 N enters natural terrestrial ecosystems from the atmosphere through deposition of reactive N species (formed by lightning or fossil-fuel combustion), as well as through biological N fixation (BNF). N is taken up from the soil by plants, eventually returned as organic litter and is incorporated into the soil organic matter or becomes mineralized, meaning that microbial activity converts organic N back to plant nutrients, namely ammonium (NH 4 + ), which can then be converted to nitrate (NO 3 -) during nitrification and be taken up by plants again. This 35 loop of plant N uptake from the soil and mineralization of organic N can be regarded as the internal N cycle.
However, soil N may also be lost from the ecosystem through gaseous and leaching loss processes (Fig.1).
Biogeosciences Discuss., https://doi.org /10.5194/bg-2018-232 Manuscript under review for journal Biogeosciences Discussion started: 4 June 2018 c Author(s) 2018. CC BY 4.0 License. Thereby, the loop of N entering ecosystems from the atmosphere and leaving them through loss processes and eventually ending up back in the atmosphere can be regarded as the external N cycle.
Numerous studies have demonstrated the effects of including N cycle dynamics when using terrestrial biosphere models (TBMs) to examine terrestrial C cycle responses to elevated atmospheric CO 2 concentrations (eCO 2 ) and climate change (Thornton et al., 2007;Sokolov et al., 2008;Zaehle et al., 2010a;Goll et al., 2012;Wania et al., 5 2012;Smith et al., 2014). The theoretical effect of previously defined, conceptual C-cycle feedbacks (Cox et al., 2000) may be altered by simulating C-N dependencies: Additional plant C assimilation under increased atmospheric CO 2 concentrations may be limited by N availability (Zaehle et al., 2010a). Increased global temperatures may not only increase ecosystem respiration and load the atmosphere with more CO 2 , but also stimulate N mineralization and provide more N to support plant C assimilation (Sokolov et al., 2008;Thornton et 10 al., 2009).
The magnitude of N effects on model predictions varies between studies, in part due to differences in how N cycle processes are formulated and included in the models (Zaehle and Dalmonech, 2011). Consequently, TBMs vary in their ability to reproduce the results of eCO 2 field experiments (Zaehle et al., 2014;Walker et al., 2015).
To gain understanding of the mechanics underlying this uncertainty, new studies have emerged that assess the 15 influence of variety in the representation of individual N cycle processes in model perturbation experiments Wieder et al., 2015;Meyerholt et al., 2016). However, a similar comparative study of N loss representation in TBMs is still lacking, although differences in N loss representation have been suspected in the past as a driving factor of model divergence in response to perturbation when different TBMs were compared (Thomas et al., 2013;Zaehle et al., 2014;Walker et al., 2015). 20 In nature, the pathways for gaseous N loss from ecosystems are manifold (Firestone and Davidson, 1989). The aerobic process of nitrification is the oxidation of NH 4 + to nitrite (NO 2 -) and then to NO 3 -. It is associated with the emission of nitric oxide (NO) and nitrous oxide (N 2 O) during the reduction of NO 2 when oxygen is limiting . Denitrification is carried out by denitrifying soil bacteria that are able to sequentially reduce oxidized forms of N (NO 3 -→ NO 2 -→ NO → N 2 O → dinitrogen (N 2 )) in the presence of organic matter to 25 produce molecular N 2 that is emitted to the atmosphere. Anaerobic conditions are a prerequisite, as oxygen depletion causes oxidized N to act as a substitute electron acceptor to denitrifying bacteria. In the process of denitrification, NO and N 2 O may also be emitted, which makes this mechanism particularly climate relevant.
Denitrification is considered the most important terrestrial N loss flux with 110 Tg N yr -1 estimated for the year 2000 (Bouwman et al., 2013). Volatilization of ammonia (NH 3 ) is of special importance in agriculture and may 30 take place when N from manure or fertilizers cannot react to form NH 4 + in the soil due to alkaline conditions or high soil temperatures and is lost to the atmosphere in its gaseous form.
As for the mechanisms of leaching N loss, the pathway of NH 4 + loss is adsorption to soil particles and the possible leaching of NH 4 + bound in this manner, especially in clay soils (Kowalenko and Cameron, 1976;Matschonat and Matzner, 1996). In contrast to NH 4 + , NO 3 in soils is prone to leaching losses due to its negative 35 charge and inability to adhere to soil particles. In scenarios of high precipitation or irrigation and high NO 3 concentrations, NO 3 can be lost to groundwater through vertical transport and thus cause pollution and reduced ecosystem productivity. This hydrological export of N can also affect dissolved organic N, a flux that has been shown to be of sizable magnitude at some sites (Perakis and Hedin, 2002;Gerber et al., 2010).
Despite such general understanding of the pertinent processes, the reason for the variety of N loss representations 40 in TBMs is the difficulty in properly characterizing N loss fluxes at large spatio-temporal scales in nature, given Biogeosciences Discuss., https://doi.org /10.5194/bg-2018-232 Manuscript under review for journal Biogeosciences Discussion started: 4 June 2018 c Author(s) 2018. CC BY 4.0 License. the strong variability in space and time of the associated trace gas and water fluxes (Boyer et al., 2006). In addition, the relevant fluxes are also very difficult to measure in the field, especially in the abundance needed to constrain global models. Therefore, modellers need to resort to educated guesses on how to represent this poorly constrained ecosystem flux. Model implementations vary between the application of generic loss terms (e.g. Thornton and Rosenbloom, 2005;Wang et al., 2010) and the explicit formulation of the constituting loss fluxes 5 by simulating the environmental conditions that are assumed to influence specific loss processes (nitrification, denitrification, leaching, fire;e.g. Xu-Ri and Prentice, 2008;Huang and Gerber, 2015). In the latter case, explicit treatment of gaseous N loss may even enable detailed estimates of the ecosystem emissions of the greenhouse gas N 2 O . Between these cases of simplified and complex formulations lie a number of TBMs that represent N loss fluxes at "intermediate" complexity (Yang et al., 2009;Goll et al., 2012;Smith et al., 10 2014). Models also differ with respect to the N cycle component from which the respective loss flux is derived.
Some models focus their simulation of gaseous N loss processes on soil N turnover, i.e. N mineralization (Thornton et al., 2007;Wang et al., 2010), while others base their calculations on the size of the soil inorganic N pool Zaehle and Friend, 2010). Some TBMs include leaching of dissolved organic N (DON) directly from the soil organic matter (SOM) N pool (Gerber et al., 2010;Smith et al., 2014). Such 15 heterogeneous representation of ecosystem level N losses is a particular limitation when attempting to estimate the effect of N limitation on terrestrial C sequestration, both at present and under future scenarios.
The aim of this study was to determine the extent to which variation between different N loss algorithms would influence simulated C sequestration responses to eCO 2 . We expected that the different paradigms of concentration-based and turnover-based N losses would lead to different predicted N loss magnitudes especially 20 under N stress. With simulated depletion of the inorganic soil N pool under eCO 2 , modelled N loss should decrease more strongly if it is concentration-based than when it is turnover-based. Assuming largely N-limited vegetation growth and fixed ecosystem N inputs, such differences in N loss rates could lead to different C sequestration responses, demonstrating that the choice of N loss formulation plays a notable role in shaping the predictions of C-N TBMs in simulated perturbation scenarios. 25 To examine the impact of different N loss algorithms in a TBM, we added two new alternative N loss modules to the O-CN TBM (Zaehle and Friend, 2010). The original O-CN N loss formulation was in part adopted from Xu-Ri and Prentice (2008) and largely bases gaseous N losses on the concentration of reactive N in the soil compartment. As alternatives, we added two more N loss algorithms that base the largest gaseous N loss flux on the N mineralization flux from soil organic matter to the soil pool of reactive N, inspired by the formulations 30 presented by Thornton and Rosenbloom (2005) and Wang et al. (2010). Thereby, our selection of N loss formulations encompassed the cases of complex and simplified algorithms mentioned above.
As a simple base scenario, we performed eCO 2 experiments with these three O-CN versions that only differed in their N loss algorithms at three temperate test sites that only differed in their climate forcing. This was done to examine how model functioning was affected in quasi-equilibirum and under eCO 2 regarding the calculated N 35 loss fluxes and the effect on C sequestration, and to illustrate the approximate climate sensitivity of these patterns. Next, we performed long-term simulations at a temperate site to gain insight into the centennial-scale effect of the three loss algorithms on the evolution of ecosystem N limitation under eCO 2 . We then performed eCO 2 simulations on a global 2°×2° grid using the three model versions to examine the effects of N loss variety in different ecosystem types that exhibited inherently different degrees of N limitation and climate regimes. 40 Biogeosciences Discuss., https://doi.org /10.5194/bg-2018-232 Manuscript under review for journal Biogeosciences Discussion started: 4 June 2018 c Author(s) 2018. CC BY 4.0 License.

Methods
In Sect. 2.1, we describe the three modular N loss algorithms. In Sect. 2.2, we describe the different eCO 2 experiments we performed.

O-CN terrestrial biosphere model and nitrogen loss formulations: 5
As a TBM framework, we used the O-CN model that was fully described in Zaehle and Friend (2010) and . We used the standard O-CN N loss formulation ("NL1") and added two more formulations as alternatives ("NL2", "NL3"), based on formulations used in other TBMs (Fig. 1). Here, we describe the three N loss formulations in detail, supplemented by Appendix B in case of the more complex "NL1" formulation. 10

NL1
The original O-CN representation of ecosystem N loss  follows the representation in the Lund-Potsdam-Jena-Dynamic-Nitrogen (LPJ-DyN) TBM (Xu-Ri and Prentice, 2008) with additions from the 15 DNDC denitrification and decomposition model . O-CN treats the nitrification and denitrification processes explicitly to determine gaseous losses of NO, NO 2 , N 2 O, and N 2 . In addition, NH 3 is subject to volatilization based on soil pH. Leaching of solute NH 4 + and NO 3 occurs in proportion to the soil water lost by soil drainage. A full description can be found in Appendix B.
The NL1 algorithm determines total ecosystem N loss (N L ; Eq. 1) as the sum of gaseous N losses from 20 nitrification (N nit ), denitrification (N denit ), and volatilization (N vol ), as well as leaching (N lea ): (1)

NL2 25
The NL2 approach, inspired by Wang et al. (2010), includes N loss fluxes based on soil N turnover and soil inorganic N concentration. A fixed fraction of the instantaneous net N mineralization flux is lost to the atmosphere to represent gaseous N losses associated with microbial soil processes of nitrification and denitrification (N L, gas ; Eq. 2). Gaseous N loss only occurs when net N mineralization is positive: 30 where f gl is the fraction (0.05) of the net N mineralization flux (N nm ) that is lost in gaseous form. To account for leaching losses (N L, lea ; Eq. 3), the total soil inorganic N pool (N min ) is reduced at every time step: where f ll is the fraction (0.5) of the soil inorganic N pool lost to leaching. The total ecosystem N loss per (N L ; Eq. 35 4) is then given by the sum of gaseous and leaching losses: Biogeosciences Discuss., https://doi.org /10.5194/bg-2018-232 Manuscript under review for journal Biogeosciences Discussion started: 4 June 2018 c Author(s) 2018. CC BY 4.0 License.

NL3
The N loss formulation NL3, inspired by Thornton and Rosenbloom (2005), describes sequential processes of gaseous loss during net N mineralization, gaseous loss from the soil inorganic N pool, and lastly leaching loss from the remaining soil inorganic N reservoir. Similar to NL2, the net N mineralization flux (N mn ) is 5 accompanied by a fractional denitrification flux (N L,g1 ; Eq. 5): where f g1 is the fraction (0.01) of the net N mineralization flux that is lost in gaseous form. Next, excess 10 inorganic N remaining in the soil after immobilization and plant N uptake (N min ) is subject to further gaseous loss representing volatilization and denitrification (N L,g2 ; Eq. 6):

15
where f g2 is the fraction (0.002) of the soil inorganic N pool lost in gaseous form. Any remaining inorganic N in the soil is then subject to fractional leaching loss in constant proportion (N L,l ; Eq. 7): 20 where f l is the fraction (0.1) of soil inorganic N lost to leaching. The total ecosystem N loss per time step (N L ; Eq. 8) is then given by the sum of gaseous and leaching losses: 25

Forcing and simulation protocol
We conducted three seperate sets of eCO 2 simulation experiments, two at the local and one at the global scale.

Local simulations I
The first set of local simulations was carried out at a representative temperate forest site ("S1") at the coordinates 6° E, 48° N. We included two more sites ("S2", "S3") that were identical to S1, with the exception that we increased air temperatures by 5 K ("S2") or doubled precipitation ("S3") relative to the climate forcing at S1, 35 thereby creating an ensemble of three "pseudosites". This was done to include the effect of climate variation, but without further confounding the results with influences from e.g. different soil and vegetation histories, keeping the effect of the N loss formulation as isolated as possible. For each model version and each pseudosite, the vegetation and soil C and N pools were spun up to equilibrium over 900 simulation years until the year 1700,  Quéré et al. (2016). Model spin-up used atmospheric CO 2 concentrations from the year 1860 (Le Quéré et al., 2016), 1850 N deposition rates according to Lamarque et al. (2010), BNF according to the "FOR" approach described by Meyerholt et al. (2016), and vegetation cover from the SYNMAP product (Jung et al., 2006). To limit the driving factors in this theoretical study, we disregarded crop vegetation and subsequently N fertilizer application, as well as land-use change. 5 The models were then run on a half-hourly time step for 313 simulation years using the climate forcing described

Local simulations II
The second set of local simulations was carried out only at the temperate "S1" site from Sect. 2.2.1, with the following modifications to the first set. Instead of 63 years as above, a different eCO 2 experiment was conducted over 300 simulation years. After spin-up, atmospheric CO 2 concentrations were kept constant at the 1860 level 15 (286 ppm) between 1700 and 1860. Between 1860 and 2006, atmospheric CO 2 increased according to ambient concentrations. Next, the 2006 concentration (380 ppm) was kept constant for the following 300 years to create the experiment control runs. For the treatment runs, atmospheric CO 2 was set to 580 ppm between the years 2006 and 2306.
Instead of the recorded climate data, randomly selected climate years for the "S1" site were used from the years 20 1901-1930 throughout all simulations. After the year 2006, atmospheric N deposition rates were kept constant to the 2006 level.

25
Global simulations were carried out on a global grid of 2° x 2° resolution. The ecosystem C and N pools were spun up to equilibrium for 1291 years until 1850, using 1700 vegetation cover (Hurtt et al., 2006), 2000 fertilizer application (Zaehle et al., 2010b), 1850 N deposition rates (Lamarque et al., 2010), 1850 ambient CO 2 concentrations, and 1901-1930 global climate data. The climate data was taken from the 5th Phase of the Since N input rates of deposition and fixation were fixed, the total N loss rates predicted by the models for the quasi-equilibrium simulation period (1700-1750) barely differed between pseudosites (Fig. 2a,e,i). At S1, NL3 calculated the lowest fraction of N loss through leaching because leaching was assumed to take place after plant N uptake and two gaseous loss pathways were accounted for, leaving less N in the inorganic soil pool for 10 leaching (Fig. 2a). Despite the differences in concept and detail of N loss representation, the proportions between gaseous loss and leaching for NL1 and NL2 were predicted to be more even. The simulations that employed the NL1 loss algorithm, however, generated more year-to-year variability than the other simulations, which generally applied across pseudosites. This was because losses in NL1 were mostly dependent on the N concentration in the inorganic soil pool, which underwent pronounced fluctuations from N plant uptake, losses 15 and inputs, and generally did not accumulate N over substantial periods of time. In contrast, NL2 and NL3 based a large portion of the calculated N losses on the mineralization flux ( Fig. 1), which was derived from the pool of soil organic matter (SOM). This pool was far larger in absolute amounts of N than the inorganic pool, making N loss a rather consistent flux in comparison as long as vegetation, i.e. substrate for N mineralization, was present.
Subjected to ten years of eCO 2 , the simulated ecosystems showed considerable variation in how N loss rates 20 were predicted to evolve, depending on the applied loss algorithm (Fig. 2b). While NL1 and NL3 predicted decreases in total loss rates, NL2 predicted increased gaseous N loss and a slight reduction of leaching rates, resulting in an overall increase of total N loss. The NL1 response was dominated by a decrease of gaseous loss.
In the O-CN TBM, eCO 2 led to increased plant growth, accompanied by increased plant N demand. When this demand was met through increased plant N uptake, the soil inorganic N pool became more strongly depleted 25 than it would have under ambient CO 2 concentrations. Gaseous N losses in NL1 decreased because they were mostly based on the soil inorganic N concentration.
In contrast, applying the NL2 model under eCO 2 at S1 resulted in an increase of gaseous N losses. This was a direct reflection of the exclusive dependence of gaseous loss on net N mineralization in NL2 (Fig. 1). In the O-CN TBM, the depletion of the soil inorganic N pool under eCO 2 led to an increased C:N ratio of SOM, which in 30 turn led to increased N release from the mineralization of organic material to steer SOM C:N back towards a target ratio (Zaehle and Friend, 2010). Thus, when NL2 was applied, the effect on gaseous N loss was an increase. Depletion of soil N also caused a decrease in (concentration-dependent) N leaching, however the net effect (total N loss) was positive (Fig. 2b).
The NL3 algorithm produced reduced total N loss under eCO 2 through reduced gaseous loss, albeit at a smaller 35 magnitude than NL1. Although NL3 featured a similar mineralization-dependence of gaseous N loss as NL2, we found that most gaseous loss change in NL3 occured independently of net N mineralization change under eCO 2 .
This meant that any gaseous loss increase that occured with increased net N mineralization was superseded here by the secondary, soil N concentration-based loss pathway that reduced gaseous N loss with soil N depletion ( Fig. 1). As with the NL1 formulation, the eCO 2 effect on leaching was negligible for NL3 at S1. 40 Biogeosciences Discuss., https://doi.org /10.5194/bg-2018-232 Manuscript under review for journal Biogeosciences Discussion started: 4 June 2018 c Author(s) 2018. CC BY 4.0 License. Despite these major model differences in predicted N loss changes under eCO 2 at S1, the model differences in predicted C sequestration changes (NPP; Fig. 2c) and ecosystem C accumulation during the experiment (Fig. 2d) were small. All model versions predicted NPP increases around 25% and growth of the total ecosystem C pool between 4 % and 5 %. Interestingly, the NL2 loss model predicted the largest increase in ecosystem C despite also predicting the only increase of ecosystem N loss. In particular, this largest C increase was simulated in the 5 SOM C pool, where subsequently the SOM C:N also increased stronger for NL2 than for the other loss algorithms.

10
For the quasi-equilibrium loss partitioning at the S2 pseudosite (Fig. 2e), the most prominent difference compared to the S1 site was that the +5 K change in air temperature caused a smaller leaching portion predicted by the NL1 algorithm. At higher temperatures, soil evaporation increased and soil drainage decreased, which for NL1 led to a reduction of leaching loss as a consequence of the coupling of drainage and leaching in this algorithm (see Appendix B). As this coupling was not applied in the NL2 and NL3 formulations, they predicted 15 barely any effect of higher air temperatures on the partitioning between gaseous loss and leaching.
Notable effects of higher temperatures on the eCO 2 responses of N loss rates (Fig. 2f) included a reduction in leaching for NL1 and a stronger gaseous loss reduction for NL3. With the reduced baseline leaching for NL1, the reduction of soil inorganic N concentrations under eCO 2 resulted in now notable reduction of leaching loss. For NL3, the gaseous loss reduction was dominated by the secondary, N concentration-based pathway. Given the 20 high temperatures at the S2 site, N mineralization rates were already at a high level and increased under eCO 2 at a lower rate than at the S1 site. The rather small fraction (0.01, see Sect. 2.1.3.) of the N mineralization flux that was lost in gaseous form in NL3 did not make for substantial loss through the primary, flux-based pathway.
Instead, gaseous losses were reduced stronger than at S1 because more inorganic N was left in the soil to be lost through the secondary loss pathway. 25 Compared to the S1 site, the temperature increase at S2 resulted in higher NPP responses (over 30%; Fig. 2g) and more ecosystem C accumulation under eCO 2 (over 5%; Fig. 2h). Higher temperatures led to higher gross primary productivity (GPP) responses to eCO 2 due to the sensitivity of modelled photosynthesis, which subsequently propagated to the NPP and C accumulation responses.

S3 site
When precipitation was doubled, the leaching portion of N loss for the NL1 formulation in quasi-equilibrium increased dramatically (Fig. 2i), owing to the dependence of leaching on drainage. Since in this state the total ecosystem N loss was essentially prescribed by the fixed rates of ecosystem N input, the gaseous loss portion 35 was minimized. This was further aided by a decreased nitrification rate in NL1 at S3 due to a decreased aerobic soil fraction (see Appendix B), which reduced the associated gaseous losses. While there was no effect of precipitation increase for NL2, the proportions of loss pathways changed for NL3 towards more leaching and less gaseous loss.
The precipitation increase brought about a number of changes to the sensitivity of N loss under eCO 2 (Fig. 2j). 40 For NL1, most of the N loss reduction was now simulated as reduced leaching, a consequence of most NL1 N Biogeosciences Discuss., https://doi.org /10.5194/bg-2018-232 Manuscript under review for journal Biogeosciences Discussion started: 4 June 2018 c Author(s) 2018. CC BY 4.0 License. loss now occuring as leaching (Fig. 2i). For NL2, the prediction of total N loss change switched from an increase to a decrease, on account of the leaching decrease now being of greater magnitude than the gaseous N loss increase. For NL3, precipitation increase led to strongly increased leaching and gaseous N loss reduction, approximately quadrupling the total N loss reduction compared to S1.
All models predicted NPP responses to eCO 2 of approximately 20 % (Fig. 2k) with even predictions between N 5 loss models when precipitation was increased. Model differences were also minimal regarding ecosystem C accumulation (Fig. 2l), however, all models predicted higher accumulation at S3 compared to S1 (approximately 5 %).

Local simulations II 10
While the previous section dealt with the short term effects of a ten year step-increase in atmospheric CO 2 concentrations, the second set of local experiments was designed to investigate long term effects of eCO 2 on N limitation of vegetation growth and how sensitive they were to the applied N loss formulation. During the 300 year eCO 2 simulation at the temperate "S1" site (300 years step increase from constant control CO 2 by 200 ppm), 15 all loss models predicted the total N loss rate to decrease (Fig. 3a). In the long term, the NL1 response was less pronounced than the NL2 and NL3 predictions. Note that the NL2 model produced a positive N loss response early on in the simulation (see also Sect. 3.1), but eventually predicted a negative N loss response close to the NL3 model prediction.
The ratio of total ecosystem N loss and net N mineralization (termed "N cycle openness") was predicted to 20 decline by all models (Fig. 3b). This meant that in all simulations, relative to the control runs, less N was lost from the system compared to new N becoming available from mineralization, i.e. the internal ecosystem N cycle became more "closed". As mentioned in Sect. 3.1.1, net N mineralization generally increased in the O-CN TBM under eCO 2 . For the models that predicted reduction of N loss at S1 in the shorter term (NL1, NL3; Fig. 2b), reduced N cycle openness was therefore an expected result, notwithstanding the slightly different experimental 25 designs. However, using the NL2 model that calculated ecosystem N loss increase early on resulted in reduced N cycle openness as well, meaning that early N mineralization increased at a higher rate than N loss did.
Model predictions of N cycle openness responses did differ in that NL1 predicted an inital sharper decline relative to the control runs than NL2 and NL3 did (Fig. 3b). This was a consequence of NL1 N loss being largely dependent on the soil inorganic N concentration that declined quickly in response to the step increase of the 30 atmospheric CO 2 concentration (Fig. 3a). In contrast, the dependencies of NL2 and NL3 on the N mineralization flux made for a more gradual decline of N loss (and N cycle openness), owing to N mineralization depending on the slower dynamics of the SOM N pool. Over the 300 simulation years, all loss models approached the same approximate absolute magnitude of N cycle openness reduction.
The magnitudes of NPP responses (Fig. 3b) were again largely unaffected by N loss differences, which was 35 expected considering the findings in Sect. 3.1.1 (Fig. 2c). There was a trend of NL1 sustaining a larger NPP response than the other models until simulation year 175, however, this trend dissolved over the following decades of simulation until the end of the experiment.

Global simulations
Having examined the eCO 2 effects of N loss variety in detail at a temperate site, we next applied the three loss algorithms in global simulations to observe the dynamics between eCO 2 , N loss and ecosystem C accumulation for different vegetation types and climate regimes. Notably, only the atmospheric CO 2 concentration was varied, 5 whereas climate and N forcing was fixed (see Sect. 2.2.3).
The three model versions were spun-up to quasi-equilibrium for the soil and vegetation C and N pools using 1850 atmospheric CO 2 levels according to the RCP8.5 scenario (285 ppm; Meinshausen et al., 2011). Yet, after 138 simulations years and having gradually reached atmospheric CO 2 levels of 350 ppm, global N loss rates were still predicted to be similar between the models (Fig. 4a,b,c). This showed that the models did not differ 10 much when responses to gradual, low magnitude eCO 2 were calculated. The hotspots of N loss were regions with high density of agricultural land use and, to a lesser extent, the tropical zone with high natural N turnover.
All loss models tended to predict more or less pronounced reductions of total ecosystem N loss in global simulations under eCO 2 by the time a 550 ppm atmospheric CO 2 concentration was reached, i.e. after 64 more years of simulation (Fig. 4d,e,f). Model predictions differed, however, in the magnitudes of N loss reductions 15 with some notable regional disagreement. Some of the regions for which all models consistently predicted sizable reductions of loss rates were arid parts of the Canadian Prairies, most northern temperate and boreal regions of Russia, as well as regions surrounding the Central Asian deserts, where vegetation cover and baseline N turnover was low, therefore not contributing much to the total global N loss flux in absolute terms. The models also predicted N loss reduction for tropical rainforests, regions with some of the highest global N stocks 20 and turnover rates. The most notable model differences could be found between the NL1 and NL2 loss models, with NL1 predicting generally large (often greater than 30%) negative N loss responses, and NL2 predicting generally smaller (mostly lower than 20%) negative responses (Fig. 4d,e,g). Approximately, the responses predicted by NL3 could be classified as close to NL2 in the temperate and boreal regions, and close to NL1 in the tropics. The large 25 negative response for NL1 in the boreal regions was a manifestation of the soil N concentration-based N loss fluxes. The boreal regions are usually considered strongly N limited in their vegetation growth, i.e. low on N available in the soil for plant uptake. Therefore, a small absolute decrease in soil N concentration due to plant uptake increase under eCO 2 was enough to result in a high relative reduction of total N loss in NL1. This mechanism did not apply as strictly in the NL2 and NL3 models, which resulted in less pronounced responses in 30 the boreal regions. NL2 predicted the smallest decrease in tropical N loss rates because its loss pathways were the least affected by the soil inorganic N concentration decreasing under eCO 2 , as NL2 used this dependence only to determine leaching. The model differences regarding N loss responses were most prominent in highly N limited regions (Fig. 4g).
Some ecosystems in these regions such as boreal forests are also known to store large amounts of C. This raised 35 the question whether model differences in global N loss responses, including in highly N limited regions, also resulted in appreciable model differences with respect to predicted C accumulation under eCO 2 .
As the N loss algorithms predicted different rates of N loss change in response to global atmospheric CO 2 increase (Fig. 4), there was also disagreement on the amounts of N that would accumulate in ecosystem over the entire global simulation period (1850-2100) when atmospheric CO 2 concentrations increased from 285 ppm to 40 936 ppm (Fig. 5). The predicted N accumulation varied between 1549.7 Tg N (NL1; corresponds to 6.2 Tg N yr - The resulting model disagreement regarding the additional C accumulation under eCO 2 was small (Fig. 5). Using the three N loss algorithms led to calculated cumulative terrestrial C uptake of 135.0 Pg C (NL1; corresponds to 0.54 Pg C yr -1 ), 142.3 Pg C (NL2; 0.57 Pg C yr -1 ), and 139.8 Pg C (NL3; 0.56 Pg C yr -1 ) in our eCO 2 scenario. In particular, there were only small model differences simulated for tropical or boreal forests, where the variety in 10 N accumulation was highest.

Discussion
The variability between N loss algorithms in predicted C accumulation for the RCP8.5 eCO 2 scenario between 15 1850 and 2100 (Fig. 5) was very low compared to the large variability in predictions of different TBMs for a similar scenario (Jones et al., 2013). This indicated that uncertainty in N loss representation was not a major driver for variability in future C sink predictions. This result was obtained in spite of some non-negligible variety in predicted ecosystem N accumulation during the global eCO 2 experiment (Fig. 5). The tendency of global results to indicate that variety in predictions of N loss change (Fig. 4) and N accumulation (Fig. 5) under eCO 2 20 did not have large impact on corresponding predictions of responses in C sequestration was in principle also found at the site level (Fig. 2, 3). We expected that the effect would be limited, because it was shown before that in the O-CN TBM, the magnitude of N loss is about 20 % of the magnitude of plant N uptake (Meyerholt et al., 2016). However, the small margin in C predictions between the N loss algorithms is still remarkable, considering the different levels of complexity with which the loss fluxes were determined. 25 The lack of direct connection regarding model variety between N predictions and C predictions appears plausible for the tropical zone, where vegetation growth is typically not considered N-limited. Therefore, N variety in loss rates and accumulation were not expected to affect C predictions strongly. The reasons for the small C effect outside the tropics are, however, less clear. The phenomenon might be explained by the concept of flexible C:N stoichiometry in organic soil and plant tissues employed in the O-CN TBM (Zaehle and Friend, 2010;Meyerholt 30 and Zaehle, 2015). The buffering capacity of flexible ecosystem C:N ratios could indeed attenuate the variety in N loss responses to eCO 2 and render the effect on C accumulation minimal as seen in our results (Fig. 5).
However, we found that when we employed fixed ecosystem C:N ratios in both organic soil and vegetation following , predicted C accumulation only became more variable in far northern latitudes due to variable productivity under strong N limitation, while becoming even more uniform everywhere 35 else (Fig. A1). The variety of global C accumulation predictions was hardly affected (143.9, 135, 138 Pg C for NL1, NL2, NL3; 1850 -2099). Thus, with the exact mechanisms remaining difficult to discern, we report that Aside from the short term predictions of the NL2 loss algorithm, this reduction effect (Figs. 2, 3, 4) was in line with the expected mechanisms in play for C-N TBMs under eCO 2 (Zaehle et al., 2014;Walker et al., 2015), with increased plant N demand causing more N uptake from the soil, leaving less residual N to be lost. However, although denitrification is being considered the dominant mechanism of ecosystem NO 3 loss (Fang et al., 2015), most field experiments have not found decreases in N 2 O emissions under eCO 2 . Positive or neutral responses 5 have been common, primarily obtained in temperate or boreal forests or grasslands (Van Groeningen et al., 2011;Dijkstra et al., 2012). It should be noted that none of these field experiments were conducted in the tropics, which may hamper comparisons (Huang and Gerber, 2015), also seeing how a substantial portion of our obtained N loss decrease was observed in these latitudes (Fig. 4). Nevertheless, a number of mechanisms have been proposed that would explain increased denitrification under eCO 2 (Butterbach-Bahl and Dannenmann, 10 2011), some of them on the level of abstraction with which the N cycle is commonly represented in TBMs. For example, eCO 2 could change the plant's water use efficiency through decreased plant transpiration, leading to higher soil water content and a higher likelihood of anaerobic soil conditions that would benefit denitrification.
Further, eCO 2 could stimulate the rhizodeposition of labile C compounds such as amino acids and sugars, increasing soil microbial activity and thereby providing benefitial conditions for denitrification through oxygen 15 depletion and provision of organic C. While such processes were in principle represented in the NL1 formulation (see Appendix B), the predicted result was still the reduction of N loss under eCO 2 , because the trends were stronger determined by the increased plant N uptake and soil mineral N depletion. It might be that these main mechanism proposed by models to reduce N loss under eCO 2 just do not apply as generally, especially in very N rich soils. This leads to the challenge for models to treat the actual size and distribution of global soil N 20 inventories, if modellers hope to draw connections between soil N content and soil N emissions on large spatial scales.
As for other model studies, Huang and Gerber (2015) only found initial reduction of soil N 2 O emissions in the tropical biome in a long-term eCO 2 simulation using a TBM, followed by a substantial increase and the responses in other ecosystems being neutral or positive throughout. Aber et al. (2002)

demonstrated that in a 25
stand-scale ecosystem model, the N loss (leaching only) response was only negative under simulated eCO 2 when the experiment was not confounded by other perturbations such as increased N deposition and climatic change.
In that model setting, other perturbations of N input and accelerated N turnover would eventually increase N losses, a result that we also largely obtained globally, when we added increased N deposition and transient climate to the eCO 2 experiment (Fig. A2). While there is still discrepancy between the mechanisms that likely 30 control the N loss response to eCO 2 in nature and the mechanisms that shape model responses, the most immediate notion here is that this effect needs to be studied in actual tropical ecosystems that are N rich and will be crucial for the global climate under future change. Also, the above comparisons are limited by differences in the typical durations of field campaigns and model simulation runs. Time scale may well be pivotal here, since the functioning of the N cycle and its sensitivity to changing climate and biogeochemistry has long been 35 hypothesized to change over longer (decadal and onwards) time scales (Aber et al., 1989;Vitousek and Howarth, 1991;Luo et al., 2004).
Performing a local temperate 300 year eCO 2 experiment, we found that initial model differences in N loss rate responses would over time approximately converge to a similar level that remained negative, i.e. N loss reduction (Fig. 3). The persistent N loss reduction over time meant that using our array of N loss algorithms 40 within the framework of the O-CN TBM never resulted in a prediction of long-term progressive N limitation Biogeosciences Discuss., https://doi.org /10.5194/bg-2018-232 Manuscript under review for journal Biogeosciences Discussion started: 4 June 2018 c Author(s) 2018. CC BY 4.0 License. (Luo et al., 2004) at the temperate site. Walker et al. (2015) conducted similar local model experiments at temperate sites, but instead of comparing N loss algorithms they compared entire TBMs. They found a variety of N loss responses, from negligible responses to initial reductions that would over time subside and approach zero or even positive responses, which suggested that some TBMs predicted onset of progressive N limitation under a long term eCO 2 regime. The different N loss algorithms applied in these TBMs (some of which were represented 5 in our study) were likely influential in producing the variety of eCO 2 responses presented in Walker et al. (2015).
However, we showed here that variety in N loss formulations alone was not enough to produce such heterogeneous responses in the long term. Rather, the results produced by different TBMs were confounded by many other assumptions about the C and N cycles that were inherent to each model. We have shown that the different loss algorithms simulated variable partitioning of N losses between gaseous 10 and leaching losses, both in quasi-equilibrium and for eCO 2 responses (Fig. 2). In reality, the partitioning of the two N loss pathways means no less than the difference between predicted air pollution or water pollution if N enriched ecosystems are considered (Aber et al., 1989). Consequently, adequate N cycle representation still mandates a better grasp of this partitioning issue, for which major model discrepancies have been shown before (Thomas et al., 2013). Some of the model differences were fairly obvious consequences of the respective 15 formulations, such as a very high leaching portion in a high precipitation environment when leaching was a function of drainage (NL1), or the virtual elimination of leaching in the hierarchical structure of NL3 when gaseous N loss was already substantial (Fig. 2). These examples showed that in TBMs that include both gaseous and leaching pathways, there has not been a consensus on the proper partitioning between fluxes that needs to be reflected. The reason for this is surely the lack of field evidence of simultaneous measurements of both pathways 20 to inform models. Some models not considered here have included leaching only (Gerber et al., 2010;Wania et al., 2012; note that the Gerber et al. formulation was updated in Huang and Gerber, 2015), which further illustrates that across current TBMs, N loss representation and in particular the partitioning between gaseous and leaching losses has not been an issue of priority. Our results that indicated little effect of the N loss representation on future C sequestration under eCO 2 should not encourage continued neglect of this aspect of N 25 cycle modelling (Walker et al., 2015).
Further, deriving gaseous N loss from the N mineralization flux as described e.g. by Thornton and Rosenbloom (2005) and Wang et al. (2010) may be too coarse of an approximation. The "Hole-in-the-pipe" model (Firestone and Davidson, 1989), frequently cited as inspiration, highlighted that the fluxes of nitrification and denitrification are the important processes from which NO and N 2 O emission are derived. At best, mineralization 30 could be considered as a "distal" indicator of nitrification activity, since it aids in providing NH 4 + substrate, one of the controlling factors of nitrification. From the ecosystem modelling perspective, denitrification is usually assumed to be controlled by the presence of NO 3 -, organic C and anaerobic conditions, not N mineralization.
Likely, the relatively simple algorithms we implemented as NL2 and NL3 were in part the result of the respective N-cycle model not discriminating between NO 3 and NH 4 + and rather calculating a generic N species 35 for convenience. This discrimination will likely be necessary if N cycle modelling and N loss algorithms in particular are to be advanced across terrestrial biosphere models.
The above example illustrates that regional scale modelling of N loss fluxes in TBMs is still developing, as current research continues to e.g. investigate the sensitivity of leaching to global change in a TBM (Braakhekke et al., 2017). Importantly, one of the most climate-relevant aspects of N loss fluxes, N 2 O emission, has become 40 the focus of a model intercomparison project aiming to understand past and present N 2 O fluxes by utilizing state-Biogeosciences Discuss., https://doi.org /10.5194/bg-2018-232 Manuscript under review for journal Biogeosciences Discussion started: 4 June 2018 c Author(s) 2018. CC BY 4.0 License.
of-the-art models along with available data (Tian et al., 2018). We contend that such efforts to consolidate the representation of ecosystem N loss processes could best be aided by field experiments that investigate N loss rates (and their partitioning between gaseous and leaching components) under perturbation in the regions that we identified as crucial with respect to modelled N loss uncertainty. Our global simulations showed consistent model disagreement in northern temperate and boreal latitudes on the magnitude of N loss rate changes under 5 eCO 2 and a small effect on C sequestration. Vegetation growth in these regions is usually thought to be strongly N limited, which is largely controlled by the ecosystem N budgets, including N loss rates. For future C sink predictions, it might be important to describe how N cycling affects the mobility of the large C reserves in these latitudes that are often considered undersampled with respect to many ecosystem variables. On the other hand, we found large N loss rate changes in the tropics, but also faced a lack of appropriate field experiments to 10 evaluate these results. While tropical vegetation growth is usually not considered N limited, the flipside is that many of these regions are N rich, with large potential of water eutrophication or the outgassing of NO and N 2 O, all environmental issues of note where new experiments are urgently needed to inform models.

Code availability
The used model code is available from the authors upon request.

Author contribution
JM and SZ designed the study. JM and SZ set up the simulations. JM ran the local experiments, SZ ran the 20 global experiments. JM and SZ analyzed the results. JM wrote the manuscript.

Figure 2:
Average fate of ecosystem nitrogen (N) input at the three pseudosites "S1", "S2", and "S3" for the 5 1700-1750 quasi-equilibrium simulation period without perturbation (a, e, i; g N m -2 yr -1 ). Average N loss rate responses (g N m -2 yr -1 ) to ten years of simulated elevated atmospheric CO 2 (eCO 2 , +200 ppm) (b, f, j). Net primary productivity (NPP) responses (%) after ten years of eCO 2 , relative to control simulations (c, g, k). Change (%) of total ecosystem carbon (C) after ten years of eCO 2 , relative to the start of the experiment (d, h, l). Colours indicate the applied N loss algorithm ("NL1", "NL2", "NL3"). For a, b, e, f, i, and j, shading lines 10 indicate leaching, unshaded colour indicates gaseous loss. Numbers on the bars indicate the standard deviation (g N m -2 yr -1 ) of the leaching and gas loss components over the 1700-1750 quasi-equilibrium period. For a, e, and i, the sums of N allocation to organic pools (vegetation biomass and soil organic matter) were in the range of -0.07 -0.01 g N m -2 yr -1 and were omitted here. External inputs of reactive N (biological N fixation + N deposition) were fixed at 1.631 g N m -2 yr -1 . For b, f, and j, black bars indicate the total N loss responses (gaseous loss + 15 leaching).

Appendix B: Nitrogen loss algorithm description for NL1
Gaseous nitrogen losses NL1 includes the adsorption of ammonium (NH 4 + ) ions to clay particles, following Li et al. (1992). Thereby, the 5 NH 4 + fraction of the soil pool of available inorganic nitrogen (N) is first reduced according to the soil clay content. The anaerobic volume fraction of the soil (anvf; Eq. B1) is estimated from the fractional soil moisture content (θ) by an empirical function: In the aerobic part of the soil (1-anvf), the fraction of ammonium (NH 4, aerob ) is subject to nitrification, with the 10 gross rate (N nit ; Eq. B2) depending on NH 4 + concentrations, as well as response functions for temperature (f nit (T)) and soil pH (g nit (pH)): where α is chosen such that at 20°C and favourable pH conditions, 10% of the ammonium is nitrified per day.
The temperature function (Eq. B3) is taken from Xu-Ri and Prentice (2008) where T soil is the soil temperature in °C. The soil pH function (Eq. B4) is taken from Zhang et al. (2002): N 2 O loss from nitrification (N 2 O nit ; Eq. B5) is estimated from the gross nitrification rate : . (B6) Loss of NO from nitrification (Eq. B7) follows the formulation by  with an additional component representing chemonitrification following Kesik et al. (2005): where β (Eq. B11) is a function that describes the regulatory influence of soil microbial activity and NO 3-5 concentrations on gross denitrification : where R mb is the microbial respiration rate, NO 3 is the soil NO 3 concentration (aerobic and anaerobic), and K R and K NO3 are half-saturation constants. The temperature sensitivity of denitrification (f denit (T); ; Eq. B12) is taken from Xu-Ri and Prentice (2008) ) . (B12) The sensitivity of denitrification to soil pH is described by g denit (pH) (Eq. B13): Gaseous N losses from denitrification (NO denit , Eq. B14; N 2 O denit , Eq. B15; N 2,denit , Eq. B16) are then estimated from the gross denitrification rate, taking into account the different sensitivities to soil pH for the respective 15 proportions of emissions of NO, N 2 O, and N 2 : where β NO and β N2O are constants and g denit,NO (Eq. B17) and g denit,N2O (Eq. B18) are functions that scale emission 20 of different gaseous N compounds with soil pH: with the diffusion coefficient d ox that depends on soil moisture and also determines volatilization losses from the soil N pool of other N gases: Nitrogen leaching Leaching of NH 4 + and NO 3 occurs in proportion to the water lost from soil drainage, calculated as described by 10 De Rosnay and Polcher (1998).  Anaerobic volume fraction of the soil (Eq. B1) Fractional soil moisture content (Eq. B1) NH 4 + in the aerobic fraction of the soil N pool (Eq. B2) Gross N nitrification rate (Eq. B2) Factor to scale nitrification activity (Eq. B2) Temperature response function of nitrification (Eq. B2) pH response function of nitrification (Eq. B2) Soil temperature N 2 O emission from nitrification (Eq. B5) Temperature function for N 2 O emission (Eq. B6) NO emission from nitrification (Eq. B7) Temperature function for chemonitrification (Eq. B8) pH function for chemonitrification (Eq. B9) Gross denitrification rate (Eq. B10) Microbe function of gross denitrification (Eq. B11) Temperature function for denitrification (Eq. B12) pH function for denitrification (Eq. B13) NO 3 in anaerobic fraction of the soil N pool (Eq. B10) Microbial respiration rate (Eq. B11) Half-saturation constant (Eq. B11) Half-saturation constant (Eq. B11) NO loss from denitrification (Eq. B14) N 2 O loss from denitrification (Eq. B15) N 2 loss from denitrification (Eq. B16) Constant (Eq. B14) Constant (Eq. B15) pH sensitivity function for NO denitrification (Eq. B17) pH sensitivity function for N 2 O denitrification (Eq. B18) pH sensitivity function for NO denitrification (Eq. B19) pH sensitivity function for N 2 O denitrification (Eq. B20) Volatilization of NH 3 from the soil N pool (Eq. B21) NH 4 + concentration in the soil N pool (Eq. B21) Soil moisture dependent diffusion coefficient (Eq. B22) Volatilization of NO from the soil N pool (Eq. B22) Volatilization of N 2 O from the soil N pool (Eq. B23) Volatilization of N 2 from the soil N pool (Eq. B24) NO concentration in the soil N pool (Eq. B22) N 2 O concentration in the soil N pool (Eq. B23) N 2 concentration in the soil N pool (Eq. B24) -g N m -2 g N m -2 dt 1.2 --°C g N m -2 dt g N m -2 dt -g N m -2 dt --g N m -2 dt dt g N m -2 g N m -2 dt g N m -2 dt g N m -2 dt 0.78 0.54 ---g N m -2 dt g N m -2 0.001-0.005 dt g N m -2 dt g N m -2 dt g N m -2 dt g N m -2 g N m -2 g N m -2 5 10 Biogeosciences Discuss., https://doi.org /10.5194/bg-2018-232 Manuscript under review for journal Biogeosciences Discussion started: 4 June 2018 c Author(s) 2018. CC BY 4.0 License.