Variability of projected terrestrial biosphere responses to elevated levels of atmospheric CO 2 due to uncertainty in biological nitrogen fixation

Including a terrestrial nitrogen (N) cycle in Earth system models has led to substantial attenuation of predicted biosphere–climate feedbacks. However, the magnitude of this attenuation remains uncertain. A particularly important but highly uncertain process is biological nitrogen fixation (BNF), which is the largest natural input of N to land ecosystems globally. In order to quantify this uncertainty and estimate likely effects on terrestrial biosphere dynamics, we applied six alternative formulations of BNF spanning the range of process formulations in current stateof-the-art biosphere models within a common framework, the O-CN model: a global map of static BNF rates, two empirical relationships between BNF and other ecosystem variables (net primary productivity and evapotranspiration), two process-oriented formulations based on plant N status, and an optimality-based approach. We examined the resulting differences in model predictions under ambient and elevated atmospheric [CO2] and found that the predicted global BNF rates and their spatial distribution for contemporary conditions were broadly comparable, ranging from 108 to 148 Tg N yr (median: 128 Tg N yr), despite distinct regional patterns associated with the assumptions of each approach. Notwithstanding, model responses in BNF rates to elevated levels of atmospheric [CO2] (+200 ppm) ranged between −4 Tg N yr (−3 %) and 56 Tg N yr (+42 %) (median: 7 Tg N yr (+8 %)). As a consequence, future projections of global ecosystem carbon (C) storage (+281 to +353 Pg C, or+13 to+16 %) as well as N2O emission (−1.6 to +0.5 Tg N yr, or −19 to +7 %) differed significantly across the different model formulations. Our results emphasize the importance of better understanding the nature and magnitude of BNF responses to change-induced perturbations, particularly through new empirical perturbation experiments and improved model representation.

To assess this uncertainty, we tested six alternative approaches to represent BNF embedded 1 within the framework of a common TBM, the O-CN model (Zaehle and Friend, 2010), which 2 comprises a comprehensive description of the terrestrial C and N cycles and their interactions 3 with the terrestrial energy and water balance. Applying all BNF schemes directly in a full 4 TBM allowed us to appraise the consequences of uncertainty in BNF representations for the 5 simulated C cycle. The BNF models included a prescribed global map of static BNF rates, 6 two simple empirical relationships between BNF and other ecosystem variables (NPP and 7 ET), two formulations based on plant N status, and an approach following a basic form of 8 optimality of plant N acquisition (Table 1). 9 We first applied these alternative BNF model versions of O-CN to simulate the pre-industrial 10 to present-day global patterns of the terrestrial C and N cycle to analyze the implied spatial 11 patterns of BNF and associated projected C and N fluxes. We then sought to test the implied 12 sensitivity of BNF, and thus the coupled C-N cycles, to changes in N limitation. We did this 13 by driving the model versions with idealized transient and step-wise eCO 2 scenarios to make 14 the functional model differences clearly apparent. The increased C availability increased plant 15 N demand, and this demand was met with a variety of approaches to determine the ecosystem 16 N input of BNF, which emphasized the different characteristics of the alternative approaches. 17 In particular, we expected a pronounced discrepancy between empirical and mechanistic BNF 18 representations, highlighting a previously unquantified source of variation in the predictions 19 of C-N terrestrial biosphere models.    8 We conducted simulations applying six alternative models of symbiotic BNF currently 9 applied in TBMs, which are described in Sect. 2.2.1 to 2.2.6 (Zaehle and Dalmonech, 2011; 10 Table 1; Appendix). Conceptually, the BNF models can be summarized as model forcing  availability), BNF in this approach is set to converge towards zero when soil inorganic N 27 concentrations exceed 2 g N m -2 . Thus, average BNF rates still vary due to any mechanics that 28 affect the soil N pool, such as seasonal variations in plant N uptake and organic matter 29 mineralization, or long-term shifts in these quantities under perturbation. Because this 30 approach does not separate between symbiotic and asymbiotic pathways, BNF in FOR is environments colder (or warmer) than 25°C, more C needs to be invested into BNF (Fisher et 1 al., 2010). The costs of root N uptake are implicitly accounted for through root turnover, 2 leading to higher uptake costs for higher investment into uptake structures (i.e. roots) to attain 3 a given rate of BNF. BNF is thus limited by the N status of the plant and its C resources.   15 The OPT model uses an optimality-based approach that follows the concept described by 16 Rastetter et al. (2001). In this model, BNF only occurs when the C cost of BNF, indicative of 17 energy (glucose) investment, is lower than the C cost of root N uptake. This cost of C 18 investment in root N uptake is evaluated as the potential plant C gain if a marginal amount of 19 C was allocated to leaves for photosynthesis, relative to the potential plant N gain if that same 20 marginal amount of C was allocated to increase fine root mass instead. This way, the C cost 21 of root N uptake is defined as the amount of C from photosynthesis the plant relinquishes in 22 favour of investment into root N uptake. If this cost is higher than the (fixed) C cost of BNF, 23 BNF occurs and is determined as a saturating function of root mass and the difference in C 24 cost between root N uptake and BNF. Notably, the occurrence and magnitude of BNF does 25 not feed back on the determination of plant root N uptake in this approach. conditions promote high photosynthetic efficiency, e.g. through high irradiation or elevated 28 atmospheric CO 2 concentrations, and increasing leaf mass is a worthwhile investment.

