JULES-CN: a coupled terrestrial Carbon-Nitrogen Scheme (JULES vn5.1)

Understanding future changes in the terrestrial carbon cycle is important for reliable projections of climate change and impacts on ecosystems. It is known that nitrogen could limit plants’ response to increased atmospheric carbon dioxide and is therefore important to include in Earth System Models. Here we present the implementation of the terrestrial nitrogen cycle in the JULES land surface model (JULES-CN). Two versions are discussed the one implemented within the UK 5 Earth System Model (UKESM1) which has a bulk soil biogeochemical model and a development version which resolves the soil biogeochemistry with depth. The nitrogen cycle is based on the existing carbon cycle in the model. It represents all the key terrestrial nitrogen processes in an efficient way. Biological fixation and nitrogen deposition are external inputs, and loss occurs via leaching and a bulk gas loss parameterisation. Nutrient limitation reduces carbon-use efficiency (CUE ratio of 10 net to gross primary productivity) and can slow soil decomposition. We show that ecosystem level limitation of net primary productivity by nitrogen is consistent with observational estimates and that simulated carbon and nitrogen pools and fluxes are comparable to the limited available observations. The impact of N limitation is most pronounced in northern mid-latitudes. The introduction of a nitrogen cycle improves the representation of interannual variability of global net ecosystem exchange 15 which was much too pronounced in the carbon cycle only versions of JULES (JULES-C). It also reduces the CUE and alters its response over the twentieth century and limits the CO2-fertilisation effect, such that the simulated current day land carbon sink is reduced by about 0.5 Pg C yr−1. The inclusion of a prognostic land nitrogen scheme marks a step forward in functionality and realism for the JULES and UKESM models. 20 1 https://doi.org/10.5194/gmd-2020-205 Preprint. Discussion started: 24 July 2020 c © Author(s) 2020. CC BY 4.0 License.


