JULES-CN: a coupled terrestrial carbon–nitrogen scheme (JULES vn5.1)

Abstract

Understanding future changes in the terrestrial carbon cycle is important for reliable projections of climate change and impacts on ecosystems. It is well known that nitrogen (N) could limit plants' response to increased atmospheric carbon dioxide and it is therefore important to include a representation of the N cycle in Earth system models. Here we present the implementation of the terrestrial nitrogen cycle in the Joint UK Land Environment Simulator (JULES) – the land surface scheme of the UK Earth System Model (UKESM). Two configurations are discussed – the first one (JULES-CN) has a bulk soil biogeochemical model and the second one is a development configuration that resolves the soil biogeochemistry with depth (JULES-CN${}_{\text{layer}}$). In JULES the nitrogen (N) cycle is based on the existing carbon (C) cycle and represents all the key terrestrial N processes in a parsimonious way. Biological N fixation is dependent on net primary productivity, and N deposition is specified as an external input. Nitrogen leaves the vegetation and soil system via leaching and a bulk gas loss term. Nutrient limitation reduces carbon-use efficiency (CUE – ratio of net to gross primary productivity) and can slow soil decomposition. We show that ecosystem level N limitation of net primary productivity (quantified in the model by the ratio of the potential amount of C that can be allocated to growth and spreading of the vegetation compared with the actual amount achieved in its natural state) falls at the lower end of the observational estimates in forests (approximately 1.0 in the model compared with 1.01 to 1.38 in the observations). The model shows more N limitation in the tropical savanna and tundra biomes, consistent with the available observations. Simulated C and N pools and fluxes are comparable to the limited available observations and model-derived estimates. The introduction of an N cycle improves the representation of interannual variability of global net ecosystem exchange, which was more pronounced in the C-cycle-only versions of JULES (JULES-C) than shown in estimates from the Global Carbon Project. It also reduces the present-day CUE from a global mean value of 0.45 for JULES-C to 0.41 for JULES-CN and 0.40 for JULES-CN${}_{\text{layer}}$, all of which fall within the observational range. The N cycle also alters the response of the C fluxes over the 20th century and limits the CO${}_{\mathrm{2}}$ fertilisation effect, such that the simulated current-day land C sink is reduced by about 0.5 Pg C yr${}^{-\mathrm{1}}$ compared to the version with no N limitation. JULES-CN${}_{\text{layer}}$ additionally improves the representation of soil biogeochemistry, including turnover times in the northern high latitudes. The inclusion of a prognostic land N scheme marks a step forward in functionality and realism for the JULES and UKESM models.

1 Introduction

Terrestrial ecosystems absorb around 25 % of anthropogenic carbon emissions , and changes in the future land carbon (C) 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 biomass, which relies on sufficient availability of nitrogen (N) within the soil–plant system. Therefore, the availability of N impacts the land C sink, both in the present and in a higher atmospheric carbon dioxide (CO${}_{\mathrm{2}}$) future.

Nitrogen exists in the terrestrial system in organic and inorganic forms and is continually cycled. In a stable climate the external inputs to the coupled vegetation and soil system – biological N fixation and N deposition – are balanced by the losses from this system – N leaching and N gas loss associated primarily with denitrification processes. Depending on the nutrient status of the vegetation and soil, changes in the balance of the inputs and outputs of N can drive adjustments in vegetation biomass and soil organic matter. Within the system, organic N is transferred from the vegetation to the soil through the production of litter and disturbance. The litter decomposes into soil organic matter and in turn is mineralised into inorganic N. Both inorganic and organic N may become available for plant uptake, although the amount of organic N uptake by plants is small and typically not included in models .

Any increase in atmospheric CO${}_{\mathrm{2}}$ drives an increase in the land C uptake and hence an increase in the gross primary productivity (GPP). This results in an extra demand for N, which could potentially limit the increase in future C stocks. For example, suggest that in some areas N could limit any future C sink by up to 70 %. N cycling also tends to reduce the sensitivity of land C uptake to temperature. Warmer conditions lead to increased plant respiration and soil respiration, which tends to reduce the land C sink. However, the increased soil respiration also leads to accelerated N mineralisation and increased N availability to plants, which may provide a counteracting increase in GPP. This latter effect is absent from models that do not include a N cycle. As a result of neglecting these important effects, land surface models without an interactive N cycle tend to overestimate both CO${}_{\mathrm{2}}$ and temperature effects on the land C sink . In addition, climate projections assessed by the IPCC using CMIP5 Earth system models that lacked terrestrial carbon cycle have been shown to exhibit a major and systematic bias in their future projection of land carbon sink . An increasing number of land surface and climate models now include constraints on the land C sink caused by N limitation . In fact, recent simulations have generated a range of estimates for the sensitivity of the C cycle to N availability . For example, used a perturbed model ensemble to show that N limitation reduces both the projected future increase in land C store due to CO${}_{\mathrm{2}}$ fertilisation and the projected loss in land C caused by climate change. The inclusion of nitrogen cycle processes in many CMIP6 models has been a major advance . show how CMIP6 models have a much reduced spread in their simulation of airborne fraction than in CMIP5 and this is attributable to the inclusion of the N cycle in about half of these latest-generation models. However, process understanding and evaluation of these models is still in its infancy .