29
Furthermore, high plant root mass or low soil inorganic N availability will increase the C cost 30 of increasing root N uptake and consequently favor BNF. This approach has not been used in Asymbiotic BNF was calculated for the fraction of the soil receiving light, thus declining with 5 increasing light interception by the vegetation. A maximum rate of 0.2 g N m -2 yr -1 was 6 assumed based on the data presented by Cleveland et al. (1999), which was modulated by soil 7 moisture availability and soil temperature to account for reduced biochemical activity in dry, 8 cold, or hot environments. All simulation experiments were repeated for each of the six BNF models described above.

11
The aim was to elucidate the effects of the alternative representations on estimates of present-12 day BNF and its impact on terrestrial C and N cycles, as well as on projections of the 13 consequences of increasing atmospheric CO 2 concentrations, a key factor in decreasing N 14 availability over time. N deposition, land-use, and fertilizer data, as well as observed changes in atmospheric CO 2 25 concentration (A; Fig. 2). We used this simulation to evaluate the differences in estimates of 26 the global C and N cycles under present-day conditions, as described in Sect. 3.1. 27 We then evaluated the effect of eCO 2 on terrestrial C and N fluxes for the different models by 28 comparing A to a simulation with a larger increase in atmospheric CO 2 concentrations (B; 29 Fig. 2), with the other forcings as in A (Sect. 3.2). To avoid a dependency of the simulations on a specific future emission pathway under a particular scenario, we applied a monotonic 1 increase of atmospheric CO 2 from 1860 conditions (286 ppm) at a rate of 0.5% yr -1 , which 2 corresponds to an average growth rate of 2.1 ppm yr -1 , approximately comparable to the 3 currently observed growth rate of atmospheric CO 2 , arriving at 600 ppm at the end of the 4 simulation. We also compared B to a simulation with CO 2 fixed at 1860 conditions (286 ppm, 5 C) to elucidate the cumulative effect of eCO 2 on the time evolution of key ecosystem fluxes 6 and stocks of C and N. of asymbiotic BNF in total BNF between 1.0% (NDS) and 1.4% (OPT).
Notwithstanding, individual BNF models differed considerably in their predictions in many 1 regions (Fig. 3b). In Europe, the eastern US, East Asia, and extratropical South America, the 2 empirical models (AET, PRO) predicted higher BNF rates than the other approaches. In these 3 regions with wide-spread human activity, fertilizer application and atmospheric N deposition 4 caused high N availability for plants, which either directly reduced BNF (FOR, OPT), or over 5 time diminished the plants' N demand and thereby BNF (NDT, NDS). These mechanisms did 6 not apply in the empirical models. Another important model difference is the large 7 discrepancy in simulated BNF in northern Russia and Canada (Fig. 3b) that mainly stems 8 from very high BNF rates predicted by the N demand-based models (NDT, NDS). In both 9 approaches, strong N limitation in these regions increased BNF beyond plausible rates 10 (Cleveland et al., 1999), occasionally in excess of 3 g N m -2 yr -1 in the case of NDS (Fig. 4b).