Introduction
Terrestrial ecosystems absorb around 25% of anthropogenic carbon emissions (Le Quéré et al., 2018;Friedlingstein et al., 2019), and changes in the future land carbon sink will feedback to climate via the proportion of the emissions remaining in the atmosphere. Under projected climate change, the primary mechanism for increased terrestrial sequestration is an increase in plant productivity and 25 biomass, which relies on sufficient availability of nitrogen within the soil-plant system. Therefore the availability of nitrogen impacts the land carbon sink, both in the present and in a higher atmospheric carbon dioxide (CO 2 ) future.
Nitrogen exists in the terrestrial system in organic and inorganic forms and is continually cycled. 30 In a stable climate the external inputs-biological fixation and nitrogen deposition -are balanced by the external losses-leaching and gas loss. Depending on the nutrient status of the vegetation and soil, changes in the balance of the inputs and outputs of nitrogen can drive adjustments in vegetation biomass and soil organic matter. Internally organic nitrogen is lost from vegetation through the production of litter and disturbance. The litter decomposes into soil organic matter and in turn is 35 mineralised into inorganic nitrogen. Both inorganic and organic nitrogen may become available for plant uptake (Weintraub and Schimel, 2005).
In a changing climate, rising atmospheric CO 2 drives an increase in the terrestrial carbon cycle and Gross Primary Productivity (GPP). This extra demand for nitrogen will limit the potential in-40 crease in future carbon stocks. For example, Zaehle (2013a) suggest that, in some areas, nitrogen could limit future carbon uptake by up to 70%. Nitrogen cycling also tends to reduce the sensitivity of land carbon uptake to temperature. Warmer conditions lead to increased plant respiration and soil respiration, which tends to reduce the land carbon sink. However, the increased soil respiration also leads to accelerated nitrogen mineralisation and increased nitrogen availability to plants, which may 45 provide a counteracting increase in GPP. This latter effect is absent from models that do not include a nitrogen cycle, As a result of neglecting these important effects, land-surface models without an interactive nitrogen cycle tend to overestimate both CO 2 and temperature effects on the land carbon sink (Wenzel et al., 2016;Cox et al., 2013). An increasing number of land surface and climate models now include constraints on the land carbon sink caused by nitrogen limitation (Zaehle et al.,50 2014; Wania et al., 2012;Smith et al., 2014). Recent simulations have generated a range of estimates for the sensitivity of the C cycle to N availability (Meyerholt et al., 2020a;Arora et al., 2019). For example, Meyerholt et al. (2020a) use a perturbed model ensemble and show that the projected future increase in land carbon store caused by CO 2 fertilisation is reduced by between 9 and 39 % due to nitrogen limitation and the projected losses in terrestrial carbon caused 55 by climate change are between a reduction of 39 % and a slight increase of 1 %. 2 https://doi.org/10.5194/gmd-2020-205 Preprint. Discussion started: 24 July 2020 c Author(s) 2020. CC BY 4.0 License. The purpose of this paper is to describe and evaluate the implementation of a coupled carbon and nitrogen cycle within the Joint UK Land-Environment Simulator (Best et al., 2011;Clark et al., 2011) (JULES at vn5.1 -http://jules-lsm.github.io/vn5.1/release_notes/JULES5.1.html). JULES is 60 the land surface component of the later generation of Hadley Centre climate models including the UK Earth System Model (UKESM) (Sellar et al., 2019). The addition of the nitrogen component described here was carried out alongside other developments such as improved plant physiology and extended plant functional types (Harper et al., 2018), an enhanced representation of surface exchange and hydrology (Wiltshire et al., 2020) and a new managed land module (Robertson and Liddicoat,65 in prep.). These separate components have been combined to make the land surface component of UKESM and were used for the most recent Global Carbon Budget annual assessment .
The philosophy behind the developments described here is to produce a parsimonious model that 70 captures the large-scale role of nitrogen limitation on carbon use efficiency (CUE -ratio of net to gross primary productivity) and net ecosystem productivity (NEP). This is achieved by extending the implicit representation of nitrogen in the existing dynamic vegetation and plant physiology modules to enable a more comprehensive nitrogen cycle within the land surface. Nutrient limitation operates through two mechanisms; the available carbon for growth and spreading is reduced, and the decom-75 position of litter carbon into the soil carbon is slowed. This is achieved by explicitly representing the demand for nitrogen within the vegetation and soil modules and then reducing plant net carbon gain to match available nutrients. In the soil module an additional decomposition rate modifier is introduced that slows decomposition to match available nutrients. The current structure of the TRIFFID dynamic vegetation model (Cox, 2001), in particular the fixed allometry and carbon allocation, is 80 largely unchanged. As the aim of this scheme is to capture the impact on terrestrial carbon stores, loss terms are aggregated and not speciated. The model's reduction of vegetation growth and spreading due to nitrogen limitation will have only a minor impact on the GPP and autotrophic respiration.
Therefore the emergent impact of the nitrogen scheme will be to reduce NPP and hence the carbon use efficiency of the vegetation. The excess carbon which cannot be used for growth goes to 85 non structural carbohydrates, root exudates and biogenic volatile organic compounds (Collalti and Prentice, 2019).
Two nitrogen model configurations are described here-JULES-CN and JULES-CN layered -both of which are directly derived from the JULES-C configuration. JULES-C is the land configuration of the HadGEM2-ES (Collins et al., 2011) Earth System Model used in CMIP5 (Taylor et al., 2012), 90 and is also used in the Global Carbon Budget annual assessments (Le Quéré et al., 2018) coupled carbon-nitrogen model based on JULES-C. The soil biogeochemistry is represented by a single level in JULES-CN whereas it varies as a function of depth in JULES-CN layered . This paper describes the additional model structure required for the two configurations; and assesses the simulated stocks and fluxes and their changes over the 20 th century. 95 2 JULES model description JULES is the land surface component of the new UK community Earth System model, UKESM1 (Sellar et al., 2019). JULES can also be run offline forced by observed meteorology globally, regionally or at a single location. A full description of JULES is provided by Best et al. (2011) and Clark et al. (2011). In particular, JULES represents the surface energy balance, a dynamic snowpack  (Cox, 2001). In this version of TRIFFID, five plant functional types (PFTs) are included: broadleaf tree, needleleaf tree, C 3 grasses, C 4 grasses and shrubs.

105
The soil carbon model in JULES-C is based on the RothC model (Clark et al., 2011). Recently, Burke et al. (2017); Chadburn et al. (2015) added a representation of permafrost soil processes to JULES, including a representation of the vertical distribution of soil carbon which we build upon here.

Model developments
What follows is a description of the extension of the carbon cycle used by JULES-C in HadGEM2-ES (Collins et al., 2011) and Global Carbon Budget annual assessment in 2018 (Le Quéré et al., 2018) to include an interactive nitrogen cycle. As standard, JULES-C includes an implicit representation of nitrogen which has been extended to be fully interactive. For clarity we include a full