The purpose of this paper is to describe and evaluate the implementation of coupled C and N cycles within the Joint UK Land-Environment Simulator (JULES at vn5.1 – http://jules-lsm.github.io/vn5.1, last access: 9 April 2021). JULES is the land surface component of the UK Earth System Model (UKESM) . The addition of the N cycle to JULES described in this paper was carried out alongside other developments such as improved plant physiology and an extended number of plant functional types , an enhanced representation of surface exchange and hydrology and a new module for land management . These separate components combine 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 to capture the established first-order emergent response of N addition on growth, which translates into leaf area index (LAI) and biomass. Our approach is to simulate the large-scale role of N limitation on vegetation C use efficiency (CUE – ratio of net to gross primary productivity) and soil C turnover. This is achieved by extending the implicit representation of N in the existing dynamic vegetation and plant physiology modules to be fully interactive. At the core of surface exchange in JULES is a coupled stomatal conductance photosynthesis scheme with a dependency on the leaf N concentration. Similarly, plant maintenance respiration has a dependency on leaf, root, and stem N concentration . Implicit within JULES, even in simulations excluding the N cycle, is the parameterisation of plant tissue level N concentrations and associated allometry (discussed further in Sect.  and by ). Simulations with an interactive C cycle therefore assume that enough N is available to meet vegetation growth and turnover. Here, we simply limit growth if not enough N is available. To do this requires a full representation of the N cycle in the land surface, including a coupled soil C and N organic and soil inorganic N scheme.

At the ecosystem level, the C and N cycles are closely coupled, with each flux of C associated with a corresponding flux of N linked via the C to N ratios. In JULES, nutrient limitation operates through two mechanisms. Firstly, the vegetation cannot uptake as much C – any C that the plants cannot uptake is denoted excess C. Secondly, the decomposition of litter C is slowed because there is insufficient N present. This is achieved by explicitly representing the demand for N within the vegetation and soil modules and then reducing plant net primary productivity to match available nutrients. In the soil module an additional decomposition rate modifier is introduced that slows respiration by microbes to match available nutrients. The current structure of the TRIFFID dynamic vegetation model , in particular the fixed allometry and C allocation, is largely unchanged. As the aim of this scheme is to capture the impact on terrestrial C stores, N loss terms are aggregated and not speciated. The model's reduced uptake of vegetation C due to N limitation is designed to have only a minor impact on the GPP. The emergent impact of the N scheme is modelled by reducing net primary productivity (NPP) and hence the carbon use efficiency (CUE) of the vegetation. In reality, the C the plants are unable to use because of insufficient N (defined as $\mathrm{\Psi }$) becomes non-structural carbohydrates, root exudates, or biogenic volatile organic compounds . However, to simplify the carbon balance in JULES-CN, it is added to the autotrophic respiration.

A key assumption in the JULES representation of vegetation and common amongst complex Dynamic Global Vegetation Models (DGVMs) is of fixed plant stoichiometry (mass ratio of C to N atoms or C : N ratio). The implication is that leaf-level photosynthetic capacity does not vary with available N. This is consistent with field experiments enhancing N fertilisation that find increases in growth but no corresponding change in photosynthetic capacity . However, more recent analyses do make the link between nutrient availability and leaf level N concentrations (e.g. ). In general, models make different assumptions about the tightness of the coupling mechanism between the C and N cycles, leading to substantial uncertainty in their projections . Within the fully coupled Earth systems models used in the Coupled Climate Carbon Cycle Model Intercomparison Project (C4MIP) for quantifying C feedbacks only 4 out of 11 models include a N cycle representation, and only 2, of which JULES is one, include both N and dynamic vegetation . The representation of the N cycle in the full complexity Earth system models remains challenging, and there is clearly a need for simple models capturing the first-order responses. This is the first time a N cycle has been incorporated in JULES and it is expected to be improved and developed with time as the knowledge of how important processes can be represented in existing frameworks improves.

2 Introduction to JULES

JULES is the land surface component of the new UK community Earth system model, UKESM . JULES can also be run offline forced by observed meteorology globally, regionally, or at a single location. A full description of the main components of JULES is provided by and . In particular, JULES represents the surface energy balance, a dynamic (one-dimensional) snowpack model, vertical heat and water fluxes, soil freezing, large-scale hydrology, and C fluxes and C storage in both vegetation and soil. Typically, JULES represents four soil layers down to a total depth of 3 m. Within JULES, C stocks and fluxes in and between the soil and vegetation, along with competition between different vegetation types, are modelled by the Top-Down Representation of Interactive Foliage and Flora Including Dynamics (TRIFFID) . In this version of TRIFFID, five plant functional types (PFTs) are represented: broadleaf trees, needleleaf trees, C${}_{\mathrm{3}}$ grasses, C${}_{\mathrm{4}}$ grasses, and shrubs. The soil C model is based on the Roth-C model . Recently, and added a representation of permafrost soil processes to JULES, including a representation of the vertical distribution of soil organic C that we build upon here. JULES-C is the standard carbon cycle configuration (a configuration defines a specific set of switches and parameters) and was used in the Global Carbon Budget annual assessment in 2018 .

What follows is a description of the extension of the C cycle process modelled by the JULES-C configuration to include an interactive N cycle. This results in two new model configurations: JULES-CN and JULES-CN${}_{\text{layer}}$. The soil biogeochemistry is represented by a single bulk layer in JULES-CN, whereas it varies as a function of depth in JULES-CN${}_{\text{layer}}$. As standard, JULES-C includes an implicit representation of N that has been extended to be fully interactive. The N cycle processes are added to the TRIFFID dynamic vegetation and Roth-C soil C models. For clarity we include a full description of the C and N cycle, including the existing TRIFFID and Roth-C models, and highlight where and how their processes have been modified.

Table 1

Default values of PFT-specific parameters for allometry, allocation, and vegetation N and C stoichiometry in the JULES-CN and JULES-CN${}_{\text{layer}}$ configurations. The subscript ($i$) is present to show that it is a PFT-specific value; $n_{l\mathrm{0},i}$ is the N concentration at the top of the canopy but is shown here as $\mathrm{1}/n_{l\mathrm{0},i}$ so that it is comparable to expected C : N ratios from the literature.

Symbol (units)	Definition	Broadleaf	Needleleaf	C${}_{\mathrm{3}}$	C${}_{\mathrm{4}}$	Shrub
		tree	tree	grass	grass
${\mathit{\sigma }}_{\text{l},i}$ (kg C m${}^{-\mathrm{2}}$)	Specific density of leaf C	0.0375	0.1000	0.0250	0.0500	0.0500
$a_{\text{wl},i}$ (kg C m${}^{-\mathrm{2}}$)	Allometric coefficient	0.65	0.65	0.005	0.005	0.10
$a_{\text{ws},i}$ (–)	Ratio total C to respiring stem C	10.0	10.0	1.0	1.0	10.0
$b_{\text{wl},i}$ (–)	Allometric exponent	1.667	1.667	1.667	1.667	1.667
${\mathit{\eta }}_{\text{sl},i}$ (kg C m${}^{-\mathrm{2}}$ per unit LAI)	Live stemwood coefficient	0.01	0.01	0.01	0.01	0.01
${\mathit{\mu }}_{\text{rl},i}$ (–)	Ratio root N to top leaf N	1.0	1.0	1.0	1.0	1.0
${\mathit{\mu }}_{\text{sl},i}$ (–)	Ratio stem N to top leaf N	0.1	0.1	1.0	1.0	0.1
$\mathrm{1}/n_{l\mathrm{0},i}$ ((kg C) (kg N)${}^{-\mathrm{1}}$)	C : N ratio at canopy top	21.7	30.3	13.7	16.67	16.67
$k_{\text{n},i}$ (–)	N profile coefficient	0.78	0.78	0.78	0.78	0.78
${\mathit{\kappa }}_{\text{r},i}$ (–)	Root N retranslocation coef.	0.2	0.2	0.2	0.2	0.2
${\mathit{\kappa }}_{\text{l},i}$ (–)	Leaf N retranslocation coef.	0.5	0.5	0.5	0.5	0.5
${\mathcal{L}}_{\text{min},i}$ (–)	Minimum balanced LAI	3.0	3.0	1.0	1.0	1.0
${\mathcal{L}}_{\text{max},i}$ (–)	Maximum balanced LAI	9.0	9.0	4.0	4.0	4.0
$f_{\text{DPM},i}$ (–)	Decomposable litter fraction	0.25	0.25	0.67	0.67	0.33

3 JULES developments

3.1 Vegetation C and N

The TRIFFID dynamic vegetation model provides the core of the vegetation module . TRIFFID represents the vegetation cover at each location in terms of the fractional area covered and the leaf area index (LAI) and canopy height of each plant functional type (PFT). In JULES the C fluxes are calculated at the model time step (typically 0.5–1 h) prior to any N limitation (if configured). These fluxes are then aggregated to the time step required for running TRIFFID (once every 10 d in the current implementation) so that allocation of C can take place. TRIFFID employs fixed allometry such that the split of vegetation carbon between leaf, root, and stem is defined by a single state prognostic variable that defines the total biomass. Biomass density increases via growth and is reduced by litter production and competition with other PFTs . Biomass can also increase by spreading through an increase in covered area. N limitation reduces growth and spreading such that the change in vegetation N cannot exceed the N uptake rate.

This section documents the vegetation model, starting with the structure of the vegetation (Sect. ) and including the additional complexity of labile N (Sect. ). The following subsection describes how growth and spreading is limited by N availability (Sect. ). The final subsection describes how vegetation C and N is turned over by disturbance and competition and aggregated from PFTs to the grid cell level (Sect. ). Biological N fixation is input directly into the soil inorganic N pool and is described later in Sect. .

3.1.1 Vegetation structure

The mean canopy height per PFT $i$ is converted via allometric equations into a maximum or balanced leaf area index for each PFT (${\mathcal{L}}_{\text{b},i}$ in m${}^{\mathrm{2}}$ m${}^{-\mathrm{2}}$). ${\mathcal{L}}_{\text{b},i}$ is the prognostic variable used in JULES to describe the vegetation and is functionally the equivalent of the potential leaf area. Given ${\mathcal{L}}_{\text{b},i}$, leaf, root, and wood pools are diagnosed for each PFT as introduced in . The balanced leaf area index is updated interactively following the C balance and is coupled to the surface exchange via surface albedo, roughness and heat capacity. This section is included to fully document the new scheme, but the equations can also be found in .

Figure 1

Stoichiometry of the vegetation nitrogen pools as a function of canopy height for individual PFTs at full leaf. Leaf N concentration are defined at the canopy level and are higher than those for the top leaf. The grey region shows the defined range of canopy height within the model. Note that both the $x$ and $y$ scales are very different for each PFT.

[Figure omitted. See PDF]

The vegetation C density per PFT ($C_{\text{v},i}$ in kg C m${}^{-\mathrm{2}}$) can be separated into leaf ($L_{\text{c},i}$ in kg C m${}^{-\mathrm{2}}$), fine root ($R_{\text{c},i}$ in kg C m${}^{-\mathrm{2}}$), and total stem plus coarse root ($W_{\text{c},i}$ in kg C m${}^{-\mathrm{2}}$) pools, each of which is related allometrically to the balanced leaf area (${\mathcal{L}}_{\text{b},i}$). Each component is then related to ${\mathcal{L}}_{\text{b},i}$. Root C is set equal to leaf C, which is itself a linear function of ${\mathcal{L}}_{\text{b},i}$, and total stem C is related to ${\mathcal{L}}_{\text{b},i}$ by a power law : $$\begin{array}{}\text{1}& {\displaystyle }& {\displaystyle }{C}_{\text{v},i}=L_{\text{c},i}+R_{\text{c},i}+W_{\text{c},i},\text{2}& {\displaystyle }& {\displaystyle }{L}_{\text{c},i}={\mathit{\sigma }}_{\text{l},i}{\mathcal{L}}_{\text{b},i},\text{3}& {\displaystyle }& {\displaystyle }{R}_{\text{c},i}=L_{\text{c},i},\text{4}& {\displaystyle }& {\displaystyle }{W}_{\text{c},i}=a_{\text{wl},i}({\mathcal{L}}_{\text{b},i}{)}^{b_{\text{wl},i}},\end{array}$$ where ${\mathit{\sigma }}_{\text{l},i}$ (kg C m${}^{-\mathrm{2}}$), $a_{\text{wl},i}$ (kg C m${}^{-\mathrm{2}}$), and $b_{\text{wl},i}$ (dimensionless) are PFT-dependent allometric parameters defined in Table . By definition ${\mathcal{L}}_{\text{b},i}$ does not have an explicit seasonal cycle but responds to changes in the vegetation C on both short (seasonal) and long (centennial) timescales. A high ${\mathcal{L}}_{\text{b},i}$ is related to a high C density and tall canopies. It should be noted that leaf seasonality is represented by a separate phenology model and is not directly affected by N availability. TRIFFID combines Eq. () with a “pipe model” approach to obtain the canopy height for PFT $i$ ($h_i$ in m): 5$$h_i={\displaystyle \frac{W_{\text{c},i}}{a_{\text{ws},i}{\mathit{\eta }}_{\text{sl},i}}}{\left({\displaystyle \frac{a_{\text{wl},i}}{W_{\text{c},i}}}\right)}^{\mathrm{1}/b_{\text{wl},i}},$$ where ${\mathit{\eta }}_{\text{sl},i}$ (kg C m${}^{-\mathrm{2}}$ per unit LAI) relates respiring stem to leaf C (Table ) and $a_{\text{ws},i}$ is the ratio of total stem C to respiring stem C. We can combine Eqs. () and () to relate (${\mathcal{L}}_{\text{b},i}$) to canopy height ($h_i$) and these two variables can be used interchangeably to describe the state of the vegetation. During a simulation the C pools are updated interactively and the canopy height and balanced leaf area diagnosed for each PFT. This representation allows changes in vegetation C to feedback to surface exchange.

The root and total stem N pools are defined using stoichiometric relationships as a function of the C 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 N pools as part of the new scheme: $$\begin{array}{}\text{6}& {\displaystyle }& {\displaystyle }{R}_{\text{n},i}={\mathit{\mu }}_{\text{rl},i}{n}_{l\mathrm{0},i}{R}_{\text{c},i},\text{7}& {\displaystyle }& {\displaystyle }{W}_{\text{n},i}={\mathit{\mu }}_{\text{sl},i}{n}_{l\mathrm{0},i}{W}_{\text{c},i},\end{array}$$ where ${\mathit{\mu }}_{\text{rl},i}$ and ${\mathit{\mu }}_{\text{sl},i}$ are dimensionless stoichiometric parameters linking the top leaf N concentration ($n_{l\mathrm{0},i}$ in kg N kg C ${}^{-\mathrm{1}}$) to the total stem and root N concentration via $n_{l\mathrm{0},i}$. The leaf N pool ($L_{\text{n},i}$ in kg N m${}^{-\mathrm{2}}$) has an additional dependency on phenological state (Sect. ) and assumed distribution of N in the canopy. Following Eq. (), the total vegetation N store per PFT ($N_{\text{v},i}$ in kg N m${}^{-\mathrm{2}}$) is given by 8$$N_{\text{v},i}=L_{\text{n},i}+R_{\text{n},i}+W_{\text{n},i}.$$

The C : N ratio of the root and stem pools are fixed 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 C for woody PFTs (Fig. ). This is due to the relatively greater proportion of stem C at higher biomass. Grasses show less variation with biomass due to their comparatively small amount of structural C relative to leaf area, which also results in woody PFTs having higher C : N ratios. Equations (1)–(8) show that the total vegetation N increases with canopy height and biomass (Fig. ).

Figure 2

Total vegetation N along with N pools of leaf, root, and wood as a function of canopy height for individual PFTs at full leaf. The grey region shows the defined range of canopy height within the model. Note that both the $x$ and $y$ scales are very different for each PFT.

[Figure omitted. See PDF]

3.1.2 Labile C and N: phenology and mobilisation

The total leaf C pool per PFT ($L_{\text{c},i}$, Eq. ) varies allometrically with the vegetation C state on both short (seasonal) and long (centennial) timescales but not with changes in phenological state. Implicit within TRIFFID is a labile leaf C pool that acts as a reserve of C during spring and a store during fall. $L_{\text{c},i}$ therefore includes a labile pool from which C can be mobilised during leaf out plus an allocated pool representing the actual LAI. The labile pool is zero at full leaf out and at the allometrically defined maximum during the no-leaf period. As part of the N coupling we introduce the ability for plants to retranslocate some of the allocated N to the labile N pool according to the phenology. The new parameterisation of retranslocation and labile N is therefore dependent on the leaf phenological state and the fixed stoichiometry. In JULES, leaf phenology is controlled by a second state variable ($p_i$), which relates the LAI (${\mathcal{L}}_i$) at any moment in time to the balanced leaf area index (${\mathcal{L}}_{\text{b},i}$):

9$${\mathcal{L}}_i=p_i{\mathcal{L}}_{\text{b},i},$$ where $p_i$ is a scalar between 0 and 1 that describes the phenological state of the system . For evergreen plants $p_i$ is a constant of 1. The two state variables ${\mathcal{L}}_{\text{b},i}$ and $p_i$ combine to define the vegetation phenological state for each PFT $i$. Using the phenological state we extend the equivalent approach to leaf C such that the leaf N pool ($L_{\text{n},i}$) has fixed allometry dependent on the phenological state and the magnitude of leaf retranslocation. We introduce this simple parameterisation under the assumption that higher leaf retranslocation during autumn implies a higher labile N store. The leaf N pool therefore becomes 10$$L_{\text{n},i}=p_i{n}_{\text{lc},i}{L}_{\text{c},i}+(\mathrm{1}-p_i)\left({\displaystyle \frac{\mathrm{1}+{\mathit{\kappa }}_{\text{l},i}}{\mathrm{2}}}\right)n_{\text{lc},i}{L}_{\text{c},i},$$ where ${\mathit{\kappa }}_{\text{l},i}$ is the dimensionless leaf N retranslocation coefficient and $n_{\text{lc},i}$ is the mean canopy N concentration (defined in Eq.  in kg N kg C ${}^{-\mathrm{1}}$). Here ${\mathit{\kappa }}_{\text{l},i}$ is set to 0.5 for all PFTs . The formulation of the labile pool in this configuration means that around half of the N required for full leaf out is taken from leaf retranslocation with a further quarter acquired during the dormant phase while the rest is acquired during the leaf-out period.

JULES assumes a process-based scaling up of leaf-level photosynthesis to the the canopy level. In both the JULES-CN and JULES-CN${}_{\text{layer}}$ configurations, to be consistent with the JULES-C model, we assume a multi-level canopy with leaf N decreasing exponentially through the canopy (CanRadMod 5). The plant physiology routines uses this assumed distribution to calculate penetration through the canopy and photosynthesis on individual layers before scaling back to the canopy . In the application here, we use this distribution to be fully consistent with the physiology. The vertical distribution of leaf N concentration ($n_{\text{lc},i}\left(d\right)$ in kg N kg C ${}^{-\mathrm{1}}$) in the canopy is described by 11$$n_{\text{lc},i}\left(d\right)=n_{\text{l}\mathrm{0},i}{e}^{-k_{\text{n},i}d},$$ where $k_{\text{n},i}$ is a constant representing the profile of N density and $d$ represents the fraction of canopy above the layer. Based on observed N profiles in the Amazon basin , a value of 0.78 for $k_{\text{n},i}$ was found . Equation () is independent of leaf area and therefore equates to a constant of proportionality relating PFT-specific top leaf N to the mean canopy N concentration.

3.1.3 Vegetation growth and allocation

The previous section describe how the vegetation C ($C_{\text{v},i}$, Eq. ) and vegetation N ($N_{\text{v},i}$, Eq. ) for each PFT vary with vegetation size and phenological state. This section describes how growth and spreading are limited by available N. Growth is the increase in C density and spreading is the increase in vegetation cover from recruitment and reproduction.

Net primary productivity (NPP) in JULES-C is simply the difference between GPP and autotrophic respiration ($R_{\text{a}}$). In JULES-CN the potential NPP or NPP${}_{\text{pot}}$ is defined in the same way as the NPP in JULES-C before the explicit N cycle was included, i.e. the potential amount of C that can be allocated to growth ($g$) and spreading by TRIFFID. In JULES-CN and in order for the NPP to achieve its potential it needs to be able to uptake sufficient inorganic N. If not enough inorganic N is available, the system is N-limited, and an additional term is required in the C balance representing C that cannot be assimilated into the plant due to lack of available N ($\mathrm{\Psi }$ in kg C m${}^{-\mathrm{2}}$). A positive $\mathrm{\Psi }$ results in a reduction of carbon use efficiency (CUE).

The C balance per PFT $i$ is given by

12$${\displaystyle \frac{\text{d}{C}_{\text{v},i}}{\text{d}t}}=(\mathrm{1}-{\mathit{\lambda }}_i){\mathrm{\Pi }}_i-{\mathrm{\Lambda }}_{\text{c},i}-{\mathrm{\Psi }}_{\text{g},i},$$ where ${\mathrm{\Pi }}_i$ is the potential NPP per unit area of PFT in (kg C m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$ (prior to nutrient limitation) and ${\mathrm{\Lambda }}_{\text{c},i}$ (kg C m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$) is the PFT-specific litterfall rate (Sect. ). Any excess C from growth (${\mathrm{\Psi }}_{\text{g},i}$) is considered an additional plant respiration term and at the end of the TRIFFID time step is used to reduce the potential NPP for each PFT to its actual value. ${\mathit{\lambda }}_i$ is the coefficient for partitioning the NPP between growth and spreading. ${\mathit{\lambda }}_i$ is utilised in increasing the fractional coverage of the vegetation, and $(\mathrm{1}-{\mathit{\lambda }}_i)$ increases the C of the existing vegetated area. ${\mathit{\lambda }}_i$ is a function of the vegetation C which itself is a function of the balanced LAI for PFT $i$ (${\mathcal{L}}_{\text{b},i}$): 13$${\mathit{\lambda }}_i=\left\{\begin{array}{ll}\mathrm{1}& {\mathcal{L}}_{\text{b},i}>{\mathcal{L}}_{\text{max},i}\\ \frac{{\mathcal{L}}_{\text{b},i}-{\mathcal{L}}_{\text{min},i}}{{\mathcal{L}}_{\text{max},i}-{\mathcal{L}}_{\text{min},i}}& {\mathcal{L}}_{\text{min},i}<{\mathcal{L}}_{\text{b},i}\le {\mathcal{L}}_{\text{max},i}\\ \mathrm{0}& {\mathcal{L}}_{\text{b},i}\le {\mathcal{L}}_{\text{min},i}.\end{array}\right.$$

The equivalent N balance per PFT is given by 14$${\displaystyle \frac{\text{d}{N}_{\text{v},i}}{\text{d}t}}=(\mathrm{1}-{\mathit{\lambda }}_i){\mathrm{\Phi }}_i-{\mathrm{\Lambda }}_{\text{n},i},$$ where ${\mathrm{\Phi }}_i$ (kg N m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$) is the PFT-specific N uptake (see Eq. ) and $(\mathrm{1}-{\mathit{\lambda }}_i){\mathrm{\Phi }}_i$ is equal to ${\mathrm{\Phi }}_{\text{g},i}$, the N uptake available for growth. ${\mathrm{\Lambda }}_{\text{n},i}$ is the PFT N litter flux after taking into account the retranslocation of N from leaves or roots. The N available for spreading is a fraction ${\mathit{\lambda }}_i$ of the total available N with a fraction $(\mathrm{1}-{\mathit{\lambda }}_i)$ available for growth.

Litter is produced by the turnover of the leaf, wood and root pools for each PFT, defined as 15$${\mathrm{\Lambda }}_{\text{c},i}={\mathit{\gamma }}_{\text{l},i}{L}_{\text{c},i}+{\mathit{\gamma }}_{\text{r},i}{R}_{\text{c},i}+{\mathit{\gamma }}_{\text{w},i}{W}_{\text{c},i},$$ and 16$${\mathrm{\Lambda }}_{\text{n},i}=(\mathrm{1}-{\mathit{\kappa }}_{\text{l},i}){\mathit{\gamma }}_{\text{l},i}{L}_{\text{n},i}+(\mathrm{1}-{\mathit{\kappa }}_{\text{r},i}){\mathit{\gamma }}_{\text{r},i}{R}_{\text{n},i}+{\mathit{\gamma }}_{\text{w},i}{W}_{\text{n},i},$$ for litter C (${\mathrm{\Lambda }}_{\text{c},i}$ in kg C m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$) and litter N (${\mathrm{\Lambda }}_{\text{n},i}$ in kg N m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$), respectively. ${\mathit{\gamma }}_{\text{r},i}$ and ${\mathit{\gamma }}_{\text{w},i}$ are turnover rates in s${}^{-\mathrm{1}}$ (Table 6 of ). The leaf turnover rate (${\mathit{\gamma }}_{\text{l},i}$) is a temperature-dependent turnover rate consistent with the phenological state and defined in . The equivalent term for N allows for retranslocation of N from leaves into the labile store and a reduced N cost of maintaining fine roots. ${\mathit{\kappa }}_{\text{l},i}$ and ${\mathit{\kappa }}_{\text{r},i}$ are the dimensionless coefficients for the retranslocation of leaf and root N shown in Table  .

In JULES-CN the N available for plant uptake for each PFT $i$ ($N_{\text{avail},i}$ in kg N m${}^{-\mathrm{2}}$) is the inorganic soil N pool ($N_{\text{in}}$ in kg N m${}^{-\mathrm{2}}$) split equitably between the PFTs assuming there is no differential ability between PFTs to acquire N and the whole pool is available for uptake during the model time step. The available N in JULES-CN${}_{\text{layer}}$ is more complicated and takes into account the soil profile. This is discussed in Sect. .

Equations () and () have two remaining unknowns for each PFT: the plant N uptake for growth (${\mathrm{\Phi }}_{\text{g},i}$) and the excess C from growth (${\mathrm{\Psi }}_{\text{g},i}$). The litter fluxes are functions of the total vegetation pool and therefore can be solved at the same time. Solving for the case where ${\mathrm{\Psi }}_{\text{g},i}=\mathrm{0.0}$ gives the total vegetation N demand for growth. If the N demand is less than the available N in a given time step ($\mathrm{\Delta }t$) (${\mathrm{\Phi }}_{\text{g},i}<(\mathrm{1}-{\mathit{\lambda }}_i)N_{\text{avail},i}/\mathrm{\Delta }t$), then growth is unlimited and the fluxes can be updated accordingly. Where N is limiting, growth N uptake is set equal to the available N (${\mathrm{\Phi }}_{\text{g},i}=(\mathrm{1}-{\mathit{\lambda }}_i)N_{\text{avail},i}/\mathrm{\Delta }t$) and the excess C for growth ${\mathrm{\Psi }}_{\text{g},i}$ can be derived. Following the solution of $\frac{\text{d}{N}_{\text{v},i}}{\text{d}t}$ the C store and balanced LAI (${\mathcal{L}}_{\text{b},i}$) are updated and the leaf, root, and wood pools for each PFT can be derived following the allometric equations (Eqs. –).

The remaining proportion (${\mathit{\lambda }}_i$) of NPP and N is allocated to spreading. The N demand for spreading is equal to the C allocated to spreading scaled by the whole plant stoichiometry: 17$${\mathrm{\Phi }}_{\text{s},i}={\displaystyle \frac{N_{\text{v},i}}{C_{\text{v},i}}}\left({\mathrm{\Pi }}_i-{\displaystyle \frac{\text{d}{C}_{\text{v},i}}{\text{d}t}}-{\mathrm{\Psi }}_{\text{s},i}\right),$$ where ${\mathrm{\Psi }}_{\text{s},i}$ (or ${\mathit{\lambda }}_i{\mathrm{\Psi }}_i$) is the excess C term from spreading and $\frac{N_{\text{v},i}}{C_{\text{v},i}}$ is the inverse of the the whole plant C : N ratio. As with growth limitation, Eqs. () is first solved to find the N demand for spreading (${\mathrm{\Psi }}_{\text{s},i}=\mathrm{0.0}$). If the arising demand is less than the available N (${\mathrm{\Phi }}_{\text{s},i}<{\mathit{\lambda }}_i{N}_{\text{avail},i}/\mathrm{\Delta }t$), spreading is unlimited. If N demand is in excess of that available, the uptake is set equal to the available N flux (${\mathrm{\Phi }}_{\text{s},i}={\mathit{\lambda }}_i{N}_{\text{avail},i}/\mathrm{\Delta }t$) and the excess (${\mathrm{\Psi }}_{\text{s},i}$) assimilate is solved for.

Total excess C per PFT $i$ (${\mathrm{\Psi }}_i$) is therefore the combination of that from growth plus spreading: 18$${\mathrm{\Psi }}_i={\mathrm{\Psi }}_{\text{s},i}+{\mathrm{\Psi }}_{\text{g},i}.$$ Similarly, total N uptake per PFT $i$ (${\mathrm{\Phi }}_i$) is therefore the combination of N uptake used for growth plus N uptake used for spreading: 19$${\mathrm{\Phi }}_i={\mathrm{\Phi }}_{\text{s},i}+{\mathrm{\Phi }}_{\text{g},i}.$$

The PFT level N uptake and excess C are weighted by the fraction of coverage of each PFT in a grid cell ($v_i$) and summed to get the grid-averaged values: $$\begin{array}{}\text{20}& {\displaystyle }& {\displaystyle }\mathrm{\Phi }=\sum _i{v}_i{\mathrm{\Phi }}_i,\text{21}& {\displaystyle }& {\displaystyle }\mathrm{\Psi }=\sum _i{v}_i{\mathrm{\Psi }}_i.\end{array}$$

This excess C ($\mathrm{\Psi }$) is considered an additional plant respiration term and at the end of the TRIFFID time step is used to reduce the potential NPP to its actual value.

The C and N allocated to spreading allow the vegetation to expand onto bare ground. Where space is limiting the PFTs compete for space. The competition is handled in the Lotka–Volterra competition routines (see , for full details). N only indirectly affects competition through the PFT-specific allometric relationships. The competition code subsequently updates the fractional coverage of model PFTs ($v_i$).

3.1.4 Vegetation turnover and total litter production

The previous sections describe how N interacts to limit both growth and spreading of vegetation in the dynamic vegetation model. This final section describes the turnover of C and N through large-scale disturbance and competition.

Turnover is aggregated across PFTs to provide a grid box mean litter flux term to the soil biogeochemistry processes which acts at a grid box level. Total litter C (${\mathrm{\Lambda }}_{\text{c}}$, kg C m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$) is made-up of the area-weighted sum of the litter C from each PFT (${\mathrm{\Lambda }}_{\text{c},i}$), along with large-scale PFT-dependent disturbance rate, and a density dependent component from intra-PFT competition for space. Large-scale disturbance is implemented in TRIFFID as a constant disturbance rate per PFT and captures processes such as wind-throw and other mortality events. Density dependent litter production arises through competition for space with increased turnover when space is limiting and plants are competing for space and light.

22$${\mathrm{\Lambda }}_{\text{c}}=\sum _i{v}_i\left({\mathrm{\Lambda }}_{\text{c},i}+{\mathit{\gamma }}_{\text{v},i}{C}_{\text{v},i}+({\mathrm{\Pi }}_i-{\mathrm{\Psi }}_i)\sum _j{c}_{i,j}{v}_j\right),$$ where $c_{i,j}$ are the competition coefficients describing the effect of PFT $i$ on PFT $j$, ${\mathit{\gamma }}_{\text{v},i}$ is a large-scale disturbance term of PFT $i$, $v_i$ is the vegetation fraction of PFT $i$, ${\mathrm{\Pi }}_i$ is defined in Eq. () and ${\mathrm{\Psi }}_i$ in Eq. (). The effect of N limitation on the litter C flux is captured in the excess C term per PFT (${\mathrm{\Psi }}_i$). Similarly to the total litter C, total litter N (${\mathrm{\Lambda }}_{\text{n}}$, kg N m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$) is given by 23$${\mathrm{\Lambda }}_{\text{n}}=\sum _i{v}_i\left({\mathrm{\Lambda }}_{\text{n},i}+{\mathit{\gamma }}_{\text{v},i}{N}_{\text{v},i}+{\mathrm{\Phi }}_i\sum _j{c}_{i,j}{v}_j\right).$$ Both ${\mathrm{\Lambda }}_{\text{c}}$ and ${\mathrm{\Lambda }}_{\text{n}}$ vary according to the vegetation type and the relative amount of stem, leaves, and roots being turned over. This means that the C : N ratio also varies in time and space.

3.2 Soil biogeochemistry

Here we describe the addition of a prognostic soil N model for JULES-CN that extends the Roth-C soil C model used in JULES-C .

The original Roth-C soil C model represents four C pools ($p$) for each grid box. Plant litter input is split between two C pools of decomposable (DPM) and resistant (RPM) plant material, with the fraction that goes to each depending on the overlying vegetation PFT and parameterised via $f_{\text{DPM},i}$. Grasses provide a higher fraction of decomposable litter input and trees provide a higher fraction of resistant litter input. The other two C pools are microbial biomass (BIO) and long-lived humified (HUM) pools. The DPM and RPM pools can be characterised as representing litter and BIO and HUM as representing soil organic matter. C from decomposition of all of the pools is partly released to the atmosphere, and the remaining fraction (${\mathit{\beta }}_R$) enters the BIO and HUM pools. The C pools are updated according to $$\begin{array}{}\text{24}& {\displaystyle }& {\displaystyle \frac{\text{d}{C}_{\text{DPM}}}{\text{d}t}}=\sum _i\left(v_i{f}_{\text{DPM},i}{\mathrm{\Lambda }}_{\text{c},i}\right)-R_{\text{DPM}},\text{25}& {\displaystyle }& {\displaystyle \frac{\text{d}{C}_{\text{RPM}}}{\text{d}t}}=\sum _i\left(v_i(\mathrm{1}-f_{\text{DPM},i}){\mathrm{\Lambda }}_{\text{c},i}\right)-R_{\text{RPM}},\text{26}& {\displaystyle }& {\displaystyle \frac{\text{d}{C}_{\text{BIO}}}{\text{d}t}}=\mathrm{0.46}{\mathit{\beta }}_R{R}_{\text{tot}}-R_{\text{BIO}},\text{27}& {\displaystyle }& {\displaystyle \frac{\text{d}{C}_{\text{HUM}}}{\text{d}t}}=\mathrm{0.54}{\mathit{\beta }}_R{R}_{\text{tot}}-R_{\text{HUM}},\end{array}$$ where $t$ is the time in s, $C_p$ are the C pools in kg C m${}^{-\mathrm{2}}$ (where $p$ is one of DPM, RPM, BIO, HUM), ${\mathrm{\Lambda }}_{\text{c},i}$ is the litter input for PFT $i$ in kg C m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$ (term in brackets in Eq. ), $f_{\text{DPM},i}$ represents the fraction of litter from each PFT $i$ that goes into DPM with the rest $(\mathrm{1}-f_{\text{DPM},i})$ going into the RPM pool (dependent on amount of woody vegetation), and $R_{\text{tot}}$ is the total turnover in kg C m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$, where the $R_p$ represent the turnover of each C pool:

28$$R_{\text{tot}}=R_{\text{DPM}}+R_{\text{RPM}}+R_{\text{BIO}}+R_{\text{HUM}}.$$

The soil respiration to the atmosphere ($r_{\text{h}}$ in kg C m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$) is given by 29$$r_{\text{h}}=(\mathrm{1}-{\mathit{\beta }}_R)R_{\text{tot}},$$ where ${\mathit{\beta }}_R$ depends on soil clay content (clay in %) and ranges from 0.25 for a soil with no clay content to 0.15 for a clay soil: 30$${\mathit{\beta }}_R={\displaystyle \frac{\mathrm{1}}{\mathrm{4.09}+\mathrm{2.67}{e}^{(-\mathrm{0.079}\text{clay})}}}.$$

For each C pool there is an equivalent N pool, with the N pools following a similar structure to the C pools: $$\begin{array}{}\text{31}& {\displaystyle }& {\displaystyle \frac{\text{d}{N}_{\text{DPM}}}{\text{d}t}}=\sum _i\left(v_i{f}_{\text{DPM},i}{\mathrm{\Lambda }}_{\text{n},i}\right)-M_{\text{DPM}},\text{32}& {\displaystyle }& {\displaystyle \frac{\text{d}{N}_{\text{RPM}}}{\text{d}t}}=\sum _i\left(v_i(\mathrm{1}-f_{\text{DPM},i}){\mathrm{\Lambda }}_{\text{n},i}\right)-M_{\text{RPM}},\text{33}& {\displaystyle }& {\displaystyle \frac{\text{d}{N}_{\text{BIO}}}{\text{d}t}}=\mathrm{0.46}{I}_{\text{tot}}-M_{\text{BIO}},\text{34}& {\displaystyle }& {\displaystyle \frac{\text{d}{N}_{\text{HUM}}}{\text{d}t}}=\mathrm{0.54}{I}_{\text{tot}}-M_{\text{HUM}}.\end{array}$$

Inputs into the litter pools (DPM, RPM) are from the litter N flux (${\mathrm{\Lambda }}_{\text{n},i}$ in kg N m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$, Eq. ), and losses are determined by the pool-specific mineralisation of organic N into inorganic N ($M_{\text{p}}$ in kg N m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$). Following the framework of the Roth-C model, input into both the BIO and HUM N pools is from the total immobilisation of inorganic N into organic N ($I_{\text{tot}}$ in kg N m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$): 35$$I_{\text{tot}}=I_{\text{DPM}}+I_{\text{RPM}}+I_{\text{BIO}}+I_{\text{HUM}}.$$

For each soil C pool ($p$), the potential turnover – i.e. the turnover rate when the N in the system is not limiting – is given by ($R_{p,\text{pot}}$) 36$$R_{p,\text{pot}}=k_p{C}_p{F}_T\left(T_{\text{soil}}\right)F_{\mathit{\theta }}\left(\mathit{\theta }\right)F_v\left(v\right),$$ where the $k_p$ are fixed constants in s${}^{-\mathrm{1}}$ . The functions of temperature ($F_T\left(T_{\text{soil}}\right)$) and moisture ($F_{\mathit{\theta }}\left(\mathit{\theta }\right)$) depend on the temperature ($T_{\text{soil}}$) and moisture content ($\mathit{\theta }$) near the soil surface. The function $F_v\left(v\right)$ depends on the vegetation cover fraction ($v$) . The potential mineralisation of organic N when the system is not N-limited ($M_{p,\text{pot}}$) is related to the potential turnover rates by the C to N ratio of each pool (${\text{CN}}_p$): 37$$M_{p,\text{pot}}={\displaystyle \frac{R_{p,\text{pot}}}{{\text{CN}}_p}}.$$

Similarly, the potential immobilisation of inorganic N into the organic N pools ($I_{p,\text{pot}}$) is related to pool potential turnover ($R_{p,\text{pot}}$), the retained fraction of respiration (${\mathit{\beta }}_R$), and the C to N ratio of the destination pool in the decomposition chain: 38$$I_{p,\text{pot}}={\mathit{\beta }}_R{\displaystyle \frac{R_{p,\text{pot}}}{{\text{CN}}_{\text{soil}}}},$$ where ${\text{CN}}_{\text{soil}}$ is a model parameter that fixes the C to N ratios of the two destination soil organic pools (HUM and BIO) and has a default value of 10. The C to N ratio of the DPM and RPM litter pools is a function of litter quality and varies temporally and spatially depending on the contributions of the different PFTs within the grid cell. Potential mineralisation ($M_p$) and potential immobilisation ($I_p$) fluxes are defined before any N limitation is applied and take values that maintain the constant C : N ratio for the HUM and BIO pools.

When N is limiting, the turnover of the two litter pools (DPM and RPM) into the soil organic matter pools is additionally limited by the availability of N. 39$$R_p=k_p{C}_p{F}_T\left(T_{\text{soil}}\right)F_{\mathit{\theta }}\left(\mathit{\theta }\right)F_v\left(v\right)F_{\text{N}},$$ where $p$ is one of RPM or DPM. The nitrogen-limited mineralisation and immobilisation of the DPM and RPM pools (Eqs.  and ) are now effectively a function of $R_p$.

$F_{\text{N}}$ is the litter decomposition rate modifier and is given by the ratio of the N available in the soil to the N required by decomposition (Eq. ). $F_{\text{N}}$ is limited to a range of 0.0 to 1.0. When $F_{\text{N}}$ is equal to 1, the decomposition, mineralisation, and immobilisation take place at the potential rate and the system is not N-limited. Where $F_{\text{N}}$ is less than 1, the availability of N limits the decomposition of litter into soil organic matter. This limitation is because respiration is carried out by microbes that require sufficient N to convert the RPM and DPM pools into BIO and HUM pools. $F_{\text{N}}$ is given by 40$$F_{\text{N}}={\displaystyle \frac{(M_{\text{BIO}}+M_{\text{HUM}}-I_{\text{BIO}}-I_{\text{HUM}})\mathrm{\Delta }t+N_{\text{in}}}{(D_{\text{DPM}}+D_{\text{RPM}})\mathrm{\Delta }t}},$$ where $N_{\text{in}}$ is the total soil inorganic N pool (in kg N m${}^{-\mathrm{2}}$; discussed in Sect.  and defined in Eq. ) and $\mathrm{\Delta }t$ is the time step. $D_{\text{DPM}}$ and $D_{\text{RPM}}$ are the net demand associated with decomposition of each of the litter pools: 41$$D_p=I_{p,\text{pot}}-M_{p,\text{pot}},$$ where $p$ is one of RPM or DPM. This demand is always positive because the C to N ratio of soil is much less than the C to N ratio of the DPM and RPM pools. When the net demand is in excess of the available inorganic N, the system is N-limited and $F_{\text{N}}<\mathrm{1.0}$. This available N is mainly the net mineralised N from the turnover of BIO and HUM pools but also from the inorganic N pool. N limitation reduces the soil respiration, mineralisation, and immobilisation of the two litter pools (RPM and DPM). The C : N ratio of these two pools are variable in time and are represented as prognostic variables. The other two organic matter pools (BIO and HUM) always respire and are mineralised and immobilised at the potential rate (so $F_{\text{N}}$ is effectively 1.0).

If the net mineralisation is positive some of the N is emitted as gas, according to 42$$N_{\text{gas}}=f_{\text{gas}}(M_{\text{tot}}-I_{\text{tot}}),$$ where $N_{\text{gas}}$ is one component of the gas emission in kg N m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$, $f_{\text{gas}}$ is a parameter that sets the fraction of the N flux that is emitted as gas to the atmosphere. Following , it is assumed that 1 % of net mineralisation is emitted as gas ($f_{\text{gas}}$ is set to 0.01). $M_{\text{tot}}$ is the total mineralisation flux in kg N m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$: 43$$M_{\text{tot}}=M_{\text{DPM}}+M_{\text{RPM}}+M_{\text{BIO}}+M_{\text{HUM}}.$$ If pool sizes become too small, $N_{\text{gas}}$ could become negative to ensure N is conserved.

3.2.1 Vertical discretisation

The vertical discretisation of the soil C and N follows . There is a set of four soil C and N pools (DPM, RPM, BIO, HUM) in every soil model layer. As in the turnover rate is determined for each soil layer depending on the temperature, moisture conditions and N availability in that layer. An extra reduction of turnover with depth ($z$) is included to account for factors that are currently missing in the model such as priming effects, anoxia, soil mineral surface, and aggregate stabilisation. The potential turnover of each layer is given by

44$$R_{p,\text{pot}}\left(z\right)=k_p{C}_p\left(z\right)F_T\left(T_{\text{soil}}\right(z\left)\right)F_{\mathit{\theta }}\left(\mathit{\theta }\right(z\left)\right)F_v\left(v\right)\exp (-{\mathit{\xi }}_{\text{resp}}z).$$ $F_T\left(T_{\text{soil}}\right(z\left)\right)$, $F_{\mathit{\theta }}\left(\mathit{\theta }\right(z\left)\right)$ and $C_p\left(z\right)$ are now all dependent on depth. $T_{\text{soil}}\left(z\right)$ and $\mathit{\theta }\left(z\right)$ are the simulated layered soil temperature and soil moisture content, and $C_p\left(z\right)$ is the simulated soil C content for each layer and pool $p$. The additional reduction of turnover with depth is exponential, with ${\mathit{\xi }}_{\text{resp}}$ an empirical parameter (in m${}^{-\mathrm{1}}$) that controls the magnitude of the reduction . The larger the value of ${\mathit{\xi }}_{\text{resp}}$, the more inhibited the respiration is with increasing depth. Here ${\mathit{\xi }}_{\text{resp}}$ was tuned to give a realistic estimate of soil C in a vertically resolved version of JULES-C as in . When N is limiting, the respiration of the DPM and RPM pools are reduced by a factor of $F_{\text{N}}\left(z\right)$, which is also now a function of depth and dependent on the available N in the relevant layer. $M_{\text{p}}$ and $I_p$ are also calculated as a function of depth based on their relationship with respiration.

The vertical mixing of each soil N pool follows that of the soil C pools. 45$$\begin{array}{rl}{\displaystyle \frac{\partial N_{\text{DPM}}\left(z\right)}{\partial t}}& ={\displaystyle \frac{\partial }{\partial z}}\left(D\left(z\right){\displaystyle \frac{\partial N_{\text{DPM}}\left(z\right)}{\partial z}}\right)\\ & +\sum _i\left(v_i{f}_{\text{DPM},i}{\mathrm{\Lambda }}_{\text{n},i}{f}_{\text{lit}}\left(z\right)\right)-M_{\text{DPM}}\left(z\right)\end{array}$$ 46$$\begin{array}{rl}{\displaystyle \frac{\partial N_{\text{RPM}}\left(z\right)}{\partial t}}& ={\displaystyle \frac{\partial }{\partial z}}\left(D\left(z\right){\displaystyle \frac{\partial N_{\text{RPM}}\left(z\right)}{\partial z}}\right)\\ & +\sum _i\left(v_i(\mathrm{1}-f_{\text{DPM},i}){\mathrm{\Lambda }}_{\text{n},i}{f}_{\text{lit}}\left(z\right)\right)-M_{\text{RPM}}\left(z\right)\end{array}$$ 47$${\displaystyle \frac{\partial N_{\text{BIO}}\left(z\right)}{\partial t}}={\displaystyle \frac{\partial }{\partial z}}\left(D\left(z\right){\displaystyle \frac{\partial N_{\text{BIO}}\left(z\right)}{\partial z}}\right)+\mathrm{0.46}{I}_{\text{tot}}\left(z\right)-M_{\text{BIO}}\left(z\right)$$ 48$$\begin{array}{rl}{\displaystyle \frac{\partial N_{\text{HUM}}\left(z\right)}{\partial t}}& ={\displaystyle \frac{\partial }{\partial z}}\left(D\left(z\right){\displaystyle \frac{\partial N_{\text{HUM}}\left(z\right)}{\partial z}}\right)\\ & +\mathrm{0.54}{I}_{\text{tot}}\left(z\right)-M_{\text{HUM}}\left(z\right)\end{array}$$ $I_{\text{tot}}\left(z\right)$ is the total immobilisation in kg N m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$ in each layer (following Eq. ). $D\left(z\right)$ is the diffusivity in m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}}$ and varies both spatially and with depth ($z$) . 49$$D\left(z\right)=\left\{\begin{array}{ccc}{D}_{\text{o}}& ;& z\le \mathrm{1}\text{m}\\ \frac{D_{\text{o}}}{\mathrm{2}}(\mathrm{3}-z)& ;& \mathrm{1}\text{m}<z<\mathrm{3}\text{m}\\ \mathrm{0.0}& ;& z\ge \mathrm{3}\text{m}\end{array}\right\}$$ $D_{\text{o}}$ (m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}}$) varies spatially depending on the freeze and thaw state of the soil. In regions without permafrost, $D_{\text{o}}$ represents a bioturbation mixing rate equivalent to 1 cm${}^{\mathrm{2}}$ yr${}^{-\mathrm{1}}$. When permafrost is present, $D_{\text{o}}$ represents cryoturbation and increases to a value equivalent to 5 cm${}^{\mathrm{2}}$ yr${}^{-\mathrm{1}}$. The parameterisation of $D\left(z\right)$ in Eq. () means that the soil organic pools can transfer between the active layer and the permanently frozen soils in a steady-state climate, albeit at a relatively slow rate. The PFT-dependent litter inputs ($f_{\text{lit}}\left(z\right){\mathrm{\Lambda }}_{\text{n},i}$ for litter N) are distributed so that they decline exponentially with increasing depth. Here $f_{\text{lit}}\left(z\right)$ is independent of the PFT type and hence the root distribution is as follows: 50$$f_{\text{lit}}\left(z\right)={\displaystyle \frac{e^{-{\mathit{\xi }}_{\text{lit}}z}}{{\sum }_{z=\mathrm{0}}^{z_{\text{max}}}{e}^{-{\mathit{\xi }}_{\text{lit}}z}\text{d}z}},$$ where ${\mathit{\xi }}_{\text{lit}}$ is the parameter to reduce the litter input with increasing depth and is set to 0.2 m or 5 m${}^{-\mathrm{1}}$ and $z$ is the mid-point of each layer.

