Nitrogen cycle impacts on CO2 fertilisation and climate forcing of land carbon stores

Anthropogenic fossil fuel burning increases atmospheric carbon dioxide (CO2) concentration, which is adjusting the climate system. The direct impact of rising CO2 levels and climate feedback alters the terrestrial carbon stores. Land stores are presently increasing, offsetting a substantial fraction of CO2 emissions. Less understood is how this human-induced carbon cycle perturbation interacts with other terrestrial biogeochemical cycles. These connections require quantification, as they may eventually suppress land fertilisation, and so fewer emissions are allowed to follow any prescribed future global warming pathway. Using the new Joint UK Land Environment Simulator-CN large-scale land model, which contributed to Coupled Model Intercomparison Project Phase 6 as the land component of the UK Earth System Model v1 climate model, we focus on how the introduction of the simulated terrestrial nitrogen (N) cycle modulates the expected evolution of vegetation and soil carbon pools. We find that the N-cycle suppresses, by approximately one-third, any future gains by the global soil pool when compared to calculations without that cycle. There is also a decrease in the vegetation carbon gain, although this is much smaller. Factorial simulations illustrate that N suppression tracks direct CO2 rise rather than climate change. The finding that this CO2-related effect predominantly influences soil carbon rather than vegetation carbon, we explain by different balances between changing carbon uptake levels and residence times. Finally, we discuss how this new generation of land models may gain further from emerging point knowledge held by the detailed ecological modelling community.


