Assessment of the Impacts of Biological Nitrogen Fixation Structural Uncertainty in CMIP6 Earth System Models

. Biological nitrogen fixation is the main source of new nitrogen into natural terrestrial ecosystems and consequently in the nitrogen cycle in many earth system models. Representation of biological nitrogen fixation varies, and 10 because of the tight coupling between the carbon and nitrogen cycles, previous studies have shown that this affects projected changes in net primary productivity. Here we present the first assessment of the performance of biological nitrogen fixation in models contributing to CMIP6 compared to observed and observation-constrained estimates of biological nitrogen fixation. We find that 9/10 models represent global total biological nitrogen fixation within the uncertainty of recent global estimates. However, 6/10 models overestimate the amount of fixation in the tropics, and therefore the extent of the latitudinal 15 gradient in the global distribution. For the SSP3-7.0 scenario of future climate change, models project increases in fixation over the 21 st century of up to 80%. However, while the historical range of biological nitrogen fixation amongst models is large (up to 140 kg N ha -1 yr -1 at the grid cell level and 43 - 208 Tg N yr -1 globally) this does not have explanatory power for variations within the model ensemble of net primary productivity or the coupled nitrogen-carbon cycle. Models with shared structures can have significant variations in both biological nitrogen fixation and other parts of the nitrogen cycle without 20 differing in their net primary productivity. This points to systematic challenges in the representation of carbon-nitrogen model structures, and the severe limitations of models using net primary productivity or evapotranspiration to project the biological nitrogen fixation response to elevated atmospheric carbon dioxide or other environmental changes.