The mineralised gas emissions are now a function of depth ($N_{\text{gas}}\left(z\right)$) and are calculated by repeating Eq. () for each soil layer. Similarly, the litter decomposition rate modifier ($F_{\text{N}}$) is calculated by repeating a slightly modified version of Eq. () for each soil layer. In the vertically resolved version of Eq. (), if the soil layer is frozen $N_{\text{in}}$ is not available so effectively zero.

3.3 Inorganic nitrogen

The changes in the inorganic N pool result from deposition, fixation, immobilisation losses, mineralisation inputs, grid box mean plant uptake, and inorganic N losses through leaching and gaseous emission. For the bulk layer case (JULES-CN), these terms are simply added together:

51$$\begin{array}{rl}{\displaystyle \frac{\text{d}{N}_{\text{in}}}{\text{d}t}}& =N_{\text{dep}}+\sum _i{v}_i{\text{BNF}}_i-\sum _i{v}_i{\mathrm{\Phi }}_i\\ & +M_{\text{net}}-N_{\text{leach}}-N_{\text{gasI}},\end{array}$$ where $N_{\text{in}}$ is the inorganic N in kg N m${}^{-\mathrm{2}}$, $N_{\text{dep}}$ is prescribed N deposition in kg N m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$ and $v_i$ the fractional cover of each PFT $i$. The biological N fixation (${\text{BNF}}_i$) for each PFT $i$ is described in Sect.  below, and plant uptake (${\mathrm{\Phi }}_i$) for each PFT $i$ is described in Sect. . $M_{\text{net}}$ is the net mineralisation, which is the difference between $M_{\text{tot}}$ (Eq. ) and $I_{\text{tot}}$ (Eq. ) reduced by $N_{\text{gas}}$ (Eq. ). The loss of N from the system via the inorganic pool is the sum of leaching ($N_{\text{leach}}$ in kg N m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$) plus an additional gas loss to the atmosphere ($N_{\text{gasI}}$ in kg N m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$): 52$$N_{\text{gasI}}={\mathit{\gamma }}_{\text{n}}{N}_{\text{in}},$$ where ${\mathit{\gamma }}_{\text{n}}$ is a tunable parameter (in s${}^{-\mathrm{1}}$). The total N gas loss is the sum of $N_{\text{gasI}}$ above and $N_{\text{gas}}$ from Eq. (), with $N_{\text{gasI}}$ representing approximately 90 % of the total gas loss. This additional gas loss term ($N_{\text{gasI}}$) represents missing processes relating to the gaseous loss of inorganic N and limits the effective mineral N pool size. Including $N_{\text{gasI}}$ ensures that available N does not increase excessively, potentially due to excessive biological N fixation in regions where the NPP is very close or equal to the ${\text{NPP}}_{\text{pot}}$. In the current model configuration ${\mathit{\gamma }}_{\text{n}}$ is set to 0.0028 d${}^{-\mathrm{1}}$ such that the whole pool turns over once every model year.