Introduction
Fossil fuel burning and land-use change are increasing atmospheric carbon dioxide (CO 2 ) and concentrations of other greenhouse gases (GHGs), which is driving climate change (IPCC 2013(IPCC , 2021. Both CO 2 increases and climate change forces change to the land surface stores of carbon. Presently, when averaged globally, the land is accumulating carbon, predominantly by direct CO 2 -fertilisation that raises photosynthesis and increases the carbon held in vegetation, C v (GtC). This increase also enhances the levels of litterfall, which increase soil carbon C S (GtC). The indirect effect of CO 2 via climate change can either enhance or suppress the fertilisation effect, depending on the location. Overall, climate change suppresses carbon stores and currently partially counteracts fertilisation (Friedlingstein et al 2014, Arora et al 2020). Therefore, both effects combine to form a net drawdown to land, and the resultant contemporary greening (Zhu et al 2016) offsets ∼25% of CO 2 emissions. It is essential to understand the changes in this value. Any decrease in CO 2 land drawdown implies fewer emissions to track a prescribed global temperature profile or more global warming for any specified emissions trajectory.
However, interactions with other land biogeochemical cycles are often overlooked, which may change as the vegetation size increases. In particular, it has long been suggested that the nitrogen (N) cycle eventually constrains vegetation growth (Hungate et al 2003, Ciais et al 2013, Thomas et al 2015, Zaehle et al 2015, Tharammal et al 2019 as the carbon uptake driven by increasing CO 2 levels outstrips the availability of nutrients for plants to assimilate new carbon. Implementing such a downregulation representation into the land components of Earth system models (ESMs) is occurring (Davies-Barnard et al 2020). In the recent collective modelling exercise, the Coupled Model Intercomparison Project Phase 6 (CMIP6), about half of the ESMs include a terrestrial N-cycle (Arora et al 2020). As expected, the nitrogen cycle has a major effect, across ESMs, of lowering the projected increases in land carbon stores in response to climate and CO 2 rise (Jones and Friedlingstein 2020). We advance the understanding of the results from CMIP6 by examining the behaviour of one of the contributing models, the Joint UK Land Environment Simulator (JULES). Although the earlier Coupled Climate Carbon Cycle Model Intercomparison Project results (Arora et al 2020) contained models with and without a nitrogen cycle representation, inter-model differences make it impossible to conclude with certainty the impact, on its own, of the role of the nitrogen cycle. By utilising a land-surface model offline, able to be run with the nitrogen cycle both included and excluded, we can more robustly identify the impact of nutrient cycling on future projected carbon sinks.
The JULES land surface model (Best et al 2011, Clark et al 2011 was recently enhanced to include an interactive nitrogen cycle (Wiltshire et al 2021), a configuration we call JULES-CN. This version of JULES can be run stand-alone ('offline') or coupled inside the next generation of ESMs such as UKESM1 (Sellar et al 2019). JULES can also operate in the Integrated Model Of Global Effects of climatic aNomalies (IMOGEN) framework (Huntingford et al 2010, Zelazowski et al 2018), based on climate change 'pattern scaling' (Huntingford and Cox 2000). IMOGEN approximates climate change as linear in radiative forcing, and for any ESM it emulates, the constant of linearity depends on the month, near-surface meteorological quantity, and location. A limitation of 'CMIPs' is that it is not possible to explore the structural uncertainty of coupling different land-surface models to alternative atmospheric circulation models. The CMIP ensembles have just one configuration of each land model for each climate model. Here we use IMOGEN to enable us to explore the response of JULES to the full range of CMIP-simulated changes in patterns of climate. Specifically, the IMOGEN patterns of linearity are currently fitted against 34 ESMs in the CMIP5 database (Taylor et al 2012), followed by mapping onto a common spatial grid.
Although IMOGEN simplifies some features of climate response, it acts as a powerful intermediate simulation framework, enabling an early understanding of the implications of introducing new land process representation. The high computational speed of IMOGEN also allows factorial land simulations, which may not be possible with a fully coupled ESM because of their high operational demands. By emulating many ESMs, IMOGEN facilitates the comparison of uncertainty in land surface response to rising atmospheric GHGs against uncertainty in the climate response that is also due to changing GHGs. ESMs differ substantially in expected global warming levels and have many regional differences for identical future scenarios of atmospheric GHG concentrations (e.g. Palmer and Stevens 2019).
We analyse the implications of introducing new nitrogen-based components into the JULES land model, operating in the IMOGEN global climate simulation framework. Our focus is on determining the overall implications of the N-cycle for global land carbon storage as climate changes. To isolate nitrogenbased effects, for comparison, we repeat all calculations with a version of JULES that is identical except without N-cycle components. We present estimates of change for a future pathway in atmospheric GHG concentrations related to high emissions and uncertainty bounds that capture inter-ESM differences, as possible with the IMOGEN system. We then discuss possible further land N-cycle model refinement as new point-process knowledge becomes available.

Methods and simulation structure
A summary of the new nitrogen process description implemented in the JULES-CN land model (Wiltshire et al 2021) is as follows. The JULES model version just before N-cycle implementation (henceforth referred to as 'JULES-C') assumes that prescribed invariant amounts of available nitrogen are sufficient to meet the requirements for vegetation growth and the turnover of soil organic matter. The JULES-CN model improves on this likely flawed assumption and instead provides a fully interactive N-cycle to simulate the influence of varying nitrogen availability on terrestrial carbon (C) stores. JULES-CN focuses on key components of the N-cycle, its coupling with the main C parts already modelled in the JULES-C framework, yet maintains sufficient simplicity and parameter sparseness to enable an eventual full coupling within the UKESM. Here, inorganic N available for biotic use supplies ecosystems through biological fixation of atmospheric N 2 and reactive N deposition. The fixation rate is calculated as a linear proportion of the net primary production (NPP) before N limitation (i.e. NPP potential) and is also dependent on prescribed global, spatially explicit N deposition rates. Soil inorganic N can be taken up by plants, be immobilised by soil microbes, or lost from the terrestrial Plants must acquire inorganic N to allocate C for growth, although the quantity taken up varies by plant functional type (PFT). If the amount of soil inorganic N is less than the amount required for plants to achieve their potential NPP, then the plant releases any excess C it cannot assimilate due to such N limitation, as respiration.
The release of N to soils is by litterfall, and the amount varies by PFT. There are two soil litter pools, decomposable and resistant, each with different turnover rates (yr −1 ). When N is not limiting, the potential mineralisation rate of the organic N in the litter pools is determined by the temperature and moisture near the soil surface, and by recalcitrant litter quality which depends on the fractional cover of each PFT. When the amount of inorganic N is less than the amount required for litter decomposition, the modelled turnover rates of the litter pools are reduced in proportion to the ratio of nitrogen supply relative to demand. This reduction, in turn, lowers the rate at which organic N is mineralised and reenters the soil inorganic pool, thus further constraining ecosystem productivity.
Emissions lose a constant proportion of mineralised N as gas to the atmosphere. Leaching causes a further loss of nitrogen and is a function of the net flux of moisture through the soil, the concentration of inorganic N, and a constant effective solubility of N. A portion of the soil inorganic N pool is also lost to microbial immobilisation. As with mineralisation, the rate of immobilisation depends on the turnover rate of the litter pools and the C:N ratio of the soil pool.
We first operate the JULES-C and JULES-CN models globally and in the IMOGEN climate projection structure. We drive these calculations by an identical prescription of atmospheric CO 2 concentrations that follow their known historical changes and broadly track the future RCP8.5 representative concentration pathway scenario (Meinshausen et al 2011) to the end of the 21st century. Adding to this CO 2 forcing of IMOGEN is an RCP8.5-based pathway in non-CO 2 GHGs and aerosols, expressed as a combined radiative forcing. We use a common CO 2 concentration scenario, rather than one of identical CO 2 emissions, to allow the isolation of direct N-cycle effects independent of their feedbacks via the carbon cycle. Hence IMOGEN does not allow land-and ocean-atmosphere fluxes to interact with atmospheric CO 2 in this configuration. Additional to these two 'all forcings' calculations, we undertake sets of factorial simulations for both model versions. The first set operates JULES-C and JULES-CN with atmospheric CO 2 (and N deposition changes; Lamarque et al 2013) driving the land surface model ('CO 2 only'). These factorial runs, with no climate change and only a physiological response, are of particular interest. Many (e.g. Wieder et al 2015), cite limits to nitrogen availability, named downregulation, as likely to become a substantial restriction on land vegetation fertilisation, growth, and carbon accumulation as CO 2 rises. We also perform simulations with the atmospheric CO 2 concentration prescribed to the JULES model as fixed at the preindustrial level. However, this second set of factorial runs is forced by climate change associated with the historical and RCP8.5 atmospheric GHG concentration scenarios (Meinshausen et al 2011) and as well as changing N deposition. These calculations capture the climatic but not the physiological and biogeochemical responses to rising CO 2 levels and other GHG concentration changes ('climate only'). For JULES-CN, we perform a further factorial simulation with all forcings, except for atmospheric nitrogen deposition levels fixed at pre-industrial levels. An atmospheric gas trace model provides historical and projected deposition levels. Specifically, where N deposition is time-varying, this is adopted from the Atmospheric Chemistry and Climate Model Intercomparison Project multi-model dataset, interpolated to annual fields (Lamarque et al 2013). In all simulations we 'spin up' the combined JULES and IMO-GEN system until all components, particularly soil carbon pools, are in equilibrium with pre-industrial (taken as year 1850) climate and CO 2 levels.

Results
Geographical changes in total land carbon stores, that is, vegetation with soil, are shown in figure 2. All calculations are averages taken across the 34 ESMs emulated in the IMOGEN framework. The left column shows the changes in the JULES-CN simulation by the year 2099 compared to pre-industrial estimates, as based on the atmospheric CO 2 rise of the RCP8.5 scenario. The right column shows the differences in the JULES-CN-based calculations of changes minus the JULES-C-based calculation of changes, where the latter model has only the carbon cycle.
The JULES-CN model projects an increase in land carbon storage at most locations (figure 2(a)). However, the nitrogen cycle suppresses such increases in many places, especially in northern latitudes (figure 2(b)). As expected, most carbon storage gains are due to the direct impact of CO 2 fertilisation on plant physiology (figure 2(c)), and it is this effect that the nitrogen cycle dampens (figure 2(d)). Climate change is detrimental to land carbon stores for large parts of the world (figure 2(e)), although the magnitude of change is generally smaller than that associated with CO 2 fertilisation forcing. Nitrogen forcing has relatively little effect on climate drivers (figure 2(f)). In figure 2(g), an additional factorial simulation operates with all forcings, except for N deposition held at pre-industrial levels. This factorial simulation has values similar to those shown in figure 2(a); therefore, historical and future changes in N deposition have comparatively little impact on carbon store changes. For comparison, we repeat panels figures 2(a)-(f) as figures S1(a)-(f) (available online at stacks.iop.org/ERL/17/044072/mmedia), but instead for a CO 2 concentration pathway that broadly follows the RCP2.6 scenario. There are strong similarities between the geographical features of figures 2 and SI-1, but with the magnitude of changes more minor in the latter, as expected for a simulated future representative of substantial constraints on fossil fuel emissions. The total difference, by year 2100, of the N-suppression of land carbon stores (compared to calculations without the N-cycle) is 117.4 GtC for our RCP8.5 calculations (title, figure 2(b)), while this value is 48.9 GtC for RCP2.6 (title, figure SI-1(b)). To follow the same concentration trajectory in atmospheric CO 2 , these values for each scenario translate to the amount to which compatible emissions are lower due to N-cycle effects.
Projected time-evolving changes in total global terrestrial carbon are shown in the top row of figure 3. For simplicity in isolating the natural N-cycle effect, we perform our simulations without the direct role of human land-use change affecting land carbon stores (although that does imply we cannot compare directly to either CMIP5 or CMIP6 simulations of vegetation and soil carbon stores). Individual changes to vegetation and soil stores are the second and third rows, respectively, as shown in figure 3. There is further disaggregation, by column, to all forcings (first column), CO 2 physiological forcing only (second column), and climate forcing only (third column). The spread of the curves corresponds to the impact of uncertainty in how climate will change, as expressed by the differences between the ESMs that IMOGEN emulates. The top row further supports the findings presented in figure 2 that terrestrial carbon stores will increase, and the direct CO 2 -forcing of the physiological response is the driver of carbon accumulation. Direct climate forcing diminishes the rate of growth. Figure 3 again illustrates that the terrestrial nitrogen cycle has a substantial effect by lowering the CO 2 fertilisation effect ( figure 3(b)). N-cycle inclusion has almost no impact on the climate only calculations ( figure 3(c)). When comparing the global changes in vegetation and soil carbon individually, the projection is that the most substantial impact of N-cycle suppression will occur in soil carbon stocks (figure 3(g) versus (d)). Figure SI-2 is identical to figure 3, but for a CO 2 rise tracking that of RCP2.6, and shows similar Hence, the second column is simulations with changes in the N-cycle minus changes with no nitrogen representation (sometimes called a 'diff-of-diff ' map). The first row is with all forcings. The second row is simulations with atmospheric CO2 concentration changing but only driving the JULES model (N deposition also varies). These simulations capture CO2-induced physiological changes, but with climate fixed at pre-industrial, i.e. 1850 levels. The third row is with climate change and N deposition changing, but the CO2 values provided directly to JULES held at pre-industrial levels. Panel (g) is the change in JULES-CN projections, with all forcings changing except N deposition, which is held at pre-industrial levels. In every panel, all values are the average across the 34 ESMs emulated with IMOGEN. The area integrated total global carbon is shown in the title of each panel. The future scenario is for CO2 concentrations rising close to those of the representative concentration pathway, RCP8.5. but smaller time-evolving global changes in carbon stores.
The standard deviation of the spread in changes to the terrestrial carbon pools in 2099 (due to the different climate projections across the ESMs; dots, figure 3(a)), is approximately 15% lower for the simulations with the N-cycle because of their smaller changes. Our calculations are based on the prescribed changes in atmospheric CO 2 . However, if we had used a common CO 2 emissions pathway to force our model framework, the lowered ability of simulations with N-cycle to accumulate carbon would place more CO 2 in the atmosphere. Higher atmospheric CO 2 will cause more climate change, which accentuates the differences between ESMs. Hence, any inclusion of that effect by modelling an emissions-driven fully interactive carbon cycle will expand the JULES-CN plume compared to that of JULES-C. However, an additional consideration is that with the nitrogen cycle suppressing the uptake of CO 2 , then for a common   emissions pathway, this will cause higher atmospheric CO 2 levels and enhance fertilisation in the JULES-CN simulations. A simple order-of-magnitude argument is that the difference between JULES-C and JULES-CN by the year 2099, and for any given ESM, is approximately 150 GtC. If this carbon is placed directly in the atmosphere, it corresponds to ∼+70 ppm of CO 2 , although the ocean and land drawdown will substantially lower this value. The RCP8.5-based CO 2 concentration changes are approximately 7 ppm yr −1 in the year 2099. Hence, for a modelled interactive carbon cycle with JULES-CN, the new climate state would occur less than a decade earlier than with JULES-C. From figure 3(a), the JULES-CN curves for 2089 are still below those for JULES-C for the year 2099. Therefore, we estimate that the terrestrial nitrogen cycle will cause more climate change for the same emissions trajectory than if this effect was not present, due to 'drawing down' less atmospheric CO 2 . Alternatively, an interpretation of figure 3(a) is that the terrestrial N-cycle will allow lower available CO 2 emissions to follow any prescribed atmospheric GHG trajectory, compared to ignoring this effect. Extra time-evolving global mean model diagnostics (figure 4) help explain the numerical findings in figure 3. We find that simulated nitrogen effects dampen projected NPP increases compared to calculations without the N-cycle ( figure 4(a)). Additionally, we show the ratio of NPP to gross primary production (GPP) with and without the N-cycle ( figure 4(b)). This ratio illustrates the fractional capability of the land surface to retain photosynthetic uptake for use in above-ground land carbon pools. Towards the end of the 21st century, the ratio decreases earlier and more rapidly with the inclusion of the N cycle. We also calculate, at each location, carbon residence time derived as the timeevolving gridbox mean carbon content divided by gridbox mean NPP. Figure 4(c) shows the global spatial mean value of these vegetation carbon residence times. Although the basis of residence time values in figure 4(c) is a spatial aggregation of their local values, a comparison with the global mean NPP values of figure 4(a) provides a key insight into why the N-cycle has relatively little impact on the temporal variations in vegetation carbon ( figure 3(d)). Specifically, there is a balance between smaller increases in NPP occurring in tandem with smaller decreases in vegetation residence time (year) (i.e. higher turnover), which weakens the impact of the N-cycle changes on the vegetation carbon stores. The relative differences between carbon residence time changes for JULES-CN and JULES-C are much smaller for soil carbon ( figure 4(d)). We repeat figure 4 for simulations following a RCP2.6-based scenario ( figure SI-3).
We split the findings of figure 4 to the individual PFT levels. This disaggregation allows an understanding of JULES-CN versus JULES-C projections for vegetation residence time at the PFT level and differences in projected changes in fractional cover. Following the format of figure 4, we present time-evolving changes of the areally averaged NPP for each PFT in supplementary information figure SI-4. We also show PFT-based changes in NPP/GPP (figure SI-5) and vegetation carbon residence time ( figure SI-6). Each gridbox has a single soil carbon pool common to all PFTs, and so there is no PFTspecific equivalent diagram to figure 4(d). Instead, we show the fractional cover changes ( figure SI-7). We summarise PFT-based changes (table 1)

Discussion and conclusions
JULES-CN represents N limitation by reducing tissue allocation within terrestrial plants once carbon has been assimilated (i.e. reductions in NPP, not GPP). Therefore, when the amount of N required for growth is not available, any excess fixed carbon is respired back to the atmosphere. Incorporating N limitation thus causes a lower carbon use efficiency relative to the C-only model ( figure 4(b)). Over time, the relative residence time of carbon in vegetation decreases less with the N-cycle. Therefore, incorporating the N-cycle results in lower above-ground productivity, but the carbon sequestered in vegetation remains in this pool for longer. Although there is some heterogeneity between PFT responses, the feature of smaller increases in NPP and smaller decreases in residence time holds for many (table 1). This finding agrees with the interpretation of figure 3(d) that these two factors combine and offset. Such offsetting causes vegetation carbon stores to have relatively similar future changes for JULES-CN and JULES-C (figure 3(d)).
Lower increases in NPP with the modelled N-cycle generate lower projected rates of litter production and thus carbon flux to soils. However, the residence time for soil carbon is more similar between the two models ( figure 4(d)). Together, these two soil factors cause much smaller estimated increases in soil carbon when accounting for the N-cycle ( figure 3(g)). Therefore, our simulations project soils rather than vegetation explain most suppressed increases in total land carbon changes when including nitrogen interactions ( figure 3(a), JULES-CN versus JULES-C).
Our finding that globally, vegetation turnover decreases more slowly into the future for simulations with the N-cycle (offsetting the slower increase in NPP due to N-cycle limitation) is an emergent aggregated property requiring detailed investigation. One driver may be vegetation competition, where lower NPP with the N-cycle reduces losses due to self-shading between PFTs, reducing residence time decreases. A second factor could be that in this configuration of JULES-CN, nutrient limitation has a relatively greater impact on the low biomass grassy PFTs. As such, PFTs with smaller vegetation carbon contents are more controlled by N effects (e.g. in the Arctic and comparing grasses with nearby other more woody vegetation). Such differentiation between PFTs may lower the overall effect of nitrogen on vegetation stores when averaged globally. We plan to extract timeseries from JULES simulations of all internal calculations, for every modelled timestep, and at representative locations and PFTs to isolate and understand these effects. From these high temporal timeseries, we will determine the precise balance of equation terms that leads to our discovered emergent global property of low N-cycle influence on total vegetation carbon stores.
JULES-CN adopts much of the current understanding of the N-cycle, guided by the recommendations of Zaehle (2013). We now highlight two aspects of the current JULES-CN model (Wiltshire et al 2021) that could be prioritised for future updates because they can improve the estimates of downregulation without introducing significant model complexity.
First, using a fixed rate of biological nitrogen fixation (BNF) as a proportion of NPP likely overestimates the supply of reactive N under elevated CO 2 . The JULES-CN model predicts BNF values (102 Tg N yr −1 ; Wiltshire et al 2021) that fall within the observed global range of 52-130 Tg N yr −1 (Davies-Barnard and Friedlingstein 2020) but is higher than other recent model-based predictions of 61.5 Tg N yr −1 (Yu and Zhuang 2020). Furthermore, BNF may depend on climate and N deposition (Peng et al 2020). Overall, BNF could be the largest source of uncertainty in ensemble predictions of carbon uptake (Meyerholt et al 2020). Thus, the assumption of a fixed rate of BNF may underestimate the extent of N limitation. Hence, allowing for greater complexity in BNF representation would likely lead to predictions of even greater JULES-C versus JULES-CN differences. Second, the model assumes that PFTs have a fixed stoichiometric relationship between C and N, which constrains the capacity of plants to maintain productivity when N supply declines. A recent metaanalysis showed that changes in plant C:N ratios were responsible for half of the additional C sequestered in forests under CO 2 enrichment (Zou et al 2020). Flexible stoichiometry allows growth and spreading for a lower nutrient cost. Therefore this second assumption may cause the current version of the JULES-CN model to overestimate the extent of N limitation. Incorporating stoichiometric flexibility may yield higher projected levels of future vegetation cover and litter production, albeit at higher C:N ratios, with as yet unknown impacts on soil processes (Camenzind et al 2021). Data to support model improvements include free-air CO 2 enrichment (FACE) experiments (Norby and Zak 2011) that subject plots of land to continuously raised atmospheric CO 2 . Analysis of changes to NPP or carbon stocks reveals the extent of potential CO 2 -induced fertilisation, including the opportunity to determine the role of the land N-cycle (Zak et al 2003). The operation of such experiments for many years or even decades is pivotal, as the N limitation of biomass gains with CO 2 enrichment, e.g. temperate grasslands, takes time to develop (Reich and Hobbie 2013). Many request increasing the geographical extent of FACE experiments (e.g. Jones et al 2014). Existing FACE measurement sites have provided a valuable benchmark of large-scale land models that include N-cycle representation (Zaehle et al 2014). FACE experiments demonstrate the need for land surface models to incorporate mycorrhizal fungal impacts on nitrogen availability and plant growth in response to elevated CO 2 (Terrer et al 2016). There are significant differences in the predicted CO 2 fertilisation effect depending on the form of microbial symbiont (Terrer et al 2019). Plants associated with ectomycorrhizal (ECM) symbionts can increase biomass without nitrogen limitation, whereas growth in species associated with arbuscular mycorrhizae (AM) are found to be N-limited in FACE experiments (Terrer et al 2016). In ECM-associated plants, this additional plant biomass links to lower soil C stocks because soil organic matter decomposes faster and soil C respiration increases, negating the fertilisation effect (Terrer et al 2021). The opposite occurs in AM-associated plants as these are more N-limited; therefore, additional fixed carbon is transferred to roots (as exudates for microbes) and remain in soils for longer. Such longer residence times are because, in the absence of additional N inputs, decomposition rates are lower under enriched CO 2 . Ignoring the different forms of mycorrhizal fungi is expected to result in underestimates of plant aboveground biomass growth and overestimates of soil C storage in typically ECM-associated forests, and the opposite in tropical grasslands (Terrer et al 2021).
Two additional mechanisms should be introduced to the JULES-CN model: geological supplies of reactive nitrogen ( At a scale of hundreds of kilometres or more, detection and attribution (D&A) methods may validate modelled N-cycle effects for the contemporary period. These methods have guided key climate change conclusions e.g. the 4th Intergovernmental Panel on Climate Change (IPCC) report (IPCC 2007) statement that warming is 'unequivocal' and 'likely' caused by GHG rises. D&A analyses are a form of regression that sum factorial simulations and fit them to a gridded dataset. For climate, these factors are the near-surface temperature response to rising atmospheric GHGs, changing aerosols, volcanos, and solar cycles (Hegerl et al 1997, Stott et al 2006, fitted to a temperature dataset, e.g. of the Climate Research Unit (Harris et al 2020). Regression variables are spatiotemporal patterns accounting for different geographical forcing fingerprints (e.g. short-lived aerosols have strong regional impacts near industrial centres). If the bounds on a regression coefficient do not include zero (and ideally include unity), then some observed changes are attributed to that driver. This approach could have applicability to our factorial simulations, to determine if projected spatial variations in the effects of the N cycle are detectable in global datasets of vegetation behaviour, e.g. the Sentinel-2 satellite mission (Berger et al 2012). However, we reiterate our finding that the N-cycle in JULES-CN has a larger impact on soil carbon stocks than vegetation. Beyond regression, machine learning algorithms may help identify any nonlinear fingerprints of changing large-scale environmental factors (Huntingford et al  2019).
We support the on-going focus on modelling a closed global carbon cycle. Early attempts to understand the global atmosphere-land CO 2 drawdown were cautious of land model predictions and instead used the residual between carbon emissions, atmospheric storage, and modelled atmosphere-ocean fluxes (Canadell et al 2007). Newer assessments of carbon cycle evolution use terrestrial ecosystem models and compare them against atmospheric inversions (e.g. Friedlingstein et al 2020). This recent assessment states that discrepancies between forward projections by land models and atmospheric inversion methods may be due to uncertainties in nitrogen-related complexities.
We have operated, globally, a modern land model with an interactive N-cycle (JULES-CN). In addition to carbon impacts, the N-cycle influences physical processes such as runoff (Yang et al 2019). Our main simulations are for the RCP8.5 scenario, corresponding to GHG emissions continuing at relatively high levels. IMOGEN spans uncertainty in climate change by emulating different ESM projections. Considering the N-cycle of one land model, driven by multiple ESMs outputs, contrasts with assessing multiple land models, but each forced by a single (but different) ESM (Davies-Barnard et al 2020). The calculation speed of IMOGEN allows factorial calculations that may be prohibitive with computationally demanding ESMs. Although IMOGEN is an offline model that does not capture local land-atmosphere feedbacks from different terrestrial model configurations, our factorial calculations remain informative. We find that the terrestrial nitrogen cycle substantially reduces future increases in land carbon stores, while changes in N deposition levels have relatively little impact. We estimate that the most significant impact of incorporating N-induced suppression is on soil carbon. Less capability of the land to draw down atmospheric CO 2 implies that lower GHG emissions are needed to fulfil any climate policy targets, such as constraining global warming to two degrees above pre-industrial levels.

Data availability statements
The JULES-CN and IMOGEN models are available upon request from the authors. The ensemble of Earth System Models, against which IMOGEN is calibrated, are held in the CMIP5 database (https://pcmdi.llnl.gov/mips/cmip5/).
The data that support the findings of this study are available upon reasonable request from the authors.