Vegetation carbon and nitrogen
At the core of the vegetation model is the TRIFFID Dynamic vegetation model (Cox, 2001). TRIF-FID represents the vegetation cover at each location in terms of the fractional area covered, and the leaf area index and canopy height of each PFT. The mean canopy height is converted via allometric 140 equations into a maximum or balanced leaf area index (L b ) for each PFT. The vegetation carbon density per PFT (C v ) can be separated into leaf (L c ), root (R c ) and total stem (W c ) pools, each of which is related allometrically to L b : Where σ l , a wl and b wl are PFT dependent allometric parameters ( where η sl relates respiring stem to leaf carbon (Table 1). We can combine equations 4 and 5 to relate (L b ) to canopy height (h) and these two variables can be used interchangeably to describe the state of the vegetation.

160
The root and total stem nitrogen pools are defined using stoichiometric relationships as a function of the carbon pools. These stoichiometric functions already exist in the model and are used in the calculation of plant maintenance respiration. We extend their use to explicitly define nitrogen pools as part of the new scheme.
where µ rl and µ sl are stoichiometric parameters linking the top leaf nitrogen concentration (n l0 ) to the total stem and root nitrogen pools (W n and R n , respectively). The leaf nitrogen pool (L n ) has an additional dependency on phenological state (3.1.2) and assumed distribution of nitrogen in the 170 canopy (3.1.3). Following Equation 1 the total vegetation nitrogen store per PFT is given by: The C:N ratio of the root and stem pools are constant in time and leaf pool C:N ratio only varies with phenological state. However, the relative proportions of each pool vary with total biomass resulting in the whole plant C:N ratio increasing with total vegetation carbon (Fig. 2). This is due to 175 the relatively greater proportion of stem carbon at higher biomass. Therefore woody vegetation has the highest C:N ratios due to the greater proportion of stem wood in comparison to grasses. The total vegetation nitrogen increases with canopy height (Fig. 3).

Biological Nitrogen Fixation
Biological nitrogen fixation (BNF) is assumed to be the largest supplier of nitrogen to the terrestrial 180 ecosystem. Following Cleveland et al. (1999), the nitrogen fixation is determined as a proportion of the net primary production before nitrogen limitation (N P P pot ).
The rate of fixation (ζ) is set such that global present day net primary productivity of approximately 60 Pg C yr −1 results in approximately 100 Tg N yr −1 fixation (0.0016 kg N kg C −1 ), within 185 the range most recent global estimate of BNF (Davies- . In JULES-CN this fixation is directly transferred into the inorganic soil nitrogen pool and becomes available as

195
The leaf carbon pool (L c , Equation 2) varies allometrically with the vegetation carbon state on both short (seasonal) and long (centennial) timescales but not with changes in phenological state. Implicit within TRIFFID is a labile leaf carbon pool that acts as a reserve of carbon during spring and a store during fall. L c is therefore the sum of a labile pool from which carbon can be mobilised during leaf out plus an allocated pool representing the actual leaf area index. The labile pool is zero at full 200 leaf out and at the allometrically defined maximum during the no leaf period. This distinction is inconsequential in the carbon only mode but is more critical when considering nitrogen interactions as the implication is that at all times the plant has enough nitrogen in reserve to maintain full leaf. 7 https://doi.org/10.5194/gmd-2020-205 Preprint. Discussion started: 24 July 2020 c Author(s) 2020. CC BY 4.0 License.
We therefore a new parameterisation of retranslocation and labile nitrogen that is dependent on the phenological state. Leaf phenology is controlled by a second state variable (p) which relates the leaf 205 area index at any moment in time to the balanced leaf area index (L b ).
where p is a scalar between 0 and 1 that describes the phenological state of the system (Clark et al., 2011). For evergreen plants p is a constant of 1. The two state variables L b and p combine to define the vegetation state.

210
Using the phenological state we extend the equivalent approach to leaf carbon to include the role of retranslocation of nitrogen from the leaves during leaf fall. The leaf nitrogen pool is the sum of allocated and labile components with additional dependencies on the distribution of nitrogen in the canopy and phenological state. This means that the stochiometry of the allocated and labile components are different. During leaf-off the labile component is the equivalent of the retranslocated 215 leaf nitrogen plus an additional store of nitrogen in preparation for the following bud burst. Higher retranslocation implies a larger labile nitrogen store. L n therefore becomes: where λ l is the leaf nitrogen retranslocation coefficient and n lc is the mean canopy nitrogen content (Eq 12). In this configuration λ l is set to 0.5 for all PFTs (Zaehle and Friend, 2010). The labile 220 pool is formulated so that around half of the nitrogen required for full leaf-out is taken from retranslocation with a further quarter acquired during the dormant phase, the rest is acquired during the active period.

Canopy nitrogen
JULES includes a process-based scaling-up of leaf level photosynthesis to the the canopy level.

225
There are two options for the canopy scaling up including the 'big-leaf' and a 'multi-layer' approach.
In the JULES-CN and JULES-CN layered configurations, to be consistent with JULES-C model, we assume a multi-level canopy with leaf nitrogen decreasing exponentially through the canopy (CanRadMod 5). The mean canopy nitrogen content is described by (Mercado et al., 2007): where k n is a constant representing the profile of nitrogen and z represents the fraction of canopy above the layer. Based on observed nitrogen profiles in the Amazon basin (Carswell et al., 2000), a value of 0.78 for k n was found (Mercado et al., 2007). Equation 12 is independent of leaf area and therefore equates to a constant of proportionality relating PFT-specific top leaf nitrogen (n l0 - Table   8 https://doi.org/10.5194/gmd-2020-205 Preprint. Discussion started: 24 July 2020 c Author(s) 2020. CC BY 4.0 License. 1) to the mean canopy nitrogen concentration. Canopy Leaf C:N ratios are resultingly44% higher 235 than top leaf ratios.

Allocation
Net Primary Productivity in JULES-C is simply the difference between GPP and autotrophic respiration. In JULES-CN the GPP, N P P pot and autotrophic respiration are calculated at the model timestep (1 hour in JULES-C) prior to any N limitation. These fluxes are then aggregated to the 240 timestep for running TRIFFID (once every 10 days in the JULES-C configuration) so that allocation of carbon can take place. N P P pot supplied to TRIFFID represents the potential amount of carbon that can be allocated to growth and spreading (spreading is the increase in PFT fractional coverage).
In order for the NPP to achieve its potential it needs to be able to uptake sufficient inorganic nitrogen.
If not enough inorganic nitrogen is available, the system is nitrogen limited and an additional term is 245 required in the carbon balance representing excess carbon which cannot be assimilated into the plant due to lack of available nitrogen (Ψ c ). A positive Ψ c results in a reduction of carbon use efficiency.
The PFT carbon balance is therefore: where Π c is the net primary productivity per unit area of PFT (prior to nutrient limitation) and Λ lc 250 is the PFT specific litterfall rate (Section 3.1.5). Any excess carbon (Ψ c ) is considered an additional plant respiration term and at the end of the TRIFFID timestep is used to reduce N P P pot to its actual value. λ is the coefficient for partitioning the NPP between growth and spreadingλ is utilised in increasing the fractional coverage of the vegetation and (1 − λ) increases the carbon content of the existing vegetated area. λ is a function of the vegetation carbon which itself is a function of the Should L b fall below L min then the carbon available for spreading is decreased and L b set equal to L min and the carbon pools re-diagnosed. If L b rises above L max then the carbon available for spreading is increased and L b set equal to L max and the carbon pools re-diagnosed.
where Φ is the plant nitrogen uptake (see below) and Λ ln is the retranslocation of nitrogen from leaves and roots into the plant labile pool (Section 3.1.5 below). The nitrogen available for spreading is a fraction λ of the total available nitrogen and (1 − λ) is available for growth.

265
In JULES-CN, on a PFT basis, the nitrogen available for plant uptake is the inorganic soil N pool, N in , split equitably between the PFTs assuming there is no differential ability between PFTs to acquire nitrogen. On a grid cell basis, since the competition for nitrogen depends on the change in carbon over the timestep, larger PFTs have an advantage. The nitrogen uptake in JULES-CN layered is more complicated and is discussed in Section 3.3.1.

270
The nitrogen available for growth is the total available nitrogen multiplied through by (1 − λ).
Equations 13 and 15 are then solved by bisection such that the nitrogen uptake for growth (Φ g ) is less than or equal to the available nitrogen and the balanced LAI (L b ) remains within the specified upper and lower limits (L min ,L max ). Following the solution of dNv dt the carbon store and balanced LAI (L b ) is updated and the leaf, root and wood pools derived following the allometric equations 275 (Equations 2-4).
The remaining proportion (λ) of NPP and nitrogen after growth is is allocated to spreading. The nitrogen demand for spreading is equal to the carbon allocated to spreading scaled by the whole plant stoichiometry: where Ψ s is the excess carbon term from spreading and Nv Cv defines the whole plant C:N ratio. The equation is solved such that (Φ s + Φ g ) ≤ N avail and Ψ s is minimised.
Total excess carbon per PFT is therefore the combination of those from growth and spreading 285 Total nitrogen uptake per PFT is therefore the combination of those from growth and spreading The total gridbox nitrogen uptake and excess carbon are therefore the PFT fraction (v i ) weighted total: This is the excess carbon (Ψ) that is considered an additional plant respiration term and at the end of the TRIFFID timestep and is used to reduce N P P pot to its actual value.
Carbon and nitrogen allocated to spreading allow the vegetation to expand onto bare ground.
Where space is limiting the vegetation competes for space with some carbon being turned over as 295 litter. This is handled in the Lotka-Volterra competition routines (see Clark et al. (2011) for full details). Nitrogen only indirectly affects competition through the PFT specific allometric relationships.
The competition code updates the vegetation fractions (v i ).

Litter Production
Litter is produced by the turnover of the leaf, wood and root pools and through the vegetation dy-300 namics due to large-scale disturbance and density-dependent litter production. The PFT specific litter production (Λ lc ) is defined as: where γ r and γ w are parameters and γ l is a temperature dependent turnover rate consistent with the phenological state (Clark et al. (2011). Total litterfall is made-up of the area-weighted sum of the 305 litterfall from each PFT, along with large-scale disturbance rate, and a density dependent component from intra-PFT competition for space.
where c ij are the competition coefficients describing the effect of PFT i on j, γ vi is a large scale disturbance term and v i is the vegetation fraction. The effect of nitrogen limitation on the litter 310 carbon flux is captured in the excess carbon term Ψ i .
The nitrogen equivalent of the PFT specific litter production (Λ ln ) allows for retranslocation of nitrogen from leaves and roots into the labile store.
where λ l and λ r are the retranslocation of leaf and root nitrogen coefficients, set at 0.5 and 0.2 315 (Zaehle and Friend, 2010). Similarly to the total litter carbon, total litter nitrogen is given by:

Soil Biogeochemistry
The soil biogeochemistry in JULES-CN follows the Roth-C soil carbon model (Jenkinson et al., 1990;Jenkinson and Coleman, 1999)   (1 − β R ) is the fraction of soil respiration that is emitted to the atmosphere. β R depends on soil texture (see Clark et al. (2011) for more details). The nitrogen pools follow a similar structure to the 340 carbon pools: For each soil carbon pool, the potential respiration -i.e. the respiration rate when the nitrogen in 355 the system is not limiting -is given by (R i,pot ): where the k i are fixed constants in s −1 (Clark et al., 2011). The functions of temperature (F T (T soil )) and moisture (F s (s)) depend on the temperature (T soil ) and moisture content (s) near the surface.
The function F v (v) depends on the vegetation cover fraction (v). The potential mineralisation of 360 organic N when the system is not N limited (M i,pot ) is related to the potential respiration rates by the C to N ratio of each pool (CN i ): The potential immobilisation of inorganic nitrogen into the organic nitrogen pools (I i,pot ) is related to pool potential respiration (R i,pot ), the respired fraction (β R ) and the C to N ratio of the 365 destination pool in the decomposition chain (CN i ): When nitrogen is limiting, the respiration of the DPM and RPM pools into the soil organic matter pools is additionally limited by the availability of nitrogen: where i is one of RP M or DP M . F N (N ) is the litter decomposition rate modifier and is given by the ratio of the nitrogen available in the soil to the nitrogen required by decomposition. F N is limited to a range of 0.0 to 1.0. When F N is equal to 1, the decomposition, mineralisation and immobilisation take place at the potential rate and the system is not nitrogen limited. Where F N is less than 1, the availability of N limits the decomposition of litter into soil organic matter. F N is 375 given by: where N avail is the soil available inorganic nitrogen in kg N m −2 . D DP M and D RP M are the net demand associated with decomposition of each of the litter pools: If the net mineralisation is positive some of the nitrogen is emitted as gas, according to: where N gas is the gas emission in kg N m −2 s −1 and f gas is the fraction of the nitrogen flux that is emitted as gas. M tot = M DP M + M RP M + M BIO + M HU M and is the total mineralisation flux 390 in kg N m −2 s −1 . Following Thomas et al. (2013), it is assumed that 1% of net mineralisation is emitted as gas (f gas is set to 0.01.)
F T (T soil (z)), F s (s(z)) and C i (z) are now all dependent on depth. T soil (z) and s(z) are the simulated layered soil temperature and soil moisture content and C i (z) is the simulated soil carbon content for each layer and pool i. The additional reduction of respiration with depth is exponential, with τ resp an empirical parameter (in m −1 ) that controls the magnitude of the reduction. The larger 405 the value of τ resp , the more inhibited the respiration is with increasing depth. When nitrogen is limiting, the respiration of the DPM and RPM pools are reduced by a factor of F N (N (z)) which is also now a function of depth and dependent on the available nitrogen in the relevant layer. M i and I i are also calculated as a function of depth based on their relationship with respiration.
The vertical mixing of each soil nitrogen pool follows that of the soil carbon pools: I tot is the total immobilisation in kg N m −2 s −1 in each layer. D(z) is the diffusivity in m 2 s −1 415 and varies both spatially and with depth (Burke et al., 2017). The litter inputs are distributed so that they decline exponentially with depth, with an e-folding depth of 0.2 m. Mineralised gas emissions are calculated for each soil layer, based on the balance of mineralisation and immobilisation in that layer (i.e. Equation 39 is repeated for every layer).

420
The inorganic nitrogen pool is the sum of deposition, fixation, immobilisation losses, mineralisation inputs, gridbox mean plant uptake and inorganic N losses through leaching and gaseous emission.
For the non-layered case (JULES-CN), these terms are simply added together: 45) where N in is the inorganic nitrogen in kg N m −2 . N dep is prescribed nitrogen deposition in kg N The loss of nitrogen from the inorganic pool is a function of leaching (N leach ) and an additional turnover (N turnover ).
where γ n is a tunable parameter. This additional term represents missing processes relating to the gaseous loss of inorganic nitrogen and limits the effective mineral N pool size. Without this additional turnover available N may increase excessively, potentially due to excessive biological fixation in regions that are generally unlimited. In the current model configuration this parameter is set to 1.0 (360 day − 1) such that the whole pool turns over once every model year. This results in an 435 effective saturation limit of 0.002 KgN m-2 consistent with Parton et al. (1993).
The leaching of nitrogen through the profile is assumed to be a function of the net flux of moisture through the soil profile, the concentration of inorganic N, and a parameter (β) representing the effective solubility of nitrogen. β has a value of 0.1. based on the sorption buffer coefficient of Ammonia (Gerber et al., 2010) although here it represents the sorption of all inorganic nitrogen species.

440
N leach = β(N in /θ 1m )Q subsurf ace (47) where θ 1m is the soil water content in the top 1m of soil in kg m −2 (so the inorganic nitrogen is assumed to occupy the top 1m of soil), and Q subsurf ace is the total subsurface runoff in kg m −2 s −1 .

Vertical discretisation of inorganic nitrogen
In JULES-CN layered , there is an inorganic nitrogen pool in each soil layer. The dynamics are very 445 similar to Equation 45, but every component now varies with depth, so: The net mineralisation flux (M net (z)) is the difference between M tot (z) and I tot (z) reduced by N gas (z) from each layer (see Section 3.2.1). Nitrogen deposition (N dep (z)) is added to the uppermost soil layer. Inputs from plant nitrogen fixation from PFT i are distributed according to the root 450 profile of the plants, i.e.
where f root (z) is the volumetric root fraction at a given depth.
Turnover (N turnover (z)) occurs in each layer, but with an exponential decay with depth as for the soil decomposition, which empirically represents the factors that reduce the soil activity with depth. The turnover term thus becomes: This leaves two terms in Equation 48: the plant uptake term ( ) and the N f lux term, which replaces the leaching term from Equation 45. These have a more process-based representation in the layered case. Plants cannot access all the inorganic nitrogen. Firstly, we assume that if a soil 460 layer is frozen then plants cannot uptake any of the nitrogen in that layer. Secondly, we assume that they only have direct access to a certain fraction of the soil, according to their root fraction, f root,i (which reduces with depth). So for each PFT, i, there is an 'available' inorganic nitrogen pool, which at equilibrium is as follows: As nitrogen is taken up from the available pool around the roots, there will be a delay in the area around the roots getting 're-filled'. We assume that it is constantly diffusing back to the equilibrium state where the concentration is constant both horizontally and vertically within each layer, and thus after the extraction on a particular TRIFFID timestep we update the available N pool according to:

470
where γ dif f is the rate of diffusion back to the equilibrium, set by default to 100 [360 day] −1 . Any fixation goes directly into the available pool, and other fluxes are simply added according to the ratio of the available to total inorganic N pools at equilibrium (thus the available pool would always follow Equation51 were it not for the fixation and uptake by plants). Plant uptake is extracted entirely from the available N pool, and the dependence on depth is according to the same profile as the available 475 N, i.e.
Leaching is now done in a process-based manner, where the inorganic N is transported through the soil profile by the soil water fluxes. Thus in Equation 48 we have the following term: where θ is the soil water content in kg m −2 and W f lux is the flow rate of the water through the soil in kg m −2 s −1 . Multiplying by dz n gives the change in N content for each layer, n. The total leaching is then the sum of all nitrogen that leaves the soil by lateral runoff or out of the bottom soil layer.

Historical simulations
Global transient simulations were carried out following the protocol for the S2 experiments in 485 TRENDY (Sitch et al., 2015). Forcing consisted of time-varying CO 2 , and climate from the CRU-NCEP data-set (v4, 1901, Viovy N. 2011. The fraction of 17 https://doi.org/10.5194/gmd-2020-205 Preprint. Discussion started: 24 July 2020 c Author(s) 2020. CC BY 4.0 License. agriculture in each grid cell (Hurtt et al., 2011) was set to the pre-industrial value. Nitrogen deposition was time-varying and was taken from a ACCMIP multi-model data set interpolated to annual fields (Lamarque et al., 2013). The model resolution was N96 (1.875 • longitude x 1.25 • latitude).

490
Results from three different JULES model configurations are presented here: -JULES-C represents the soil and vegetation carbon cycle in a manner comparable with HadGEM2-ES (Jones et al., 2011).
-JULES-CN is the nitrogen enabled version of JULES-C.
-JULES-CN layered is a version of JULES-CN with vertically discretised soil biochemistry.

495
In each case five PFTs were used: broadleaf trees, needleleaf trees (NET), C3 and C4 grass and shrubs. Plant competition was allowed, with TRIFFID updating vegetation fractions on a 10 day time step. The sole difference between JULES-C and JULES-CN is the inclusion of the nitrogen cycle. JULES-CN layered additionally has vertically discretised soil biogeochemistry. There are no differences in any of the shared model parameters.

500
The simulations were initialised using pre-industrial conditions. They were spun up by repeating the time period 1860-1870 until the change was less than 0.01 % decade −1 globally. The soil carbon distribution in JULES-CN layered is particularly slow to reach equilibrium. Therefore the 'modified accelerated decomposition' technique (modified-AD) described by Koven et al. (2013) was used to spin the soil carbon in these versions up to an initial pre-industrial equilibrium distribution (Burke 505 et al., 2017). Further spin up was then carried out for these layered models using repeated preindustrial conditions until the change in soil carbon was less than 0.01 % decade −1 globally. It should be noted that neither transient land-use change or fertiliser were included in any of these simulations.

510
This paper focuses on the differences between selected model configurations introduced by including the explicit nitrogen cycle within JULES. When available, we additionally use any observational based estimates to evaluate the quality of the simulations. First a broad-brush comparison between JULES-CN, JULES-C and JULES-CN layered is made. This is followed by a more of N inputs in European forests (Dise et al., 2009) and range between 59 and 118 T gN/yr −1 in the available observations (Boyer et al., 2006;Galloway et al., 2004;Gruber and Galloway, 2008).
Nitrogen uptake and net mineralisation are relatively high and are fairly similar in magnitude imply- In the model the soil carbon decomposition can be limited when the nitrogen available in the soil is less than the nitrogen required by decomposition. This process does not play a major role in our simulations.

Ecosystem residence times 560
The zonal total soil and vegetation carbon stocks are shown in Figure 7. The vegetation carbon is very similar for both JULES-C and JULES-CN as expected from Figure 4 and is consistent with the available observations. There are some differences in the soil carbon in the northern high latitudes with JULES-CN having slightly less soil carbon than JULES-C. This is a consequence of the higher nitrogen limitation on JULES-CN leading to less litter fall and subsequently less soil carbon. The 565 ecosystem residence time is defined as the total ecosystem carbon divided by the GPP. This is shown in Figure 8 for the biomes introduced in Figure 5. These residence times have been estimated annually on a grid cell by grid cell basis and then aggregated to a multi-annual mean per biome. The observational values were derived in a similar way using spatial data from Carvalhais et al. (2014).
In general the residence times are fairly similar for JULES-C and JULES-CN except for the tundra 570 biome where the residence time for JULES-CN is much longer than that for JULES-C. This is an over estimation of the residence time for this biome, however, JULES-CN is missing some key permafrost processes which will lead to an improvement (see Section 3.2.1).
The soil organic nitrogen ( Figure 9) shows a similar distribution to the soil carbon ( Figure 7) 575 reflecting the relatively consistent C to N ratio of the soil within the model. The observed soil organic nitrogen content is slightly higher at all latitudes than simulated by JULES-CN particularly in the northern tundra region.

Carbon-use efficiency
Carbon use efficiency (CUE) is defined as the ratio of net carbon gain to gross carbon assimila-580 tion during a given period (NPP/GPP). This represents the capacity of the plants to allocate carbon from photosynthesis to the terrestrial biomass. In the model nitrogen limitation restricts the ability of plants to allocate carbon and reduces the carbon use efficiency. Figure  an increase in CUE in both configurations, mainly caused by CO 2 fertilisation, but this is limited by nitrogen in the JULES-CN configuration.

Net ecosystem exchange
A key measure of a land carbon cycle model is how well it simulates the temporal variation of the 605 land carbon sink, which is the difference between Net Ecosystem Exchange (NEE) and the flux of carbon to the atmosphere from land-use change. The interannual variability in the sink is dominated by the variability of NEE, which is itself correlated with the magnitude of the temperature-carbon cycle feedback in the tropics (Cox et al., 2013). As a result, simulation of NEE variability is highly relevant to climate-carbon cycle projections (Wenzel et al., 2016).
The observations and both of the models show an upward trend in NEE but with very significant 620 interannual variability (Figure 12). Due to nitrogen limitations on CO 2 fertilization, mean NEE in JULES-CN (1.66 Pg C/yr) is lower than in JULES-C (2.06 Pg C/yr), and also lower than the 21 https://doi.org/10.5194/gmd-2020-205 Preprint. estimate from GCP (2.11 Pg C/yr). However, JULES-CN outperforms JULES-C on all of the other key metrics of the NEE variation. JULES-CN produces a smaller but much more realistic trend in NEE, and a smaller and more realistic interannual variability about that trend (see Table 5.2.3).

625
The correlation coefficient for NEE between the JULES-CN and GCP estimates (r=0.71) is also improved compared to JULES-C (r=0.63). There remains a significant underestimate of NEE in the years following the Pinatubo volcanic eruption, most likely due to the neglect of diffuse-radiation fertilization in these versions of JULES (Mercado et al., 2009). However, it is especially notable that JULES-CN significantly reduces the systematic overestimate of NEE seen in JULES-C during 630 extended La Nina periods, such as the years centred around 1974 and 2000 ( Figure 12).

Impact of vertical discretisation of soil biochemistry
Over the tropics and southern latitudes, JULES-CN layered is very comparable to JULES-CN. The main differences occur in the northern regions where there is soil freezing -adding vertically discretised soil biogeochemistry to JULES-CN has the most impact in the northern high latitudes. The 635 soil in JULES-CN layered has more inorganic nitrogen ( Figure 4) but it is not all accessible for the plants to uptake because the nitrogen uptake is limited by frozen soil. This means that in regions with frozen soil JULES-CN layered is slightly more nitrogen limited than JULES-CN ( Figure 5). Globally JULES-CN layered is possible also slightly more limited than the observations suggest ( Figure 5).
The biggest difference between JULES-CN layered and JULES-CN is for the boreal and coniferous 640 biome where the response ratio (potential NPP/ achieved NPP) for JULES-CN is 1.32 of that of JULES-CN layered is 1.48.
This additional limitation of nitrogen uptake caused by frozen soils means that JULES-CN layered has less total vegetation. However it also has more soil organic carbon and soil organic nitrogen

Conclusions
In this paper we have documented a model to quantify the impact of coupling the nitrogen cycle with the carbon cycle in a fully dynamic vegetation model. In this model, nitrogen limitation affects NPP 655 and how the carbon is allocated but it only indirectly affects the photosynthesis via leaf area development. This enables the carbon use efficiency (ratio of net carbon gain to gross carbon assimilation) to respond to changing nitrogen availability. Since the CUE affects the ability of the land surface to uptake carbon in a changing climate, this will impact carbon budgets under future projections of climate change. This scheme (based on JULES-CN) has been included within UKESM1 (Sellar et al., 2019).

660
Overall the nitrogen enabled configuration of JULES -JULES-CN -produces a more realistic trend in the net ecosystem exchange (NEE) and the interannual variability of NEE about that trend.
It also produces an improved estimate of NPP in the northern high latitudes. For other regions and diagnostics the simulation of present-day state and behaviour is not substantially different between 665 JULES-C and the nitrogen-enabled configuration, JULES-C. This is largely because JULES-C has been tuned to replicate observed carbon stores and fluxes and therefore implicitly includes a level of nitrogen availability. What JULES-C lacks is a mechanism for this to change substantially in time -either under more limiting conditions as elevated CO 2 outpaces demand for nutrients (e.g.Zaehle (2013b)), or under conditions of increased nitrogen availability due to anthropogenic deposition or 670 climate-induced mineralisation (Meyerholt et al., 2020b;Zaehle and Dalmonech, 2011). The response of the nitrogen cycle in JULES under changes in climate and CO 2 conditions-which will be affected by nutrient limitations-will be quantified and assessed in subsequent work.
An extended version of the nitrogen-enabled model model additionally includes the vertical dis-675 cretisation of the soil biogeochemistry model. This configuration improves the ecosystem residence times in the tundra. This more detailed representation of permafrost biogeochemistry in the northern high latitudes will used to understand the impact of the coupled carbon and nitrogen cycle on the permafrost carbon feedback.