The leaching of N ($N_{\text{leach}}$ in kg N m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$) 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 ($\mathit{\alpha }$, dimensionless) representing the effective solubility of N. $\mathit{\alpha }$ is assumed to have a value of 0.1 and in JULES-CN represents the combined sorption of all inorganic N species . 53$$N_{\text{leach}}=\mathit{\alpha }(N_{\text{in}}/{\mathit{\theta }}_{\mathrm{1}\text{m}})Q_{\text{subs}},$$ where ${\mathit{\theta }}_{\mathrm{1}\text{m}}$ is the soil water content in the top 1 m of soil in kg m${}^{-\mathrm{2}}$ (so the inorganic N is assumed to occupy the top 1 m of soil) and $Q_{\text{subs}}$ is the total subsurface runoff in kg m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$.

3.3.1 Biological nitrogen fixation (BNF)

Biological nitrogen fixation (BNF) is the largest natural supplier of N to the terrestrial ecosystem. Following the secondary model of , N fixation is determined as a linear proportion of the net primary production before N limitation of each PFT $i$ (${\text{NPP}}_{\text{pot},i}$):

54$${\text{BNF}}_i=\mathit{\zeta }{\text{NPP}}_{\text{pot},i}.$$ ${\text{NPP}}_{\text{pot},i}$ is defined in the same way as the net primary productivity in JULES before the explicit N cycle was included, i.e. before the excess carbon ($\mathrm{\Psi }$) is removed. BNF as a function of NPP is an established method used and assessed in other models . While some models utilise more complex BNF representations , a lightweight approach is preferred here while the benefits of extra computational expense on BNF are not yet established, and evidence is lacking that a different simple representation (e.g. evapotranspiration) would perform better . However, changes in NPP may not accurately reflect changes in BNF with forcings such as elevated atmospheric carbon dioxide or additional N .