11
The lack of temperature control on BNF in NDS resulted in notably higher predicted BNF 12 rates in the boreal zone than in NDT, which led to substantial alleviation of N limitation (Figs.
14 All models simulated the highest cumulative BNF rates for tropical forests and global 15 grasslands (Fig. 4) The model uncertainty in BNF did not cause large uncertainty in the predicted global gross 1 and net primary productivity (GPP and NPP; Table 2). Notably, the inclusion of respiration 2 costs of BNF in NDT, NDS, and OPT did not result in a significant reduction in C-use 3 efficiency, potentially because of the reduced severity of N limitation, which reduced excess 4 respiration. The spatial patterns of simulated rates of NPP were also very similar for large with high simulated NPP, the differences between models in BNF barely affected NPP.

15
The between-model difference in N input rates was, however, reflected in the other branches

27
We next analyzed the effect of increasing N stress through CO 2 fertilization by comparing the 28 final 13 years of the simulations B and A (Fig. 5). For an average atmospheric CO 2 29 concentration difference of 211 ppm, the predicted total global BNF response to eCO 2 ranged 30 between a 4 Tg N yr -1 reduction (AET) and an increase of 56 Tg N yr -1 (NDS) (median 31 increase of 7 Tg N yr -1 ), corresponding to -4 and 38 % (median 6%) of the average BNF rates under ambient CO 2 (Fig. 3a), respectively. The median predicted responses of global BNF 1 rates to eCO 2 ( Fig. 5a and b) indicated a substantial increase in N fixation in many regions. In 2 the N-demand based approaches, increased C availability increased global plant N demand, 3 having a strong relative effect in boreal and northern temperate regions that were already 4 strongly N limited (Figs. 5b and B3). The eCO 2 experiment also resulted in predicted global 5 NPP increases ( Fig. 5c and d). The predictions ranged between 15 and 21 Pg C yr -1 (median 6 17 Pg C yr -1 ), with all models simulating the highest NPP increases in the tropics (Fig. B4).

7
The increase in BNF rates in responses to eCO 2 was by far strongest in the N-demand based 8 models (Fig. 6). The increased C fixation under eCO 2 temporarily increased the simulated 9 labile reserve of allocatable C, which in NDT was directly connected to predicted BNF rates. was driven by the reduction of stomatal conductance in response to eCO 2 . In OPT, eCO 2 led 15 to more efficient photosynthesis, which reduced C allocation to roots for N uptake and 16 thereby increased global BNF rates moderately.

17
The above variation between models in BNF response magnitudes did not translate into  Table 2) also predicting high C sequestration. These ecosystem C 29 storage responses correspond to a range of C-concentration interactions in the sense of derived from these studies are not comparable, because the increment of gradual CO 2 increase 1 was only half in our study compared to Gregory et al. (2009).

2
The choice of BNF model also had substantial effects on other quantities relevant for 3 biogeochemistry-climate effects, in particular the predicted responses of N 2 O emissions to 4 eCO 2 (Fig. 7d). In the larger group of models suggesting moderate changes in global and 5 regional BNF, global N 2 O emission rates were simulated to decrease with eCO 2 . With

Discussion
Given the large variation in approaches used to calculate BNF in this study, ranging from 1 empirical correlation to process-oriented models, our simulations resulted in surprisingly 2 similar estimates of BNF for the contemporary period over large parts of the terrestrial 3 biosphere, despite very notable regional differences. The predicted range of global present-4 day BNF rates of 108-148 Tg N yr -1 compared reasonably well with the conservative end of One of the prominent regions for which simulated BNF was highly uncertain were high-10 latitude ecosystems (Fig. 3). Open vegetation in these ecosystems contributed to very high 11 BNF in the NDS scheme in boreal forests and grasslands ( Fig. 4b), which made this scheme 12 distinct from the others in this region. We also found a strong heterogeneity of predicted BNF Attempting to incorporate process hypotheses rather than empirical relationships is expedient 21 and also led to lower N losses relative to ecosystem N accumulation in comparison with other 22 approaches (Table 2), which heuristically appears to be more plausible. Yet, the behaviour of 23 the plant N status-based models NDT and NDS was likely implausible in other aspects, 24 particularly the strong, quasi-instantaneous increase of BNF under the scenario of a step-25 increase in atmospheric CO 2 (Fig. 7).