Introduction
In a key innovation to CMIP5, the majority of earth system models (ESMs) of the latest generation that contribute to CMIP6 25 (Eyring et al., 2016) include a nitrogen cycle to better represent the terrestrial carbon cycle (Arora et al., 2020;Davies-Barnard et al., 2020). Nitrogen is a key nutrient requirement for plants to plants to take up carbon, and in its bioavailable inorganic form, is highly liable to losses via gaseous and water processes (Thomas et al., 2013;Vitousek and Howarth, 1991). Over the last few decades, terrestrial carbon uptake has sequestered around a quarter of anthropogenic carbon emissions (Friedlingstein et al., 2020). However, previous assessments of ESMs have suggested that future projections of 30 terrestrial carbon storage are decreased by 37 -58% if nitrogen availability is accounted for (Wieder et al., 2015;Zaehle et al., 2014). Therefore, the accuracy of ESMs, which help guide policy on preventing further climate change, is partly determined by the functioning of the nitrogen cycles within them.
The uptake of new carbon by plants is reliant on new sources of nitrogen, as existing nitrogen may not be bioavailable. The sources of this new input of nitrogen vary by biome, including anthropogenic inputs via addition of fertiliser 70 -108 Tg yr -1 35 (Lu and Tian, 2017;Potter et al., 2010) and increased deposition, and natural sources such as lightning 3.5 -7 TgN yr -1 (Tie et al., 2002), atmospheric N deposition 63 TgN yr -1 (Lamarque et al., 2013), weathering (Holloway and Dahlgren, 2002), and biological nitrogen fixation (BNF) 40 -141 TgN yr -1 (Davies- Barnard and Friedlingstein, 2020;Vitousek et al., 2013). In many natural ecosystems BNF is likely the largest natural or anthropogenic source of new nitrogen to the terrestrial biosphere. But because of the intricate processes that control fixation, and the lack of global estimates from observations, 40 also the most uncertain (Meyerholt et al., 2016;Reed et al., 2011). Therefore, continued carbon sequestration in critical natural ecosystems that are present day and future carbon sinks is reliant on BNF. We need to know how well models are representing current the quantity and distribution of BNF to assess the reliability of the functions and therefore the robustness of future projections of terrestrial carbon uptake. Studies of individual models suggest differences in representation of BNF can lead to widely differing future terrestrial carbon sequestration (Meyerholt et al., 2016;Peng et al., 45 2 Methods

ESM Simulations
We use results from 10 ESMs: CMCC-CM2, TaiESM1, CESM2, NorESM2, UKESM1, AWI-ESM, MPI-ESM, ACCESS, 65 EC-Earth, and MIROC. The simulations used were the historical runs from CMIP6 deck simulations (Eyring et al., 2016) WCRP for the period 1950 -2014. A list of the reference ids of the simulations used can be found in the SI.

Land Surface Model Simulations
As an additional check on the performance of the ESMs, we also looked at the BNF of a number of land surface models (LSMs) in offline simulations include this in the SI. The LSMs used were CLM5, CLM4.5, JSBACH, JULES, and LPJ-70 GUESS, and are all used within ESMs considered here. These simulations and their methodology are fully described in Davies-Barnard et al., (2020).

BNF in the Models
A summary of the models' can be found in Table 1. Although there appears to be a range of approaches to BNF, every model considered here is partially or entirely based on Cleveland et al., (1999). 75

CABLE and CASACNP -Used in ACCESS
The Nitrogen cycle in the CABLE model (Law et al., 2017) of the ACCESS ESM relies on the CASACNP model, as described by Wang and Houlton, (2009);and Wang et al., (2007). Symbiotic BNF is calculated as a function of soil moisture, soil temperature, soil N availability, and NPP. Free-living BNF is calculated using biome level observational averages adapted from Cleveland et al., (1999) with a range of 0.7 -9.2 kg N ha -1 yr -1 (tropical forest highest, needleleaf 80 forest lowest) (Wang and Houlton, 2009).

CLM4.5 -Used in CMCC-CM2 and TaiESM1
The Community Land Model version 4.5 (CLM4.5; Koven et al., 2013;Oleson et al., 2010) is used in the Euro-Mediterranean Centre on Climate Change coupled climate model (CMCC-CM2; Cherchi et al., 2019) and TaiESM1 (Wang et al., 2021). The N component is described in Koven et al., (2013). 85 BNF is calculated as an exponential saturating function of NPP based on Thornton et al., (2007), which is based on Cleveland et al., (1999) with a 7 day lag to match seasonal BNF to NPP. There is no differentiation between symbiotic and free-living BNF.

CLM5 -Used in CESM2 and NorESM2
The Community Land Model version 5 (CLM5; Lawrence et al., 2019) is used in The Community Earth System Model 90 Version 2 (CESM2; Danabasoglu et al., 2020) and the Norwegian Earth System Model version 2 (NorESM2; Seland et al., 2020). CLM5 is the latest version of CLM and represents a suite of developments on top of CLM4.5. The N component is described in Fisher et al., (2010);and Shi et al., (2016).
Symbiotic BNF is calculated on a carbon cost basis for acquiring N, derived from the Fixation and Uptake of Nitrogen (FUN) approach (Fisher et al., 2010). Free-living BNF in CLM5 is calculated separately as a function of evapotranspiration 95 based on Cleveland et al., (1999) (Lawrence et al., 2019).
There is no differentiation between symbiotic and free-living BNF.
BNF is calculated as a linear function of NPP, 0.00016 kg N per kg C NPP (Wiltshire et al., 2021), based on Cleveland et al., (1999). There is no differentiation between symbiotic and free-living BNF.
BNF is a linear function of ET, 0.0102 ET (mm yr -1 ) + 0.524 (Smith et al., 2014), based on Cleveland et al., (1999). There is no differentiation between symbiotic and free-living BNF. The amount of BNF is capped at 20 kg N ha -1 yr -1 . 115

VISIT-e -Used in MIROC
VISIT-e is used in the Model for Interdisciplinary Research on Climate, Earth System version 2 for Long-term simulations (MIROC-ES2L) (Hajima et al., 2020). The nitrogen component is described in Hajima et al., (2020). Cleveland et al., (1999). Symbiotic and free-living BNF are calculated using the same function and distinguished by symbiotic BNF being directly available to plants, whereas free-living BNF is assumed to 120 be part of the litter. Symbiotic BNF represents 50% of BNF. In cropland a higher level of BNF occurs for nitrogen fixing crops, but non-fixing crops have the same BNF as natural vegetation (Hajima et al., 2020).

Observations
Following the methods of Davies-Barnard & Friedlingstein, (2020) we reviewed the BNF literature to find observational data that covered all, or close to all, BNF at a field site (i.e. including symbiotic and free-living fixation of as many BNF types as 125 are present). The locations of the site observations used can be found in SI Figure 1. Further details of the observations are in SI Table 1. Few measurements are available, with studies usually focussing on either symbiotic or free-living BNF. Since recent meta-analysis suggests that free-living is approximately a third of total BNF, and higher in some regions, we only consider data that includes explicitly both symbiotic and free-living BNF or states that all sources of BNF are measured.

Present day BNF
The majority of the models have total global BNF between 80 to 130 TgN yr -1 (Figure 1 a), within the uncertainties of two recent observation-based BNF estimates (Davies-Barnard and Friedlingstein, 2020; Vitousek et al., 2013). The range between CMCC, TaiESM1, UKESM1, MPI-ESM, and AWI-ESM which all calculate BNF from NPP is just 36 TgN yr -1 . There is, in some instances, as much variation in global total BNF within models that share components than between different models 135 (see Methods). For instance, CESM2 and NorESM2 share the same land surface model and the modelled BNF is still 43 TgN yr -1 different. However, there is little relationship between BNF function and total global BNF, with the two models using BNF based on ET encompassing the lowest and second highest values. This is suggestive of a substantial role for climate in modelling of BNF and a deliberate clustering to the most common BNF estimate (Davies-Barnard and (2020) and Vitousek et al., (2013), the majority of models predict too much BNF in the tropics. In the low latitudes (30°N to 30°S) 6 of the 10 models are above the observation-based estimate (Figure 1 c), but in the mid latitudes only 1 model is above (Figure1 d) and in the high latitudes none (Figure 1 b). The multi-model mean of BNF from CMIP6 ESMs compared to an observation-based estimate ( Figure 2) shows a broad agreement in spatial patterns, although there are clear weakness of 145 the ESMs' BNF estimates in tropical forests, where BNF is overestimated. This is to be expected, as most of the models are based on the data and linear modelling presented in Cleveland et al., (1999). However, that linear modelling has been superseded and recent studies have shown much lower BNF in tropical forests (Davies-Barnard & Friedlingstein, 2020;Sullivan et al., 2014;Vitousek et al., 2013). Although there are sources of error in the models, notably differing climate in the models' historical simulation compared to observed, these errors persist in the land surface model components of the ESMs 150 when driven with observed data (see SI Figure 2).
The pattern of high BNF in the tropics is partly due to a small number of models with very high BNF (SI Figure 3).
ACCESS has areas of anomalously high BNF in the tropics of up to 139 kg N ha -1 yr -1 . MIROC also has grid-cells of up to 193 kg N ha -1 yr -1 . Whereas in other models the tropical peak is below 41 kg N ha -1 yr -1 . While measurements of BNF from individual nitrogen fixating plants can be over 100 kg N ha -1 yr -1 , these rarely occur at a density of more than 30% cover 155 (Davies- Barnard and Friedlingstein, 2020). At the field scale BNF rarely exceeds around 15 kg N ha -1 yr -1 for free-living BNF and 20 kg N ha -1 yr -1 for symbiotic BNF (Davies- Barnard and Friedlingstein, 2020). Therefore, values much above 35 kg N ha -1 yr -1 at an ESM grid-cell level seem improbable.
Comparison of models to individual BNF field-scale observations of all BNF (free-living and symbiotic) ( Figure 3) show similar differences in latitudinal variation as the global and averaged data comparisons (Figures 1 and 2). The models 160 underestimate mid latitude wetland and peatland BNF (Massachusetts and S Germany, Figure 3 b) (Schwintzer, 1983;Waughman and Bellamy, 1980) and desert BNF (Negev Desert Israel, Figure 3 b) (Russow et al., 2008). These locations show the systemic problem with BNF predicated on NPP and focused on symbiotic BNF. Although the NPP is relatively low, the BNF is high due to the presence of free-living BNF (Russow et al., 2008;Schwintzer, 1983;Waughman and Bellamy, 1980). Free-living BNF is less likely to adhere to the assumption of being related to plant productivity, as by 165 definition it is not directly associated with plants. Symbiotic BNF represents only 0.11 kg N ha -1 yr -1 in the Negev desert measurements, but the biological crusts fix 9 -13 kg N ha -1 yr -1 (Russow et al., 2008). Considering only symbiotic BNF the models are on the correct order of magnitude. Unlike other observation sites, where some discrepancies between models and observations can be partially attributed to differences in land cover, the models are capturing desert as a low productivity environment, with the NPP and ET based models all having very low BNF. The error therefore, is in the assumption that low 170 productivity equates to low BNF. However, it is unclear at what time-scale free-living BNF becomes available to plants, and it remains uncertain to what extent it contributes to future carbon sequestration. Given that free-living BNF makes up 34 -49% of total BNF, this suggests that in terms of BNF that is bioavailable to plants contributing significantly to NPP, the modelled values ought to be lower than the global estimates shown above.
The low latitudes have a similar observed distribution of BNF to the mid latitudes, but the models generally have higher 175 BNF, with 3 stark examples (Figure 3). These are all forest locations (Tierney et al., 2019;Zheng et al., 2016) with either tropical broadleaf or pine species and relatively high productivity environments. We can see from these locations, as well as the tropical forests of S Costa Rica (Sullivan et al., 2014) that the NPP based models are particularly liable to overestimations of BNF in the tropics.

BNF in Future under SSP3-7.0 180
All the models simulate an increase in NPP over the 21 st century in SSP3-7.0 due to the combined effects of rising atmospheric carbon dioxide and climate change. Given the constrained stoichiometric ratios within plants and soils, such an increase in productivity requires additional nitrogen to sustain growth. Work on the structural uncertainty to the carbon cycle caused by BNF in individual models (Meyerholt et al., 2016Wieder et al., 2015) indicates that changes in the representation of BNF and its assumed dependency on NPP, ET, or plant N demand lead to significant variation in carbon 185 storage under elevated atmospheric carbon dioxide within the same model structure. In the context of these results and the large range of present-day BNF simulated between CMIP6 models, it would be a logical corollary if the magnitude of simulated changes in NPP was associated with the magnitude of simulated change in BNF in the SSP3-7.0 scenario.
However, in this ensemble, the increases in BNF are not proportional to those in NPP, suggesting that other model structural differences affecting the carbon cycle response to atmospheric carbon dioxide and climate change superseded the direct 190 impact of the BNF parameterisation (Fig, 4 a and b).
The models with BNF as a function of NPP should have BNF increases approximately commensurate with their increase in NPP (Figure 4 b). This is true for JSBACH (in MPI-ESM and AWI-ESM) and UKESM1, where relative changes in NPP and BNF fall nearly onto the 1:1 line. CMCC, which employs a similar representation, deviates from parity, because in parts of the tropics the simulated BNF is at the saturation-level of NPP and has reached the model prescribed maximum (see Table  195 1). TaiESM1, which uses the same underlying land model as CMCC, shows a closer relationship between NPP and BNF. This is due to the lower tropical NPP in this model leading to the BNF being further from saturation point compared to CMCC. All these models suggest little change (relative to the whole model ensemble) in net N mineralisation or N loss and MIROC + 17.5, in a model ensemble range of -29.3 to +17.5), probably because of a large fraction of vegetation carbon increase in woody biomass. This C:N ratio change effectively decreases the relative increase in demand for nitrogen associated with the increase in NPP, and illustrates that stoichiometry and BNF together affects the magnitude of the nitrogen constraint on terrestrial carbon storage . 210 The two models (CESM2 and NorESM2) using the FUN (Fisher et al., 2010) carbon cost function for BNF have almost identical absolute changes in NPP, BNF, N dep, Nfert, and N loss (Figure 4 e and SI Figure 4), despite their divergent results for present day (section 3.1). They have the largest increase in BNF but proportionally less NPP change than BNF (they are above the 1:1 line in Figure 4 b), which is probably at least partially related to the fact that these models account for the carbon costs of increasing BNF in the calculation of NPP. In effect, the extra supply of nitrogen via BNF in these models is 215 not converting to an increase in NPP as efficiently as in the rest of the model ensemble. Despite their similar land model, CESM2 has the largest N uptake increase of all the models in the ensemble, whereas NorESM2 is the only model projecting a decrease in plant N uptake (346 TgN yr -1 and -51 TgN yr -1 respectively). This difference is likely related to diverging projections in net N mineralisation, which is 788 TgN yr -1 for CESM2 and 267 TgN for NorESM2 (ensemble range 75.8 TgN yr -1 -385.4 TgN yr -1 ). In contrast, we see only 0.9 PgC yr -1 difference in increase of NPP between these two models, in 220 a model ensemble range of 1.3 -22.6 PgC yr -1 (Figure 4 a). From this we can see that in the underlying model, CLM5, nitrogen limitation plays only a small role in determining NPP and future terrestrial carbon sequestration. The large increase in N uptake in CESM2 compared to NorESM2 is not leading to proportional differences in NPP. Figure   4 e). Like CESM2 and NorESM2, ACCESS has proportionally less NPP increase for the amount of BNF increase, though 225 the absolute levels are much lower for both (Figure 4 b). It is possible that this is due to ACCESS including the phosphorus cycle in addition to the nitrogen cycle, and therefore, increases in NPP are not only constrained by the magnitude and increase in BNF, but also phosphorus availability. ACCESS is the only model where the C:N ratio change is not approximately proportional compared to the rest of the ensemble to the change in NPP (Figure 4 e). Similarly, it has the largest increase in N loss. 230