The rate of fixation ($\mathit{\zeta }$) is set such that global present-day net primary productivity of approximately 60 Pg C yr${}^{-\mathrm{1}}$ results in approximately 100 Tg N yr${}^{-\mathrm{1}}$ fixation (0.0016 kg N kg C${}^{-\mathrm{1}}$), within the range of recent global observation-based estimates of BNF . The parameterisation based on NPP results in a latitudinal gradient with the highest rates of fixation in the tropics and lowest in boreal forests and arctic tundra, which is consistent with some estimates of BNF though not recent observational meta-analyses .

In JULES-CN, which has a bulk soil biogeochemistry parameterisation, the BNF is directly transferred into the single inorganic soil N pool and becomes available as inorganic N. However, in JULES-CN${}_{\text{layer}}$ the BNF is distributed vertically in the soil depending on the fraction of roots in each layer. If a soil layer is frozen, there is no BNF into that layer. If the whole soil is frozen, fixed N goes into the inorganic N pool in the top layer.

3.3.2 Vertical discretisation of inorganic nitrogen

In JULES-CN${}_{\text{layer}}$ there is an inorganic N pool in each soil layer. The dynamics are very similar to Eq. (), but most of the components now vary with depth:

55$$\begin{array}{rl}{\displaystyle \frac{\text{d}{N}_{\text{in}}\left(z\right)}{\text{d}t}}& =\sum _i{v}_i{\text{BNF}}_i{f}_{R,i}\left(z\right)-\sum _i{v}_i{\mathrm{\Phi }}_i{f}_{I,i}\left(z\right)\\ & +M_{\text{net}}\left(z\right)-N_{\text{flux}}\left(z\right)-N_{\text{gasI}}\left(z\right).\end{array}$$ Any N deposition ($N_{\text{dep}}$) is added to the top layer of the soil only. The modifications to each term to ensure they vary appropriately with depth are discussed below. The additional parameters in Eq. () are $f_{R,i}\left(z\right)$ – the fraction of roots in each layer for PFT $i$ (Eq. ), $f_{I,i}\left(z\right)$ – the fraction of available inorganic N in each layer for PFT $i$ (Eq. ), and $N_{\text{flux}}\left(z\right)$ – the transport of inorganic N through the layer by the soil water fluxes (Eq. ).