12
The effect of the alternative BNF process representations was significant also for predictions 13 on other contemporary key N fluxes (   contemporary NPP is further explained by the fact that despite regional differences in N 2 limitation evidenced by moderate regional differences in foliar stoichiometry, on global 3 average, the simulated vegetation growth was not strongly N limited for any BNF approach and maximal values (Table A2) and approximately similar between BNF approaches (average 8 global ratios between 30 (AET) and 33 (NDS, OPT)). 9 Unlike the small effect under contemporary conditions, the uncertainty in predicted BNF rates 10 under eCO 2 had a sizeable effect on the predicted NPP and C sequestration, resulting from the 11 differences in gradual ecosystem N accumulation (Fig. 7). The ecosystem N input from BNF  We contend that improving the representation of BNF in TBMs will be greatly aided by a   Table A1. where j is the fraction of C labile available for investment into BNF (as C labile also contains the 4 assimilated C available for allocation to plant growth). The ξ function sets C inv to zero at 5 extreme temperatures: The η function scales C inv with the plants' N status, represented by their leaf C:N ratios:  where C Leaf is the leaf C pool size and BNF L is the BNF rate per unit leaf C, described in 17 differential form: where SLA is the specific leaf area. The establishment of BNF is controlled by the plants' 1 local N demand ψ per unit leaf C, which in turn is determined by the plant N deficit (D) and a 2 function (κ) that scales the advantageousness of BNF with the plants' N status: We define D as the difference between the N that is required to build new biomass from 5 newly acquired C and the N that is available to the plant for allocation to new biomass: where NPP pot is the allocatable C after respiration costs are satisfied, f cost is a dimensionless To determine the instantaneous C gain per unit leaf area (k), we consider the relationship of 1 gross primary productivity (GPP) and the fraction of absorbed photosynthetically active 2 radiation, which depends on the specific leaf area and leaf mass: We then derive the marginal C gain with C investment into leaves, gc, from the difference in as the instantaneous C Root -specific N uptake: 11 We then evaluate the C cost of N uptake (r Nup ) as: If r Nup is larger than the C cost of BNF (r Fix , assumed constant), BNF is calculated as a 14 saturating function of (r Nup -r Fix ) and root mass: 16 where v max,Fix is a maximum BNF rate and k Fix is a half-saturation constant. In case the C cost 17 of BNF is higher than the cost of root N uptake, no symbiotic BNF occurs. The asymbiotic BNF rate scales with the same temperature function applied in the NDT 21 approach, but rather than the surface temperature, the function ts involves the soil temperature Asymbiotic BNF is only calculated for the fraction of the soil surface receiving solar energy. 1 We consider light limitation by applying the simple shading function vf, causing BNF to 2 converge towards zero with canopy closure: where SLA is the specific leaf area of the respective PFT and C Leaf is the leaf C pool size.

5
Also, the limiting effect of drought conditions on heterotrophic BNF is taken into account by 6 including the soil moisture function Φ: where σ is the current amount of water stored in the soil, z is the total depth of the soil 9 reservoir, and σ max is the amount of water stored in a water saturated soil column. The     PFT-specific parameters are given in Table A2.   Table A2. PFT-specific parameters. The CN parameters were used in all models, the λ 0 and σ 1 parameters were used in the NDS model (see Table A1). The PFT classes are defined in Table   2 B1.     For the modeled BNF rates, markers indicate the mean value over all grid cells that included 28 the respective biome type, error bars indicate the corresponding standard deviation. The black line is the one-to-one line. Details on the classification of vegetation types from the data 1 source into the plant functional types applied in O-CN can be found in Table B1. Relative model-median NPP responses. Figures B3 and B4 provide BNF and NPP maps for 10 each model separately.  (g N m -2 yr -1 ) to variation in the root N uptake cost r Nup (g C g -1 N) and the half-saturation 20 constant k Fix (g C g -1 N) according to Eq. A19. C root was fixed at 200 g C m -2 , v max,Fix was 21 fixed at 0.0225 g N g -1 C yr -1 , and r Fix was fixed at 9 g C g -1 N. The arrow indicates that BNF 22 is zero when r Nup =r Fix , therefore variation in r Fix would shift the functions in x-direction. (c)

23
OPT: Sensitivity of BNF to variation in r Nup and the maximum BNF per unit root C, v max,Fix , 24 according to Eq. A19. C root was fixed at 200 g C m -2 , k Fix was fixed at 50 g C g -1 N, and r Fix 25 was fixed at 9 g C g -1 N.