Discussion
The historical simulations compared to data for BNF reveal a mixed message, with ensemble members generally performing well in high latitudes and at the global total scale, but poorly in the mid-latitudes and tropics. For the simulations under the future scenario SSP3-7.0, the impact of the projected change in BNF over time under elevated atmospheric carbon dioxide and other environmental changes is more difficult to assess because of covarying drivers and multiple interactions of 235 biospheric processes (including the water cycle and vegetation dynamics) that confound the imprint of changing BNF on the terrestrial carbon cycle response.
Limitations in methodology of both observations and model simulations partly account for the lack of agreement between them. For models, the simulations are not specific to the site, but rather taken from the closest grid-cell of a global simulation. Were site-level simulations with observed driving climate data available and the correct land cover (particularly 240 vegetation) prescribed, it is possible models would perform better. For the data, the comparison is made with simple upscaled measurements grouped by biome, which is vulnerable to skew in the underlying data (Davies- Barnard and Friedlingstein, 2020). The underlying BNF data for the historical comparison also has substantial limitations. For instance, the most commonly used method of measuring BNF, acetylene reduction assay, requires calibration to avoid variation of up to two orders of magnitude, that ~70% of studies fail to do (Soper et al., 2021). The literature is also biased away from null 245 results, making an accurate understanding of the processes underlying BNF more difficult. Thus, the problems with model representation of BNF are symptomatic of wider uncertainties in BNF observations and upscaling from single soil cores or plants to ecosystem level.
The challenge for progressing BNF modelling is what would be a suitable replacement for the functions currently used.
Symbiotic fixation is around two thirds of total BNF (Davies- Barnard and Friedlingstein, 2020) and is the focus of the more 250 process based models of BNF (as used in ACCESS, CESM2, and NorESM2). However, work with herbaceous legumes suggests that fixers may have little variation in whole plant biomass whether the nitrogen is fixed or provided as fertiliser, such that the carbon cost of acquiring nitrogen symbiotically may be lower than previously thought (Wolf et al., 2017).
Therefore the process based attempts to establish an empirical relationship between BNF and climate or soil properties at macro scale have not indicated any robust relationship, which is probably primarily the consequence of lack of data 255 availability, and biome based upscaled values have low data support and high uncertainties (Davies- Barnard and Friedlingstein, 2020). Abundance of fixers is an important parameter in the CESM2 and NorESM2 models and has a large impact on total fixation and response to fertilisation (Fisher et al., 2018), but in observations it is poorly constrained (Davies-Barnard and Friedlingstein, 2020) and not well correlated with total fixation rates, to the point of being anti-correlated (Taylor et al., 2019). 260 Free-living BNF makes up around a third of all BNF, is comprised of a heterogeneous set of organisms (Reed et al., 2011), making a single process based model challenging. Thus, the three CMIP6 models that account separately for free-living BNF use static biome level upscaling (ACCESS) or an empirical relationship with evapotranspiration (NorESM2 and CESM2) (see Table 1). The presence of separately defined processes for symbiotic and free-living BNF is not in itself a sufficient condition for a model with better representation of BNF, as shown by the contrasting performances of CESM2, NorESM2, 265 and ACCESS.
In the future scenarios, the multiple sources of uncertainty as to how and to what extent BNF will change make any definitive statements about the capacity of models to capture BNF changes difficult. While increased atmospheric carbon dioxide tends to increase BNF (Liang et al., 2016), nitrogen addition in the form of deposition or fertilisation tends to supress BNF (Zheng et al., 2019), and effects from land use change (Zheng et al., 2020), increased temperature, reduced precipitation 270 and other climate change as well as the potential effects of climate induced land cover change that may alter the composition and location of biomes. It is challenging to predict which of these factors will predominate over the coming century.
Regardless of the change in BNF in future, it is revealing that while single parameter perturbation experiments suggest BNF significantly affects terrestrial carbon storage (Meyerholt et al., 2016;Wieder et al., 2015) when in a dynamic system the effects of BNF are subsumed by structural differences in the nitrogen and carbon models, as well as the larger effects of 275 increasing carbon dioxide. In terms of confidence in model results, the process-based models have clear advantages.
However, that increased complexity does not necessarily translate into increased fidelity in the representation of model BNF.
This could be due to issues with the process-based representation of BNF, systematic problems with the model representation of the wider nitrogen cycle which BNF previously compensated for, or inaccuracies in the observational upscaled data. 280