As in Eq. (), the net mineralisation flux ($M_{\text{net}}\left(z\right)$) is the difference between $M_{\text{tot}}\left(z\right)$ and $I_{\text{tot}}\left(z\right)$ reduced by $N_{\text{gas}}\left(z\right)$ for each layer (see Sect. ). Inputs from biological N fixation from PFT $i$ are distributed according to the root profile of the PFT under consideration ($f_{R,i}\left(z\right)$): 56$$f_{R,i}\left(z\right)={\displaystyle \frac{f_{\text{root},i}\left(z\right)}{{\sum }_{z=\mathrm{0}}^{z_{\text{max}}}{f}_{\text{root},i}\left(z\right)\text{d}z}},$$ where $f_{\text{root},i}\left(z\right)$ is the volumetric root fraction of PFT $i$ at a given soil level and $z_{\text{max}}$ is the maximum depth of the soil in metres. Gas loss from the inorganic N ($N_{\text{gasI}}\left(z\right)$) occurs in each layer but with an additional exponential decay term that is a function of depth (similar to that used in Eq.  for the soil decomposition). This term empirically represents the factors that reduce soil activity with depth. The additional gas loss term thus becomes 57$$N_{\text{gasI}}\left(z\right)={\mathit{\gamma }}_{\text{n}}{N}_{\text{in}}\left(z\right)e^{-{\mathit{\xi }}_{\text{resp}}z}.$$ This leaves two terms in Eq. (): the plant uptake term ($\sum _i{v}_i{\mathrm{\Phi }}_i{f}_{I,i}\left(z\right)$), which is PFT-dependent, and the $N_{\text{flux}}\left(z\right)$ term, which replaces the leaching term from Eq. (). These have a more process-based representation in the layered case. When calculating the plant uptake term we assume that plants cannot access all the inorganic N. Firstly, if a soil layer is frozen then plants cannot uptake any of the N 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_{\text{root},i}\left(z\right)$ (which reduces with depth). So for each PFT, $i$, the available amount of the inorganic N pool ($N_{\text{avail},i}\left(z\right)$ in kg N m${}^{-\mathrm{2}}$) at equilibrium that could potentially be extracted by the vegetation is given by 58$$N_{\text{avail},i}\left(z\right)=f_{\text{root},i}\left(z\right)N_{\text{in}}\left(z\right)T\left(z\right),$$ where $T\left(z\right)$ is zero when the soil temperature is 0 ${}^{\circ }$C or colder and 1 when it is above 0 ${}^{\circ }$C. However, the system is not necessarily in equilibrium – as N is taken up from this available pool around the roots, there will be a delay in this volume getting “re-filled”. We assume that the inorganic N 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 of inorganic N by the plants on each TRIFFID time step we additionally update the available N pool according to 59$${\displaystyle \frac{N_{\text{avail},i}\left(z\right)}{\text{d}t}}={\mathit{\gamma }}_{\text{dif}}\left(f_{\text{root},i}\right(z\left)N_{\text{in}}\right(z)-N_{\text{avail},i}(z\left)\right),$$ where ${\mathit{\gamma }}_{\text{dif}}$ is the rate of diffusion back to the equilibrium, set by default to 0.28 d${}^{-\mathrm{1}}$ or approximately 100 yr${}^{-\mathrm{1}}$. $N_{\text{avail},i}\left(z\right)$ is then multiplied by $T\left(z\right)$ to incorporate the frozen soil effect. Any biological N fixation goes directly into the available pool. Plant uptake is extracted entirely from the available N pool, and the dependence on depth is according to the same profile as the available N, i.e. 60$$f_{I,i}\left(z\right)={\displaystyle \frac{N_{\text{avail},i}\left(z\right)}{{\sum }_{z=\mathrm{0}}^{z_{\text{max}}}{N}_{\text{avail},i}\left(z\right)\text{d}z}}.$$ All of the other fluxes are simply added in such a manner so as to maintain the ratio of the available to total inorganic N pools that would be present if the available and total pools were in equilibrium. Therefore the only two processes which take the available and total pools out of equilibrium are biological N fixation and uptake.

Leaching is now done in a process-based manner, where the inorganic N is transported through the soil profile by the soil water fluxes. For any given soil layer of thickness $\mathit{\delta }z$, the inorganic N flux ($N_{\text{flux}}$) is given by 61$$N_{\text{flux}}\left(z\right)=\mathit{\alpha }\mathit{\delta }z{\displaystyle \frac{\text{d}}{\text{d}z}}\left(W_{\text{flux}}\left(z\right){\displaystyle \frac{N_{\text{in}}\left(z\right)}{\mathit{\theta }\left(z\right)}}\right),$$ where $\mathit{\theta }\left(z\right)$ is the soil water content of the layer in kg m${}^{-\mathrm{2}}$ and $W_{\text{flux}}\left(z\right)$ is the flow rate of the water through the layer in kg m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$. Multiplying by $\mathit{\delta }z$ gives the change in N content for each layer. The total leaching is then the sum of all N that leaves the soil both laterally from each layer or from the bottom of the soil profile.

Table  summarises the extra parameters required for the soil biogeochemistry component of JULES-CN and JULES-CN${}_{\text{layer}}$ alongside their values.

Table 2

A summary of the extra parameters required for the soil biogeochemistry component of JULES-CN and JULES-CN${}_{\text{layer}}$.

Variable	Value	Description	Equation
Bulk soil nitrogen
$\mathit{\zeta }$	0.0016 kg N kg C ${}^{-\mathrm{1}}$	Rate of BNF	Eq. ()
${\text{CN}}_{\text{soil}}$	10 kg C kg N${}^{-\mathrm{1}}$	CN ratio of BIO and HUM pools	Eq. ()
$f_{\text{gas}}$	0.01 (proportion)	Fraction of net mineralisation emitted as gas to atmosphere	Eq. ()
${\mathit{\gamma }}_{\text{n}}$	$\mathrm{3.215}\times {\mathrm{10}}^{-\mathrm{8}}$ s${}^{-\mathrm{1}}$	Imposed turnover coefficient to determine $N_{\text{gasI}}$ release from $N_{\text{in}}$	Eq. ()
$\mathit{\alpha }$	0.1 (proportion)	Effective solubility of nitrogen in water	Eq. ()
Vertically resolved soil carbon
${\mathit{\xi }}_{\text{resp}}$	0.8 m${}^{-\mathrm{1}}$	Parameter to control reduction of respiration with depth	Eq. ()
$D_{\text{o}}$	bioturbation – 0.001 m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}}$	Soil carbon and nitrogen mixing rate	Eq. ()
	cryoturbation – 0.005 m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}}$
${\mathit{\xi }}_{\text{lit}}$	5 m${}^{-\mathrm{1}}$	Parameter to control reduction of litter input with depth	Eq. ()
Vertically resolved soil carbon and nitrogen
${\mathit{\gamma }}_{\text{dif}}$	100 per 360 d	Rate of diffusion transferring the inorganic nitrogen from $N_{\text{in}}$ to $N_{\text{avail}}$	Eq. ()

4 Historical simulations

Global transient simulations were carried out following the protocol for the S2 experiments in TRENDY , which include time varying climate, CO${}_{\mathrm{2}}$, and N deposition but pre-industrial land use. Time-varying CO${}_{\mathrm{2}}$ and climate were from the from the CRU-NCEP dataset (v7, 10.5065/PZ8F-F017, ). The fraction of agriculture in each grid cell was set to the pre-industrial value defined by . N deposition was taken from a ACCMIP multi-model dataset interpolated to annual fields . The model resolution was N96 (1.25${}^{\circ }$ latitude $\times$ 1.875${}^{\circ }$ longitude).

Results from three different JULES model configurations are presented here.

• JULES-CN includes the newly developed soil and vegetation coupled C and N cycle.

• JULES-C is shown for comparison purposes and represents the soil and vegetation C cycle as used in .

• JULES-CN${}_{\text{layer}}$ is a version of JULES-CN that has identical above-ground processes to JULES-CN but additionally includes vertically discretised soil biochemistry.

In each case five PFTs were used: broadleaf trees, needleleaf trees, C${}_{\mathrm{3}}$ and C${}_{\mathrm{4}}$ grass, and shrubs. Plant competition was allowed, with TRIFFID updating vegetation fractions on a 10 d time step. These three configurations of JULES adopt the standard four layer soils with a maximum depth of 3 m. However, it should be noted that and adopt a configuration that increases both the maximum soil depth and number of soil layers – a configuration that is recommended for detailed scientific study of northern high latitudes. The sole difference between JULES-C and JULES-CN is the inclusion of the N cycle. JULES-CN${}_{\text{layer}}$ additionally has vertically discretised soil biogeochemistry. There are no differences in any of the shared model parameters that were initially tuned for the JULES-C configuration. This enables a direct comparison between the different configurations.

The simulations were initialised using pre-industrial conditions. The models were spun up by using the meteorological data for the period 1860–1879 repeatedly until the change in the carbon stocks was less than 0.01 % decade${}^{-\mathrm{1}}$ globally. The soil C distribution in JULES-CN${}_{\text{layer}}$ is particularly slow to reach equilibrium. Therefore the “modified accelerated decomposition” technique (modified AD) described by was used to spin the soil C in these versions up to an initial pre-industrial equilibrium distribution . Further spin-up was then carried out for these layered models using repeated pre-industrial conditions until the change in soil C was again less than 0.01 % decade${}^{-\mathrm{1}}$ globally. It should be noted that neither transient land use change or fertiliser were included in any of these simulations.

5 Results

This paper mainly focuses on the differences in JULES output when including the N cycle in the model configuration. When available, we additionally use any observation-based estimates to evaluate the quality of the simulations. First, a broad-brush comparison between JULES-CN, JULES-C, and JULES-CN${}_{\text{layer}}$ is made. This is followed by a more complete discussion of the impact of the N cycle on the carbon stocks and fluxes and their changes over time. Then we show spatial distributions and time series of the N stocks and fluxes. Finally, the extra processes modelled in JULES-CN${}_{\text{layer}}$ are assessed. For completeness, figures often include both JULES-CN and JULES-CN${}_{\text{layer}}$, but JULES-CN${}_{\text{layer}}$ is only discussed at the end of the results.

Figure 3

Total area covered by each vegetation type for the three different JULES configurations. The observations are derived from the European Space Agency (ESA) Climate Change Initiative (CCI) Land Cover data for 2010 converted to JULES PFTs by .

[Figure omitted. See PDF]

Figure 4

Carbon and nitrogen stocks and fluxes for JULES-CN, JULES-C, and JULES-CN${}_{\text{layer}}$ for the period 1996–2005 (after ). C is carbon, N is nitrogen; $r_{\text{h}}$ is heterotrophic respiration, $r_{\text{a}}$ is autotrophic respiration, TER is total ecosystem respiration; GPP is gross primary productivity, NPP is net primary productivity, SOM is soil organic matter, BNF is biological N fixation, and N gas is the sum of $N_{\text{gas}}$ and $N_{\text{gasI}}$, with $N_{\text{gasI}}$ representing approximately 90 % of the total gas loss. The black numbers are the observation-constrained values from the literature, where observation-based values are not available JULES is compared with other global models: (a) heterotrophic respiration: ; (b) TER: ; (c) TER: ; (d) GPP: ; (e) vegetation carbon and SOM + litter carbon: ; (f) BNF ; (g) BNF ; (h) Vegetation nitrogen: ; (i) NPP: ; (j) soil organic nitrogen: ; (k) soil organic nitrogen: ; (l) nitrogen losses including nitrogen leaching: ; (m) nitrogen leaching: ; (n) nitrogen leaching: ; (o*) organic nitrogen immobilisation and mineralisation and plant uptake ; (p*) nitrogen uptake, vegetation nitrogen, and nitrogen emissions ; (q*) nitrogen uptake and inorganic nitrogen content ; and (r*) nitrogen uptake and total nitrogen emissions . (o*), (p*), (q*), and (r*) are model-derived estimates.