Conclusions
BNF is an important part of the nitrogen cycle, and previous work has shown how nitrogen availability (Zaehle et al., 2014) and BNF in particular impacts terrestrial carbon storage (Meyerholt et al., 2016;Wieder et al., 2015). Here we have shown that although there are considerable shortcomings in the representation of BNF in CMIP6 models, BNF is not a dominant source of uncertainty in future carbon uptake when considered in the light of the diversity of process representations 285 currently encoded in Earth System Models. While some models have a strong relationship between NPP and BNF, others utilise changes in equally uncertain parts of the nitrogen cycle to enable the increases in NPP. Even models with structural similarity can have almost identical NPP with different levels of BNF and other nitrogen variables. Therefore, the weaknesses of BNF representation in most models need to be seen in the context of the entire nitrogen cycle and the need for simultaneous improvement in process understanding and representation. 290 The two models performing best are CESM2 and CMCC-CM2, as judged by being within the global estimates by Vitousek et al., (2013) and Davies-Barnard and Friedlingstein, (2020) as well as within the uncertainty for all three latitudinal bands (low, mid and high latitudes) for the Davies-Barnard and Friedlingstein, (2020) observation-based estimates. These models have little in common in terms of BNF model processes, encompassing the simplest and most complex modelling processes discussed above. That models using a weak physical basis for BNF can simulate BNF with similar skill to more sophisticated 295 models reveals the need for greater process understanding of BNF to enable more reliability in models. the observationally-constrained ranges by Davies-Barnard and Friedlingstein, (2020) and (Vitousek et al., 2013). The four panels represent latitude bands and global. High latitudes: more than 60°N or 60°S. Mid latitudes: less than 60°N and more than 30°N less than 60°S and more than 30°S. Low latitudes: less than 30°N and less than 30°S.  to the latitude and longitude. Top panel: more than 60°N or 60°S. Middle panel: less than 60°N and more than 30°N less than 60°S and more than 30°S. Lower panel: observations less than 30°N and less than 30°S. The black lines represent single 525 values, or the confidence range as reported by the paper the observational data comes from (see SI Table 1). Figure 4: a and b show the change in BNF and NPP between the first and last decade of SSP3-7.0. c, d, and e show the normalised changes in nitrogen components in the models. Each dot represents the normalised change in [the variable] during the 21st Century, in SSP3-7.0. Each line represents a model, and each plot is a group of models that deal with BNF similarly (from the top left: carbon cost or mechanistic, ET, and NPP). 535 Table 1. Summary of the model's BNF representations. The theoretical maximum BNF value refers to any limit imposed by the equations in the model, e.g. a saturation point.