[Figure omitted. See PDF]

It should be noted that the addition of the N cycle in JULES is only one component of the recent developments. In the UKESM configuration of JULES the N cycle was combined with a new competition scheme and additional PFTs, both of which modify the global vegetation distribution. Therefore, we are most interested in the changes in the vegetation distribution between the different versions which will be caused by the N cycle. Figure  shows the total area covered by each type of vegetation. The Climate Change Initiative (CCI) land cover observations are added for completeness. In general, the models all tend to overestimate the shrubs and underestimate the grass. However, shows that once the additional PFTs and new competition scheme are included the model does a good job of representing the vegetation distribution.

As expected, the configurations with the N cycle have more bare soil and less vegetation than JULES-C. This is mainly observed as a decrease in the shrub and grass regions in JULES-CN. As we shall see later (Fig. ), this is because the tropical forests dominate the tree region and their growth is not limited by N in the model. JULES-CN${}_{\text{layer}}$ has a reduction in trees compared to JULES-CN, which is focused in the boreal region where it is more likely to simulate grass or shrubs.

5.1 Summary of C and N stocks and fluxes

Figure  provides an overview of the stocks and fluxes of C and N in JULES-CN and JULES-CN${}_{\text{layer}}$ and compares them with JULES-C. As expected for a present-day simulation, the majority of C stocks and fluxes are very similar for JULES-C and JULES-CN. The main difference is the present-day NPP, which is $\sim$ 12 % higher in JULES-C than in JULES-CN. There is also a small reduction in the GPP of $\sim$ 4 % caused by some differences in the vegetation fractional cover distribution (Fig. ) and indirectly resulting from the N limitation.

Figure 5

Zonal total values of (a) net primary productivity (NPP) and (b) gross primary productivity (GPP) for JULES-C, JULES-CN, and JULES-CN${}_{\text{layer}}$ simulations for the period 1996–2005 in Pg C (degree latitude)${}^{-\mathrm{1}}$ yr${}^{-\mathrm{1}}$. The observational constraint for NPP is from MODIS and that for GPP is from . The zonal mean carbon use efficiency (CUE $=$ NPP $/$ GPP) is shown in (c). The CUE observational constraint was digitised from . Also shown are changes in (d) NPP, (e) GPP, and (f) CUE over the historical period with respect to the multi-annual mean period of 1860–1899.

[Figure omitted. See PDF]

Soil organic N and vegetation N are both consistent with the available observation-based estimates of stocks. The biological N fixation is tuned to be approximately 100 Tg N ${\mathrm{yr}}^{-\mathrm{1}}$ in the present day and the N deposition is prescribed. The majority of N losses from the land surface occur via the gaseous pathway with total losses of 111 Tg N ${\mathrm{yr}}^{-\mathrm{1}}$ for JULES-CN. Leaching is fairly low at 7 Tg N ${\mathrm{yr}}^{-\mathrm{1}}$ compared to estimates of leaching, which are as high as $\sim$ 25 %–55 % of N inputs in European forests and range between 59 and 118 Tg N ${\mathrm{yr}}^{-\mathrm{1}}$ in the available observations . There is no N fertiliser applied in the model, which might partially explain why the leaching is so low. In reality there is $\sim$ 200 Tg N applied annually as either manure or fertiliser , a proportion of this will be leached, resulting in an increase of global leaching. N uptake and net N mineralisation are relatively high and are fairly comparable in magnitude implying a largely closed cycling of nutrients between vegetation and soil. These N stocks and fluxes are also consistent with results from other models such as , , , and .

5.2 Comparing C stocks and fluxes

Carbon use efficiency (CUE) is defined as the ratio of net C gain to gross C assimilation during a given period (NPP $/$ GPP). Plants with a higher CUE have a lower autotrophic respiration and allocate more C from photosynthesis to the terrestrial biomass and vice-versa. In JULES-CN there is less C available to be allocated because it is constrained by the amount of N present. This reduces the C use efficiency. Figure  shows the zonal total GPP and NPP for JULES-CN and JULES-C. As expected from Fig. , the NPP and GPP have very similar latitudinal profiles for the two model configurations. Both JULES-C and JULES-CN have a higher GPP in the tropics than the observations, but they are more comparable in the extra-tropical latitudes where the GPP tends to be smaller. The NPP in JULES-CN is less than JULES-C and generally closer to the MODIS observations particularly in the tropics. Figure  also shows the zonal mean CUE. JULES-CN has a lower CUE than JULES-C for all latitudes. On average it is 0.44 for JULES-CN and 0.49 for JULES-C. JULES-CN is consistently low compared to the observation-based dataset with a bias of $\sim$ 0.09. This bias is relatively constant with latitude. However, considerable uncertainties remain in these estimates.

Figure also shows the changes in these C fluxes for the period 1860–2007 with respect to the multi-annual mean period of 1860–1899. Changes over time are shown to enable the differences between the two different model configurations to be more easily compared. Apparently small differences between JULES-C and JULES-CN in the NPP and GPP become more noticeable in the CUE. The small differences between JULES-C and JULES-CN in GPP are mainly caused by small changes in the vegetation distribution and a slight increase in bare soil in JULES-CN. In the case of NPP, JULES-C increases quicker than JULES-CN because JULES-CN becomes progressively more N-limited. The change in CUE shows the impact of the N cycle on the uptake of C by the vegetation in JULES-CN over the 20th century. There is an increase in CUE in both configurations, mainly caused by CO${}_{\mathrm{2}}$ fertilisation, but this is limited by N in the JULES-CN configuration.

The zonal total soil and vegetation C stocks are shown in Fig. . The vegetation C is very similar for both JULES-C and JULES-CN, as expected from Fig. , and is consistent with the available observations. There are some differences in the soil C in the northern high latitudes, with JULES-CN having slightly less soil C than JULES-C. This is a consequence of the higher N limitation on JULES-CN leading to less litter fall and subsequently less soil C. The corresponding N-limitation-induced reduction in soil organic matter decomposition is not strong enough to offset the decrease in C input leading to a smaller pool size.

Figure 6

Zonal total values of (a) vegetation and (b) soil C for JULES-C, JULES-CN, and JULES-CN${}_{\text{layer}}$ simulations for the period 1996–2005 in Pg C (degree latitude)${}^{-\mathrm{1}}$. For the vegetation C the observation-based constraints are Saatchi: ; GEOCARB: ; and biomass: . The observation-based constraints for the soil carbon are IGBP-DIS: ; WISE: ; and Carvahlais: .

[Figure omitted. See PDF]

5.2.1 Net ecosystem exchange

A key measure of a land C cycle model is how well it simulates the temporal variation of the land C sink, which is the difference between net ecosystem exchange (NEE) and the flux of C 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 . As a result, simulation of NEE variability is highly relevant to climate–carbon cycle projections .

Figure 7

Evaluation of global annual mean NEE from JULES-CN, JULES-C, and JULES-CN${}_{\text{layer}}$ compared with observations based on estimates from Global Carbon Project (GCP) over the period from 1960 to 2009 inclusive. Positive values represent the land surface as a net sink of carbon. The solid lines are the data and the dashed–dotted lines represent a linear fit of the data against time.

[Figure omitted. See PDF]

Figure compares global annual mean values of net ecosystem exchange (NEE; defined as NPP minus heterotrophic respiration) for JULES-C and JULES-CN to observation-based estimates from the Global Carbon Project. We specifically focus on the years from 1960 to 2009, which is the maximum overlap period between the model simulations and the GCP annual budget data . To avoid the circularity of using GCP estimates of NEE, which are themselves derived from land surface models, we instead calculate the GCP estimates of NEE as the residual of the best estimates of the total emissions from fossil fuel (FF) plus land use change (LU) and the rate of increase of the carbon content of the atmosphere ($F_{\text{a}}$) plus the ocean ($F_{\text{o}}$): 62$${\text{NEE}}_{\text{gcp}}=\text{FF}+\text{LU}-F_{\text{a}}-F_{\text{o}}.$$

Figure 8

Biome-based (a) ecosystem turnover times and (b) soil carbon turnover times calculated on a grid cell by grid cell basis then aggregated temporally to biome level. JULES-C, JULES-CN, and JULES-CN${}_{\text{layer}}$ are shown for the period 1996–2005. The land surface is split into biomes based on the 14 World Wildlife Fund terrestrial ecoregions and characterised by . The observed ecosystem residence times are derived from the global dataset and the observed soil residence times are from the WISE soil carbon combined with the soil respiration.

[Figure omitted. See PDF]

The observations and both of the models show an upward trend in NEE but with very significant interannual variability (Fig. ). Due to N limitations on CO${}_{\mathrm{2}}$ fertilisation, mean NEE in JULES-CN (1.66 Pg C ${\mathrm{yr}}^{-\mathrm{1}}$) is lower than in JULES-C (2.06 Pg C ${\mathrm{yr}}^{-\mathrm{1}}$) and also lower than the estimate from GCP (2.11 Pg C ${\mathrm{yr}}^{-\mathrm{1}}$). This absolute value will be sensitive to the vegetation cover, which is much improved by including the height-based competition as has been done in UKESM1 . 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 ). The correlation coefficient for NEE between the JULES-CN and GCP estimates ($r=\mathrm{0.71}$) is also improved compared to JULES-C ($r=\mathrm{0.63}$). There remains a significant underestimate of NEE in the years following the Pinatubo volcanic eruption in 1991, most likely due to the neglect of diffuse-radiation fertilisation in these versions of JULES . However, it is especially notable that JULES-CN significantly reduces the systematic overestimate of NEE seen in JULES-C during extended La Niña periods, such as the years centred around 1974 and 2000 (Fig. ).

Table 3

Statistics of NEE from JULES-CN, JULES-C, JULES-CN${}_{\text{layer}}$, and the GCP observation-based estimates , over the period from 1960 to 2009 inclusive. Columns 2–4 respectively show the mean, linear trend, and the interannual variability (IAV; standard deviation) around that trend. Column 5 shows the correlation coefficient between each model NEE time series and the GCP time series.

	Mean	Trend	IAV
	(Pg C	(Pg C ${\mathrm{yr}}^{-\mathrm{1}}$	(Pg C
	${\mathrm{yr}}^{-\mathrm{1}}$)	${\mathrm{yr}}^{-\mathrm{1}}$)	${\mathrm{yr}}^{-\mathrm{1}}$)	$r$
JULES-CN	1.66	0.025	0.86	0.71
JULES-C	2.06	0.034	1.31	0.63
JULES-CN${}_{\text{layer}}$	1.75	0.026	0.83	0.64
GCP(residual)	2.11	0.027	1.01

5.2.2 Residence times

In general, carbon residence times of the soil and ecosystem are given by the stocks divided by the fluxes. These are emergent properties of the model and thus a valuable metric to evaluate. Figure  shows the ecosystem residence time and the soil C residence times for different biomes. Here, the land surface is split into biomes based on the 14 World Wildlife Fund terrestrial ecoregions and characterised by . The ecosystem residence time defined as the total ecosystem C divided by the GPP is shown in Fig. a. These residence times have been estimated from a multi-annual mean on a grid cell by grid cell basis and then aggregated to biomes. The observational values were derived in a similar way using spatial data from . In general the ecosystem residence times are slightly reduced in JULES-CN compared with JULES-C, both of which are generally lower than in the observations. The largest difference between observed and modelled ecosystem residence time occurs in the tundra and boreal regions and the grasslands where the observed residence times are much longer than either JULES-C or JULES-CN. The soil carbon residence time is shorter than the observation-based measure in the tundra and the boreal regions but longer in the grassland regions. Overall, this leads to the global soil carbon residence time in the model being too short. When vertical discretisation, including additional permafrost processes, is added in JULES-CN${}_{\text{layer}}$, the residence times in the boreal and tundra increase notably (see Sect.  for further discussion).

Figure 9

The spatial distribution of the response ratio, defined as the potential amount of carbon that can be allocated to growth and spreading of the vegetation (${\text{NPP}}_{\text{pot}}$), as a fraction of the NPP achieved in the natural state for (a) JULES-CN and (b) JULES-CN${}_{\text{layer}}$. A value greater than 1 means that the addition of nitrogen will enhance NPP. Any grid cells with an annual NPP of less than 0.016 g C m${}^{-\mathrm{2}}$ are set to missing. This is the spatial distribution of the metric shown in Fig. . Panel (c) shows the change in the response ratio over the historical period with respect to the multi-annual mean from the period of 1860–1899.

[Figure omitted. See PDF]

Figure 10

The response ratio is the ratio of the potential amount of carbon that can be allocated to growth and spreading of the vegetation (${\text{NPP}}_{\text{pot}}$) compared with the actual amount achieved in the natural state (NPP). As in Fig. , any grid cells with an annual NPP of less than 0.016 g C m${}^{-\mathrm{2}}$ are set to missing. The median of JULES-CN and JULES-CN${}_{\text{layer}}$ are shown for each biome for the period 1996–2005. The biomes are discussed in more detail in Fig. . The observational constraint is taken from Table 1 in , which summarises a meta-analysis of nitrogen addition experiments. The black bars show the mean of the observations, and the red lines show the uncertainty.

[Figure omitted. See PDF]

5.3 Impact of N limitation

IN JULES-CN and JULES-CN${}_{\text{layer}}$ the N limitation mainly acts through reducing the NPP. This can be quantified using the response ratio, which is defined as the ratio of the potential amount of C that can be allocated to growth and spreading of the vegetation (${\text{NPP}}_{\text{pot}}$) compared with the actual amount achieved in the natural state (NPP). Both of these diagnostics are output from the JULES simulations. Figure  shows the spatial distribution of the model-simulated response ratio. Green areas are not very N-limited, with a response ratio close to 1.0, and yellow areas are more N-limited, with a larger response ratio. There are distinct regions of N limitation – in Australia and South Africa, the Sahel, western Europe, and parts of Siberia. However much of the global land surface, and particularly forested regions, has relatively weak N limitation. Figure c also shows the JULES-CN response ratio has obvious inter-annual variability superimposed on an increasing trend over the 20th century, indicating increasing N limitation, which will limit the increase in carbon use efficiency shown in Fig. f.

Figure shows the biome-based response ratio of net primary productivity. All biomes have a response ratio of greater than 1 in both the model and observations, which means that adding extra N to the system will enhance the NPP achieved. Globally the response ratio is lower than the observations but for the majority of the biomes including the tropical forests and the tundra the model response ratios fall within the range of uncertainties of the observations. However, suggests the tropical forest is somewhat N-limited, whereas in JULES-CN tropical forest is not a N-limited biome. Phosphorus has long been considered the most limiting nutrient in tropical regions ; therefore, we expect JULES to simulate a larger response ratio in the future once a phosphorus cycle is added.

In the model, the soil C decomposition can be limited when the N available in the soil is less than the N required by decomposition. This process does not play a major role in our simulations.

5.4 Nitrogen stocks and fluxes

The zonal profile of soil organic nitrogen (Fig. ) shows a similar distribution to the soil organic C (Fig. ), reflecting the relatively consistent C to N ratio of the soil within the model. CN${}_{\text{soil}}$ – the C to N ratio of the HUM and BIO pools – is a spatially constant parameter set to 10 in these simulations. The observed soil N content is slightly higher at all latitudes than simulated by JULES-CN, particularly in the northern tundra region. In contrast to the zonal distribution of soil organic nitrogen, the soil inorganic nitrogen in JULES-CN is larger in the tropics than in the northern high latitudes (Fig. ). Figure  shows the net soil N mineralisation fluxes are large in the tropics and smaller in the northern regions. This is reflected in the spatial distribution of the N uptake. As might be expected, the spatial distribution of the N uptake as a fraction of N demand is similar to the N limitation shown in Fig. . Biological N fixation and N gas losses are an order of magnitude smaller than the N uptake and net N mineralisation. However, again the spatial patterns are very comparable. N leaching is generally very small, except in parts of South America and South-East Asia. Figure  shows a slight increase in the N demand and N uptake over the 20th century associated with the increase in vegetation growth (Fig. ). Similarly, there is an increase in the BNF that is parameterised such that it is proportional to the NPP.

Figure 11

The zonal total soil organic and inorganic nitrogen stocks in Pg N (degree of latitude)${}^{-\mathrm{1}}$. JULES-CN${}_{\text{layer}}$ shows the stocks for the top 1 m of soil. The observations of nitrogen stocks are from .

[Figure omitted. See PDF]

Figure 12

Spatial distribution of N fluxes for JULES-CN for the period 1996–2005. JULES-CN${}_{\text{layer}}$ is not shown because the spatial patterns are very similar to those for JULES-CN.

[Figure omitted. See PDF]

5.5 Impact of vertical discretisation of soil biochemistry

This section discusses the differences between JULES-CN and JULES-CN${}_{\text{layer}}$. In general, over the tropics and southern latitudes, JULES-CN${}_{\text{layer}}$ is very comparable to JULES-CN. The majority of the differences occur in the northern regions where there is soil freezing, i.e. either permafrost or seasonally frozen soils. The reduction in global mean tree-covered area seen in Fig.  is caused by a reduction in the boreal regions, which have a larger proportion of shrubs and grasses in JULES-CN${}_{\text{layer}}$. In the higher latitudes the soil in JULES-CN${}_{\text{layer}}$ also has more organic C (Fig. ). This increase in soil organic C represents a store of permafrost carbon more comparable to the carbon found by and . This build-up of carbon in JULES-CN${}_{\text{layer}}$ occurs because the decomposition deeper in the soil is reduced with the lower soil temperatures at depth – the soil C in JULES-CN only responds to the soil temperatures near the surface which are warmer. This also causes an increase in the residence time of the soil carbon shown in Fig. b. The modelled soil C residence time in JULES-CN${}_{\text{layer}}$ is now much longer and more comparable to that observed.

The spatial distributions of N fluxes in JULES-CN${}_{\text{layer}}$ (not shown) are very similar to those of JULES-CN. In addition, the time series of changes in N fluxes over the 20th century are also comparable (Fig. ). The main differences are in the N gas loss, which is larger in JULES-CN${}_{\text{layer}}$, and the N leaching, which is larger in JULES-CN. Figure shows an increase in both organic and inorganic N in JULES-CN${}_{\text{layer}}$ over that in JULES-CN in the northern high latitudes similar to that seen in the organic C. As is the case for soil organic C, in the colder regions the soil N builds up within the frozen soil because of the limitation of the decomposition rates by cold temperatures; therefore, larger pools deeper in the soil are maintained in an equilibrium climate. The parameterisation of the vertically resolved soil biogeochemistry means that once JULES-CN${}_{\text{layer}}$ is spun-up there is inorganic N within the soil profile that cannot be taken up by the vegetation, either because the soil is frozen or because the roots cannot readily access it. This means that the extra inorganic N in JULES-CN${}_{\text{layer}}$ (Fig. ) is mainly stored deeper in the soil profile and within the permafrost itself and is typically unavailable in the current climate. This improved representation of the soil biogeochemistry will have implications for simulations of climate change feedbacks from the northern high latitudes.

Figure 13

N fluxes for JULES-CN and JULES-CN${}_{\text{layer}}$ over the historical period.

[Figure omitted. See PDF]

6 Discussion

This study presents the first implementation of nutrient cycles into the UK land and Earth system models. The scheme is parsimonious in that it captures the first-order and large-scale effects of interacting carbon and nitrogen on the land surface in the simplest way possible. One important assumption is that of fixed plant stoichiometry and that a plant strives to achieve stoichiometric homeostasis to maintain ecosystem structure, function, and stability under change environments . This assumption has some support in the literature e.g. and is a common approach amongst complex DGVMs . However, recent meta analyses of field observations show a distinct increase in foliar N to additional N availability and a modelling study found that assuming fixed C : N ratios and/or scaling leaf N concentration changes to other tissues, as employed here, were not supported by available evaluation data . Employing flexible stoichiometry has the potential to significantly affect the modelled biogeochemical feedbacks. For instance, nutrient limitation tends to limit productivity and thus the production of litter, which is the input to soil organic matter, leading to a reduction in soil carbon that the nutrient limitation in soil turnover is too weak to oppose. Allowing for flexible stoichiometry may lead to a lower litter quality but a comparable amount of litter. This reduction in litter quality will strengthen the soil turnover response, possibly leading to an overall increase in soil organic matter. Plant stoichiometric relationships are therefore a key uncertainty in assessing the carbon cycle feedbacks to climate change. Future versions of this model will explore the use of plant trait information to parameterise leaf, root, and wood C : N ratios for individual PFTs and further developments to allow for flexible stoichiometry.

While the total BNF in JULES-CN is in the range of and , the spatial distribution of BNF more heavily favours the tropics than recent observations suggest . The response of BNF to the multiple factors likely to occur in future varies between factor (warming, elevated atmospheric carbon dioxide, drought, N deposition, etc.), biome, and BNF type (nodulating, bryophyte, litter, etc.) . Therefore how BNF will change is spatially variable and not controlled by a single factor. A move from an empirical to a process-driven BNF function may provide better fit to present-day BNF distribution and more robust future projections.

Further work is required to explore the impact of a spatially varying soil C to N ratio, which can vary widely depending on the amount and decomposition of organic matter within the soil. For example, peat soils have relatively high C to N ratios up to 30–40 . This type of soil is not yet included within JULES. In addition, N leaching is very low in the model, notwithstanding the lack of N fertiliser. One reason for this could be that too much mineral N is assumed to be sorped within the soil. This requires further evaluation and potential modifications to the scheme.

In this paper we have not explicitly separated the impact of CO${}_{\mathrm{2}}$ fertilisation from climate change or from the impact of N deposition. However, this was explored by , who put the response of JULES in context by comparing it with the responses from four additional land surface models and a meta-analysis of site observations. used a slightly different configuration of JULES (JULES-ES), which is the configuration used in UKESM1 with a bulk soil biogeochemistry . They found that JULES-ES has a relatively small increase in NPP caused by the addition extra N in the form of deposition compared with both the meta-analyses and CLM/LPJ-GUESS. However, it is comparable to that found in JSBACH. This small response is, in part, caused by the smaller initial N limitation in JULES-ES. However, JULES' increase in NPP in response to CO${}_{\mathrm{2}}$ fertilisation is aligned with the majority of the models and the meta-analyses.

7 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 the model, N limitation affects NPP and how much C 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 N 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) is only one of the new components of JULES that has been included within UKESM1 . Relevant additions to the JULES-ES configuration used in UKESM1 includes more plant functional types with improved plant physiology and vegetation dynamics plus a new land use module .

Overall, the N-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 JULES-C and the N-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 N availability. What JULES-C lacks is a mechanism for this to change substantially in time – either under more limiting conditions as elevated CO${}_{\mathrm{2}}$ outpaces demand for nutrients e.g. or under conditions of increased N availability due to anthropogenic deposition or accelerated soil decomposition caused by climate change leading to increased mineralisation rates . The response of the N cycle in JULES under changes in climate and CO${}_{\mathrm{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 additionally includes the vertical discretisation of the soil biogeochemistry model. This configuration improves the ecosystem residence times in the tundra and boreal regions. 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.

Code availability

The JULES code used in these simulations is available from the Met Office Science Repository Service (registration required) at https://code.metoffice.gov.uk/trac/jules (last access: 9 April 2021). To access the code a freely available non-commercial research licence is required (https://jules-lsm.github.io/, last access: 9 April 2021). The suites required for running JULES are available here: https://code.metoffice.gov.uk/trac/roses-u (last access: 9 April 2021). JULES-CN is u-ah896, JULES-C is u-ah932, and JULES-CN${}_{\text{layer}}$ is u-ai571.

Data availability

JULES output data are available at 10.5281/zenodo.4692460 .

Author contributions

AW designed the model in collaboration with the rest of the co-authors and wrote the first draft of the text. EB and SC added the layered soil biogeochemistry. PC performed the analysis of inter-annual variability. All authors reviewed the paper and proposed improvements.

Competing interests

The authors declare that they have no conflict of interest.

Financial support

This research has been supported by the European Union Horizon 2020 (grant no. 641816 and 647204), the Joint UK BEIS/Defra Met Office Hadley Centre Climate Programme (grant no. GA01101), the Natural Environment Research Council (grant nos. NE/M01990X/1 and NE/R015791/1), and the European Research Council (grant no. 742472).

Review statement

This paper was edited by Jatin Kala and reviewed by William Wieder, David Wårlind, and one anonymous referee.