Articles | Volume 14, issue 4
Model description paper
27 Apr 2021
Model description paper |  | 27 Apr 2021

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

Andrew J. Wiltshire, Eleanor J. Burke, Sarah E. Chadburn, Chris D. Jones, Peter M. Cox, Taraka Davies-Barnard, Pierre Friedlingstein, Anna B. Harper, Spencer Liddicoat, Stephen Sitch, and Sönke Zaehle

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-CNlayer). 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-CNlayer, 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 CO2 fertilisation effect, such that the simulated current-day land C sink is reduced by about 0.5 Pg C yr−1 compared to the version with no N limitation. JULES-CNlayer 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 (Le Quéré et al.2018; Friedlingstein et al.2020), 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 (CO2) 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 (Weintraub and Schimel2005).

Any increase in atmospheric CO2 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, Zaehle (2013) 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 CO2 and temperature effects on the land C sink (Wenzel et al.2016; Cox et al.2013). In addition, climate projections assessed by the IPCC using CMIP5 Earth system models that lacked terrestrial carbon cycle (Ciais et al.2014) have been shown to exhibit a major and systematic bias in their future projection of land carbon sink (Zaehle et al.2015; Wieder et al.2015b). An increasing number of land surface and climate models now include constraints on the land C sink caused by N limitation (Zaehle et al.2014; Wania et al.2012; Smith et al.2014). In fact, recent simulations have generated a range of estimates for the sensitivity of the C cycle to N availability (Meyerholt et al.2020a; Davies-Barnard et al.2020; Arora et al.2020). For example, Meyerholt et al. (2020a) used a perturbed model ensemble to show that N limitation reduces both the projected future increase in land C store due to CO2 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 (Arora et al.2020). Jones and Friedlingstein (2020) 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 (Davies-Barnard et al.2020).

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 (Best et al.2011; Clark et al.2011) (JULES at vn5.1 –, last access: 9 April 2021). JULES is the land surface component of the UK Earth System Model (UKESM) (Sellar et al.2019). 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 (Harper et al.2018), an enhanced representation of surface exchange and hydrology (Wiltshire et al.2020) and a new module for land management (Burton et al.2019). These separate components combine to make the land surface component of UKESM and were used for the most recent Global Carbon Budget annual assessment (Friedlingstein et al.2020).

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 (Cox et al.1998, 1999; Cox2001; Clark et al.2011). 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. 3.1.3 and by Wiltshire et al.2020). 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 (Cox2001), 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 Ψ) becomes non-structural carbohydrates, root exudates, or biogenic volatile organic compounds (Collalti and Prentice2019). 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) (Meyerholt and Zaehle2015) 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 (Brix and Ebell1969; Wang et al.2012; Field and Mooney1986; McGuire et al.1995). However, more recent analyses do make the link between nutrient availability and leaf level N concentrations (e.g. Mao et al.2020). 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 (Zaehle and Dalmonech2011). 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 (Arora et al.2020). 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 (Sellar et al.2019). JULES can also be run offline forced by observed meteorology globally, regionally, or at a single location. A full description of the main components of JULES is provided by Best et al. (2011) and Clark et al. (2011). 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) (Cox2001). In this version of TRIFFID, five plant functional types (PFTs) are represented: broadleaf trees, needleleaf trees, C3 grasses, C4 grasses, and shrubs. The soil C model is based on the Roth-C model (Clark et al.2011). Recently, Burke et al. (2017) and Chadburn et al. (2015) 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 (Le Quéré et al.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-CNlayer. The soil biogeochemistry is represented by a single bulk layer in JULES-CN, whereas it varies as a function of depth in JULES-CNlayer. 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 1Default values of PFT-specific parameters for allometry, allocation, and vegetation N and C stoichiometry in the JULES-CN and JULES-CNlayer configurations. The subscript (i) is present to show that it is a PFT-specific value; nl0,i is the N concentration at the top of the canopy but is shown here as 1/nl0,i so that it is comparable to expected C : N ratios from the literature.

Download Print Version | Download XLSX

3 JULES developments

3.1 Vegetation C and N

The TRIFFID dynamic vegetation model provides the core of the vegetation module (Cox2001). 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 (Clark et al.2011). 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. 3.1.1) and including the additional complexity of labile N (Sect. 3.1.2). The following subsection describes how growth and spreading is limited by N availability (Sect. 3.1.3). 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. 3.1.4). Biological N fixation is input directly into the soil inorganic N pool and is described later in Sect. 3.3.1.

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 (b,i in m2 m−2). b,i is the prognostic variable used in JULES to describe the vegetation and is functionally the equivalent of the potential leaf area. Given b,i, leaf, root, and wood pools are diagnosed for each PFT as introduced in Cox (2001). 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 Clark et al. (2011).

Figure 1Stoichiometry 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.


The vegetation C density per PFT (Cv,i in kg C m−2) can be separated into leaf (Lc,i in kg C m−2), fine root (Rc,i in kg C m−2), and total stem plus coarse root (Wc,i in kg C m−2) pools, each of which is related allometrically to the balanced leaf area (b,i). Each component is then related to b,i. Root C is set equal to leaf C, which is itself a linear function of b,i, and total stem C is related to b,i by a power law (Enquist et al.1998):


where σl,i (kg C m−2), awl,i (kg C m−2), and bwl,i (dimensionless) are PFT-dependent allometric parameters defined in Table 1. By definition 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 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. (4) with a “pipe model” approach (Shinozaki et al.1964a, b) to obtain the canopy height for PFT i (hi in m):

(5) h i = W c , i a ws , i η sl , i a wl , i W c , i 1 / b wl , i ,

where ηsl,i (kg C m−2 per unit LAI) relates respiring stem to leaf C (Table 1) and aws,i is the ratio of total stem C to respiring stem C. We can combine Eqs. (4) and (5) to relate (b,i) to canopy height (hi) 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 (Clark et al.2011). We extend their use to explicitly define N pools as part of the new scheme:


where μrl,i and μsl,i are dimensionless stoichiometric parameters linking the top leaf N concentration (nl0,i in kg N kg C −1) to the total stem and root N concentration via nl0,i. The leaf N pool (Ln,i in kg N m−2) has an additional dependency on phenological state (Sect. 3.1.2) and assumed distribution of N in the canopy. Following Eq. (1), the total vegetation N store per PFT (Nv,i in kg N m−2) is given by

(8) N v , i = L n , i + R n , i + W 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. 1). 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. 2).

Figure 2Total 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.


3.1.2 Labile C and N: phenology and mobilisation

The total leaf C pool per PFT (Lc,i, Eq. 2) 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. Lc,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 (pi), which relates the LAI (i) at any moment in time to the balanced leaf area index (b,i):

(9) L i = p i L b , i ,

where pi is a scalar between 0 and 1 that describes the phenological state of the system (Clark et al.2011). For evergreen plants pi is a constant of 1. The two state variables b,i and pi 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 (Ln,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 n , i = p i n lc , i L c , i + ( 1 - p i ) 1 + κ l , i 2 n lc , i L c , i ,

where κl,i is the dimensionless leaf N retranslocation coefficient and nlc,i is the mean canopy N concentration (defined in Eq. 11 in kg N kg C −1). Here κl,i is set to 0.5 for all PFTs (Zaehle and Friend2010). 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-CNlayer 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 (Clark et al.2011). In the application here, we use this distribution to be fully consistent with the physiology. The vertical distribution of leaf N concentration (nlc,i(d) in kg N kg C −1) in the canopy is described by Mercado et al. (2007)

(11) n lc , i ( d ) = n l 0 , i e - k n , i d ,

where kn,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 (Carswell et al.2000), a value of 0.78 for kn,i was found (Mercado et al.2007). Equation (11) 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 (Cv,i, Eq. 1) and vegetation N (Nv,i, Eq. 8) 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 (Ra). In JULES-CN the potential NPP or NPPpot 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 (Ψ in kg C m−2). A positive Ψ results in a reduction of carbon use efficiency (CUE).

The C balance per PFT i is given by

(12) d C v , i d t = ( 1 - λ i ) Π i - Λ c , i - Ψ g , i ,

where Πi is the potential NPP per unit area of PFT in (kg C m−2 s−1 (prior to nutrient limitation) and Λc,i (kg C m−2 s−1) is the PFT-specific litterfall rate (Sect. 3.1.4). Any excess C from growth (Ψ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. λi is the coefficient for partitioning the NPP between growth and spreading. λi is utilised in increasing the fractional coverage of the vegetation, and (1−λi) increases the C of the existing vegetated area. λi is a function of the vegetation C which itself is a function of the balanced LAI for PFT i (b,i):

(13) λ i = 1 L b , i > L max , i L b , i - L min , i L max , i - L min , i L min , i < L b , i L max , i 0 L b , i L min , i .

The equivalent N balance per PFT is given by

(14) d N v , i d t = ( 1 - λ i ) Φ i - Λ n , i ,

where Φi (kg N m−2 s−1) is the PFT-specific N uptake (see Eq. 19) and (1−λii is equal to Φg,i, the N uptake available for growth. Λ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 λi of the total available N with a fraction (1−λi) available for growth.

Litter is produced by the turnover of the leaf, wood and root pools for each PFT, defined as

(15) Λ c , i = γ l , i L c , i + γ r , i R c , i + γ w , i W c , i ,


(16) Λ n , i = ( 1 - κ l , i ) γ l , i L n , i + ( 1 - κ r , i ) γ r , i R n , i + γ w , i W n , i ,

for litter C (Λc,i in kg C m−2 s−1) and litter N (Λn,i in kg N m−2 s−1), respectively. γr,i and γw,i are turnover rates in s−1 (Table 6 of Clark et al.2011). The leaf turnover rate (γl,i) is a temperature-dependent turnover rate consistent with the phenological state and defined in Clark et al. (2011). 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. κl,i and κr,i are the dimensionless coefficients for the retranslocation of leaf and root N shown in Table 1 (Zaehle and Friend2010).

In JULES-CN the N available for plant uptake for each PFT i (Navail,i in kg N m−2) is the inorganic soil N pool (Nin in kg N m−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-CNlayer is more complicated and takes into account the soil profile. This is discussed in Sect. 3.3.2.

Equations (12) and (14) have two remaining unknowns for each PFT: the plant N uptake for growth (Φg,i) and the excess C from growth (Ψ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 Ψg,i=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 (Δt) (Φg,i<(1-λi)Navail,i/Δ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 (Φg,i=(1-λi)Navail,i/Δt) and the excess C for growth Ψg,i can be derived. Following the solution of dNv,idt the C store and balanced LAI (b,i) are updated and the leaf, root, and wood pools for each PFT can be derived following the allometric equations (Eqs. 24).

The remaining proportion (λ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) Φ s , i = N v , i C v , i Π i - d C v , i d t - Ψ s , i ,

where Ψs,i (or λiΨi) is the excess C term from spreading and Nv,iCv,i is the inverse of the the whole plant C : N ratio. As with growth limitation, Eqs. (17) is first solved to find the N demand for spreading (Ψs,i=0.0). If the arising demand is less than the available N (Φs,i<λiNavail,i/Δt), spreading is unlimited. If N demand is in excess of that available, the uptake is set equal to the available N flux (Φs,i=λiNavail,i/Δt) and the excess (Ψs,i) assimilate is solved for.

Total excess C per PFT i (Ψi) is therefore the combination of that from growth plus spreading:

(18) Ψ i = Ψ s , i + Ψ g , i .

Similarly, total N uptake per PFT i (Φi) is therefore the combination of N uptake used for growth plus N uptake used for spreading:

(19) Φ i = Φ s , i + Φ g , i .

The PFT level N uptake and excess C are weighted by the fraction of coverage of each PFT in a grid cell (vi) and summed to get the grid-averaged values:


This excess C (Ψ) 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 Clark et al.2011, 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 (vi).

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 (Λc, kg C m−2 s−1) is made-up of the area-weighted sum of the litter C from each PFT (Λ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) Λ c = i v i Λ c , i + γ v , i C v , i + ( Π i - Ψ i ) j c i , j v j ,

where ci,j are the competition coefficients describing the effect of PFT i on PFT j, γv,i is a large-scale disturbance term of PFT i, vi is the vegetation fraction of PFT i, Πi is defined in Eq. (12) and Ψi in Eq. (18). The effect of N limitation on the litter C flux is captured in the excess C term per PFT (Ψi). Similarly to the total litter C, total litter N (Λn, kg N m−2 s−1) is given by

(23) Λ n = i v i Λ n , i + γ v , i N v , i + Φ i j c i , j v j .

Both Λc and Λ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 (Jenkinson et al.1990; Jenkinson and Coleman1999).

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 fDPM,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 (βR) enters the BIO and HUM pools. The C pools are updated according to


where t is the time in s, Cp are the C pools in kg C m−2 (where p is one of DPM, RPM, BIO, HUM), Λc,i is the litter input for PFT i in kg C m−2 s−1 (term in brackets in Eq. 22), fDPM,i represents the fraction of litter from each PFT i that goes into DPM with the rest (1-fDPM,i) going into the RPM pool (dependent on amount of woody vegetation), and Rtot is the total turnover in kg C m−2 s−1, where the Rp represent the turnover of each C pool:

(28) R tot = R DPM + R RPM + R BIO + R HUM .

The soil respiration to the atmosphere (rh in kg C m−2 s−1) is given by

(29) r h = ( 1 - β R ) R tot ,

where β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) β R = 1 4.09 + 2.67 e ( - 0.079 clay ) .

For each C pool there is an equivalent N pool, with the N pools following a similar structure to the C pools:


Inputs into the litter pools (DPM, RPM) are from the litter N flux (Λn,i in kg N m−2 s−1, Eq. 23), and losses are determined by the pool-specific mineralisation of organic N into inorganic N (Mp in kg N m−2 s−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 (Itot in kg N m−2 s−1):

(35) I tot = I DPM + I RPM + I BIO + I 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 (Rp,pot)

(36) R p , pot = k p C p F T ( T soil ) F θ ( θ ) F v ( v ) ,

where the kp are fixed constants in s−1 (Clark et al.2011). The functions of temperature (FT(Tsoil)) and moisture (Fθ(θ)) depend on the temperature (Tsoil) and moisture content (θ) near the soil surface. The function Fv(v) depends on the vegetation cover fraction (v) (Clark et al.2011). The potential mineralisation of organic N when the system is not N-limited (Mp,pot) is related to the potential turnover rates by the C to N ratio of each pool (CNp):

(37) M p , pot = R p , pot CN p .

Similarly, the potential immobilisation of inorganic N into the organic N pools (Ip,pot) is related to pool potential turnover (Rp,pot), the retained fraction of respiration (βR), and the C to N ratio of the destination pool in the decomposition chain:

(38) I p , pot = β R R p , pot CN soil ,

where CNsoil 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 (Mp) and potential immobilisation (Ip) 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 ( T soil ) F θ ( θ ) F v ( v ) F N ,

where p is one of RPM or DPM. The nitrogen-limited mineralisation and immobilisation of the DPM and RPM pools (Eqs. 41 and 37) are now effectively a function of Rp.

FN 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. 40). FN is limited to a range of 0.0 to 1.0. When FN is equal to 1, the decomposition, mineralisation, and immobilisation take place at the potential rate and the system is not N-limited. Where FN 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. FN is given by

(40) F N = ( M BIO + M HUM - I BIO - I HUM ) Δ t + N in ( D DPM + D RPM ) Δ t ,

where Nin is the total soil inorganic N pool (in kg N m−2; discussed in Sect. 3.3 and defined in Eq. 51) and Δt is the time step. DDPM and DRPM are the net demand associated with decomposition of each of the litter pools:

(41) D p = I p , pot - M p , 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 FN<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 FN is effectively 1.0).

If the net mineralisation is positive some of the N is emitted as gas, according to

(42) N gas = f gas ( M tot - I tot ) ,

where Ngas is one component of the gas emission in kg N m−2 s−1, fgas is a parameter that sets the fraction of the N flux that is emitted as gas to the atmosphere. Following Thomas et al. (2013a), it is assumed that 1 % of net mineralisation is emitted as gas (fgas is set to 0.01). Mtot is the total mineralisation flux in kg N m−2 s−1:

(43) M tot = M DPM + M RPM + M BIO + M HUM .

If pool sizes become too small, Ngas could become negative to ensure N is conserved.

3.2.1 Vertical discretisation

The vertical discretisation of the soil C and N follows Burke et al. (2017). There is a set of four soil C and N pools (DPM, RPM, BIO, HUM) in every soil model layer. As in Burke et al. (2017) 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 , pot ( z ) = k p C p ( z ) F T ( T soil ( z ) ) F θ ( θ ( z ) ) F v ( v ) exp ( - ξ resp z ) .

FT(Tsoil(z)), Fθ(θ(z)) and Cp(z) are now all dependent on depth. Tsoil(z) and θ(z) are the simulated layered soil temperature and soil moisture content, and Cp(z) is the simulated soil C content for each layer and pool p. The additional reduction of turnover with depth is exponential, with ξresp an empirical parameter (in m−1) that controls the magnitude of the reduction (Burke et al.2017). The larger the value of ξresp, the more inhibited the respiration is with increasing depth. Here ξresp was tuned to give a realistic estimate of soil C in a vertically resolved version of JULES-C as in Burke et al. (2017). When N is limiting, the respiration of the DPM and RPM pools are reduced by a factor of FN(z), which is also now a function of depth and dependent on the available N in the relevant layer. Mp and Ip 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) N DPM ( z ) t = z D ( z ) N DPM ( z ) z + i v i f DPM , i Λ n , i f lit ( z ) - M DPM ( z )

(46) N RPM ( z ) t = z D ( z ) N RPM ( z ) z + i v i ( 1 - f DPM , i ) Λ n , i f lit ( z ) - M RPM ( z )

(47) N BIO ( z ) t = z D ( z ) N BIO ( z ) z + 0.46 I tot ( z ) - M BIO ( z )

(48) N HUM ( z ) t = z D ( z ) N HUM ( z ) z + 0.54 I tot ( z ) - M HUM ( z )

Itot(z) is the total immobilisation in kg N m−2 s−1 in each layer (following Eq. 35). D(z) is the diffusivity in m2 s−1 and varies both spatially and with depth (z) (Burke et al.2017).

(49) D ( z ) = D o ; z 1 m D o 2 ( 3 - z ) ; 1 m < z < 3 m 0.0 ; z 3 m

Do (m2 s−1) varies spatially depending on the freeze and thaw state of the soil. In regions without permafrost, Do represents a bioturbation mixing rate equivalent to 1 cm2 yr−1. When permafrost is present, Do represents cryoturbation and increases to a value equivalent to 5 cm2 yr−1. The parameterisation of D(z) in Eq. (49) 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 (flit(zn,i for litter N) are distributed so that they decline exponentially with increasing depth. Here flit(z) is independent of the PFT type and hence the root distribution is as follows:

(50) f lit ( z ) = e - ξ lit z z = 0 z max e - ξ lit z d z ,

where ξlit is the parameter to reduce the litter input with increasing depth and is set to 0.2 m or 5 m−1 and z is the mid-point of each layer.

The mineralised gas emissions are now a function of depth (Ngas(z)) and are calculated by repeating Eq. (42) for each soil layer. Similarly, the litter decomposition rate modifier (FN) is calculated by repeating a slightly modified version of Eq. (40) for each soil layer. In the vertically resolved version of Eq. (40), if the soil layer is frozen Nin 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) d N in d t = N dep + i v i BNF i - i v i Φ i + M net - N leach - N gasI ,

where Nin is the inorganic N in kg N m−2, Ndep is prescribed N deposition in kg N m−2 s−1 and vi the fractional cover of each PFT i. The biological N fixation (BNFi) for each PFT i is described in Sect. 3.3.1 below, and plant uptake (Φi) for each PFT i is described in Sect. 3.1.3. Mnet is the net mineralisation, which is the difference between Mtot (Eq. 43) and Itot (Eq. 35) reduced by Ngas (Eq. 42). The loss of N from the system via the inorganic pool is the sum of leaching (Nleach in kg N m−2 s−1) plus an additional gas loss to the atmosphere (NgasI in kg N m−2 s−1):

(52) N gasI = γ n N in ,

where γn is a tunable parameter (in s−1). The total N gas loss is the sum of NgasI above and Ngas from Eq. (42), with NgasI representing approximately 90 % of the total gas loss. This additional gas loss term (NgasI) represents missing processes relating to the gaseous loss of inorganic N and limits the effective mineral N pool size. Including NgasI 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 NPPpot. In the current model configuration γn is set to 0.0028 d−1 such that the whole pool turns over once every model year.

The leaching of N (Nleach in kg N m−2 s−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 (α, dimensionless) representing the effective solubility of N. α is assumed to have a value of 0.1 and in JULES-CN represents the combined sorption of all inorganic N species (Wania et al.2012).

(53) N leach = α ( N in / θ 1 m ) Q subs ,

where θ1 m is the soil water content in the top 1 m of soil in kg m−2 (so the inorganic N is assumed to occupy the top 1 m of soil) and Qsubs is the total subsurface runoff in kg m−2 s−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 Cleveland et al. (1999), N fixation is determined as a linear proportion of the net primary production before N limitation of each PFT i (NPPpot,i):

(54) BNF i = ζ NPP pot , i .

NPPpot,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 (Ψ) is removed. BNF as a function of NPP is an established method used and assessed in other models (Meyerholt et al.2016; Wieder et al.2015a; Thomas et al.2013b). While some models utilise more complex BNF representations (Fisher et al.2010), 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 (Davies-Barnard and Friedlingstein2020). However, changes in NPP may not accurately reflect changes in BNF with forcings such as elevated atmospheric carbon dioxide (Liang et al.2016) or additional N (Thomas et al.2013b; Ochoa-Hueso et al.2013).

The rate of fixation (ζ) is set such that global present-day net primary productivity of approximately 60 Pg C yr−1 results in approximately 100 Tg N yr−1 fixation (0.0016 kg N kg C−1), within the range of recent global observation-based estimates of BNF (Davies-Barnard and Friedlingstein2020; Vitousek et al.2013). 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 (Houlton et al.2008; Cleveland et al.1999) though not recent observational meta-analyses (Davies-Barnard and Friedlingstein2020).

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-CNlayer 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-CNlayer there is an inorganic N pool in each soil layer. The dynamics are very similar to Eq. (51), but most of the components now vary with depth:

(55) d N in ( z ) d t = i v i BNF i f R , i ( z ) - i v i Φ i f I , i ( z ) + M net ( z ) - N flux ( z ) - N gasI ( z ) .

Any N deposition (Ndep) 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. (55) are fR,i(z) – the fraction of roots in each layer for PFT i (Eq. 56), fI,i(z) – the fraction of available inorganic N in each layer for PFT i (Eq. 60), and Nflux(z) – the transport of inorganic N through the layer by the soil water fluxes (Eq. 61).

As in Eq. (51), the net mineralisation flux (Mnet(z)) is the difference between Mtot(z) and Itot(z) reduced by Ngas(z) for each layer (see Sect. 3.2.1). Inputs from biological N fixation from PFT i are distributed according to the root profile of the PFT under consideration (fR,i(z)):

(56) f R , i ( z ) = f root , i ( z ) z = 0 z max f root , i ( z ) d z ,

where froot,i(z) is the volumetric root fraction of PFT i at a given soil level and zmax is the maximum depth of the soil in metres. Gas loss from the inorganic N (NgasI(z)) occurs in each layer but with an additional exponential decay term that is a function of depth (similar to that used in Eq. 44 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 gasI ( z ) = γ n N in ( z ) e - ξ resp z .

This leaves two terms in Eq. (55): the plant uptake term (iviΦifI,i(z)), which is PFT-dependent, and the Nflux(z) term, which replaces the leaching term from Eq. (51). 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, froot,i(z) (which reduces with depth). So for each PFT, i, the available amount of the inorganic N pool (Navail,i(z) in kg N m−2) at equilibrium that could potentially be extracted by the vegetation is given by

(58) N avail , i ( z ) = f root , i ( z ) N in ( z ) T ( z ) ,

where T(z) is zero when the soil temperature is 0 C or colder and 1 when it is above 0 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) N avail , i ( z ) d t = γ dif ( f root , i ( z ) N in ( z ) - N avail , i ( z ) ) ,

where γdif is the rate of diffusion back to the equilibrium, set by default to 0.28 d−1 or approximately 100 yr−1. Navail,i(z) is then multiplied by T(z) 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 ( z ) = N avail , i ( z ) z = 0 z max N avail , i ( z ) 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 δz, the inorganic N flux (Nflux) is given by

(61) N flux ( z ) = α δ z d d z W flux ( z ) N in ( z ) θ ( z ) ,

where θ(z) is the soil water content of the layer in kg m−2 and Wflux(z) is the flow rate of the water through the layer in kg m−2 s−1. Multiplying by δ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 2 summarises the extra parameters required for the soil biogeochemistry component of JULES-CN and JULES-CNlayer alongside their values.

Table 2A summary of the extra parameters required for the soil biogeochemistry component of JULES-CN and JULES-CNlayer.

Download Print Version | Download XLSX

4 Historical simulations

Global transient simulations were carried out following the protocol for the S2 experiments in TRENDY (Sitch et al.2015), which include time varying climate, CO2, and N deposition but pre-industrial land use. Time-varying CO2 and climate were from the from the CRU-NCEP dataset (v7,, Viovy2008). The fraction of agriculture in each grid cell was set to the pre-industrial value defined by (Hurtt et al.2011). N deposition was taken from a ACCMIP multi-model dataset interpolated to annual fields (Lamarque et al.2013). The model resolution was N96 (1.25 latitude × 1.875 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 Le Quéré et al. (2018).

  • JULES-CNlayer 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, C3 and C4 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 Burke et al. (2017) and Chadburn et al. (2015) 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-CNlayer 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−1 globally. The soil C distribution in JULES-CNlayer is particularly slow to reach equilibrium. Therefore the “modified accelerated decomposition” technique (modified AD) described by Koven et al. (2013) was used to spin the soil C in these versions up to an initial pre-industrial equilibrium distribution (Burke et al.2017). 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−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-CNlayer 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-CNlayer are assessed. For completeness, figures often include both JULES-CN and JULES-CNlayer, but JULES-CNlayer is only discussed at the end of the results.

Figure 3Total 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 Poulter et al. (2015) converted to JULES PFTs by Hartley et al. (2017).


Figure 4Carbon and nitrogen stocks and fluxes for JULES-CN, JULES-C, and JULES-CNlayer for the period 1996–2005 (after Davies-Barnard et al. (2020)). C is carbon, N is nitrogen; rh is heterotrophic respiration, ra 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 Ngas and NgasI, with NgasI 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: Hashimoto et al. (2015); (b) TER: Li et al. (2018); (c) TER: Ballantyne et al. (2017); (d) GPP: Jung et al. (2011); (e) vegetation carbon and SOM + litter carbon: Carvalhais et al. (2014); (f) BNF Davies-Barnard and Friedlingstein (2020); (g) BNF Vitousek et al. (2013); (h) Vegetation nitrogen: Schlesinger (1997); (i) NPP: Zhao and Running (2010); (j) soil organic nitrogen: Batjes (2014); (k) soil organic nitrogen: Global Soil Data Task Group (2000); (l) nitrogen losses including nitrogen leaching: Gruber and Galloway (2008); (m) nitrogen leaching: Boyer et al. (2006); (n) nitrogen leaching: Galloway et al. (2004); (o*) organic nitrogen immobilisation and mineralisation and plant uptake von Bloh et al. (2018); (p*) nitrogen uptake, vegetation nitrogen, and nitrogen emissions Zaehle et al. (2010); (q*) nitrogen uptake and inorganic nitrogen content Xu-Ri and Prentice (2008); and (r*) nitrogen uptake and total nitrogen emissions Wania et al. (2012). (o*), (p*), (q*), and (r*) are model-derived estimates.


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 Sellar et al. (2019) the N cycle was combined with a new competition scheme Harper et al. (2018) 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 3 shows the total area covered by each type of vegetation. The Climate Change Initiative (CCI) land cover observations Hartley et al. (2017) are added for completeness. In general, the models all tend to overestimate the shrubs and underestimate the grass. However, Sellar et al. (2019) 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. 10), this is because the tropical forests dominate the tree region and their growth is not limited by N in the model. JULES-CNlayer 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 4 provides an overview of the stocks and fluxes of C and N in JULES-CN and JULES-CNlayer 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  12 % higher in JULES-C than in JULES-CN. There is also a small reduction in the GPP of  4 % caused by some differences in the vegetation fractional cover distribution (Fig. 3) and indirectly resulting from the N limitation.

Figure 5Zonal total values of (a) net primary productivity (NPP) and (b) gross primary productivity (GPP) for JULES-C, JULES-CN, and JULES-CNlayer simulations for the period 1996–2005 in Pg C (degree latitude)−1 yr−1. The observational constraint for NPP is from MODIS (Zhao and Running2010) and that for GPP is from Jung et al. (2011). The zonal mean carbon use efficiency (CUE = NPP / GPP) is shown in (c). The CUE observational constraint was digitised from Kim et al. (2018). 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.


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 yr−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 yr−1 for JULES-CN. Leaching is fairly low at 7 Tg N yr−1 compared to estimates of leaching, which are as high as  25 %–55 % of N inputs in European forests (Dise et al.2009) and range between 59 and 118 Tg N yr−1 in the available observations (Boyer et al.2006; Galloway et al.2004). There is no N fertiliser applied in the model, which might partially explain why the leaching is so low. In reality there is  200 Tg N applied annually as either manure or fertiliser Potter et al. (2010), 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 Xu-Ri and Prentice (2008), Smith et al. (2014), Zaehle (2013), and von Bloh et al. (2018).

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 5 shows the zonal total GPP and NPP for JULES-CN and JULES-C. As expected from Fig. 4, 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 5 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 Kim et al. (2018) observation-based dataset with a bias of  0.09. This bias is relatively constant with latitude. However, considerable uncertainties remain in these estimates.

Figure 5 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 CO2 fertilisation, but this is limited by N in the JULES-CN configuration.

The zonal total soil and vegetation C stocks are shown in Fig. 6. The vegetation C is very similar for both JULES-C and JULES-CN, as expected from Fig. 4, 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 6Zonal total values of (a) vegetation and (b) soil C for JULES-C, JULES-CN, and JULES-CNlayer simulations for the period 1996–2005 in Pg C (degree latitude)−1. For the vegetation C the observation-based constraints are Saatchi: Saatchi et al. (2011); GEOCARB: Avitabile et al. (2016); and biomass: Ruesch and Gibb (2008). The observation-based constraints for the soil carbon are IGBP-DIS: Global Soil Data Task Group (2000); WISE: Batjes (2016); and Carvahlais: Carvalhais et al. (2014).


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 (Cox et al.2013). As a result, simulation of NEE variability is highly relevant to climate–carbon cycle projections (Wenzel et al.2016).

Figure 7Evaluation of global annual mean NEE from JULES-CN, JULES-C, and JULES-CNlayer compared with observations based on estimates from Global Carbon Project (GCP) (Friedlingstein et al.2020) 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 7 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 (Friedlingstein et al.2020). 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 (Fa) plus the ocean (Fo):

(62) NEE gcp = FF + LU - F a - F o .

Figure 8Biome-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-CNlayer are shown for the period 1996–2005. The land surface is split into biomes based on the 14 World Wildlife Fund terrestrial ecoregions (Olson et al.2001) and characterised by Harper et al. (2018). The observed ecosystem residence times are derived from the Carvalhais et al. (2014) global dataset and the observed soil residence times are from the WISE (Batjes2016) soil carbon combined with the Hashimoto et al. (2015) soil respiration.


The observations and both of the models show an upward trend in NEE but with very significant interannual variability (Fig. 7). Due to N limitations on CO2 fertilisation, mean NEE in JULES-CN (1.66 Pg C yr−1) is lower than in JULES-C (2.06 Pg C yr−1) and also lower than the estimate from GCP (2.11 Pg C yr−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 Sellar et al. (2019). 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 3). The correlation coefficient for NEE between the JULES-CN and GCP estimates (r=0.71) is also improved compared to JULES-C (r=0.63). There remains a significant underestimate of NEE in the years following the Pinatubo volcanic eruption in 1991, most likely due to the neglect of diffuse-radiation fertilisation in these versions of JULES (Mercado et al.2009). However, it is especially notable that JULES-CN significantly reduces the systematic overestimate of NEE seen in JULES-C during extended La Niña periods, such as the years centred around 1974 and 2000 (Fig. 7).

Table 3Statistics of NEE from JULES-CN, JULES-C, JULES-CNlayer, and the GCP observation-based estimates (Friedlingstein et al.2020), 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.

Download Print Version | Download XLSX

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 8 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 (Olson et al.2001) and characterised by Harper et al. (2018). The ecosystem residence time defined as the total ecosystem C divided by the GPP is shown in Fig. 8a. 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 Carvalhais et al. (2014). 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-CNlayer, the residence times in the boreal and tundra increase notably (see Sect. 3.2.1 for further discussion).

Figure 9The spatial distribution of the response ratio, defined as the potential amount of carbon that can be allocated to growth and spreading of the vegetation (NPPpot), as a fraction of the NPP achieved in the natural state for (a) JULES-CN and (b) JULES-CNlayer. 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−2 are set to missing. This is the spatial distribution of the metric shown in Fig. 10. 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 10The response ratio is the ratio of the potential amount of carbon that can be allocated to growth and spreading of the vegetation (NPPpot) compared with the actual amount achieved in the natural state (NPP). As in Fig. 9, any grid cells with an annual NPP of less than 0.016 g C m−2 are set to missing. The median of JULES-CN and JULES-CNlayer are shown for each biome for the period 1996–2005. The biomes are discussed in more detail in Fig. 8. The observational constraint is taken from Table 1 in LeBauer and Treseder (2008), 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.


5.3 Impact of N limitation

IN JULES-CN and JULES-CNlayer 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 (NPPpot) compared with the actual amount achieved in the natural state (NPP). Both of these diagnostics are output from the JULES simulations. Figure 9 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 9c 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. 5f.

Figure 10 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, LeBauer and Treseder (2008) 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 (Yang et al.2014); 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. 11) shows a similar distribution to the soil organic C (Fig. 6), reflecting the relatively consistent C to N ratio of the soil within the model. CNsoil – 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. 11). Figure 12 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. 9. 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 13 shows a slight increase in the N demand and N uptake over the 20th century associated with the increase in vegetation growth (Fig. 5). Similarly, there is an increase in the BNF that is parameterised such that it is proportional to the NPP.

Figure 11The zonal total soil organic and inorganic nitrogen stocks in Pg N (degree of latitude)−1. JULES-CNlayer shows the stocks for the top 1 m of soil. The observations of nitrogen stocks are from Global Soil Data Task Group (2000).


Figure 12Spatial distribution of N fluxes for JULES-CN for the period 1996–2005. JULES-CNlayer is not shown because the spatial patterns are very similar to those for JULES-CN.

5.5 Impact of vertical discretisation of soil biochemistry

This section discusses the differences between JULES-CN and JULES-CNlayer. In general, over the tropics and southern latitudes, JULES-CNlayer 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. 3 is caused by a reduction in the boreal regions, which have a larger proportion of shrubs and grasses in JULES-CNlayer. In the higher latitudes the soil in JULES-CNlayer also has more organic C (Fig. 6). This increase in soil organic C represents a store of permafrost carbon more comparable to the carbon found by Batjes (2014) and Carvalhais et al. (2014). This build-up of carbon in JULES-CNlayer 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. 8b. The modelled soil C residence time in JULES-CNlayer is now much longer and more comparable to that observed.

The spatial distributions of N fluxes in JULES-CNlayer (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. 13). The main differences are in the N gas loss, which is larger in JULES-CNlayer, and the N leaching, which is larger in JULES-CN. Figure 11 shows an increase in both organic and inorganic N in JULES-CNlayer 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-CNlayer 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-CNlayer (Fig. 11) 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 13N fluxes for JULES-CN and JULES-CNlayer over the historical period.


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 (Sterner and Elser2002). This assumption has some support in the literature (e.g. Brix and Ebell1969; Wang et al.2012) and is a common approach amongst complex DGVMs (Meyerholt and Zaehle2015). However, recent meta analyses of field observations show a distinct increase in foliar N to additional N availability (Mao et al.2020) 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 (Meyerholt and Zaehle2015). 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 (Harper et al.2016) 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 Davies-Barnard et al. (2020) and Vitousek et al. (2013), the spatial distribution of BNF more heavily favours the tropics than recent observations suggest (Sullivan et al.2014; Davies-Barnard et al.2020). 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.) (Zheng et al.2020). 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 Hugelius et al. (2020). 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 CO2 fertilisation from climate change or from the impact of N deposition. However, this was explored by Davies-Barnard et al. (2020), 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. Davies-Barnard et al. (2020) used a slightly different configuration of JULES (JULES-ES), which is the configuration used in UKESM1 with a bulk soil biogeochemistry (Sellar et al.2019). 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 CO2 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 (Sellar et al.2019). Relevant additions to the JULES-ES configuration used in UKESM1 includes more plant functional types with improved plant physiology and vegetation dynamics (Harper et al.2016) plus a new land use module (Burton et al.2019).

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 CO2 outpaces demand for nutrients (e.g. Zaehle2013) or under conditions of increased N availability due to anthropogenic deposition or accelerated soil decomposition caused by climate change leading to increased mineralisation rates (Meyerholt et al.2020b; Zaehle and Dalmonech2011). The response of the N cycle in JULES under changes in climate and CO2 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 (last access: 9 April 2021). To access the code a freely available non-commercial research licence is required (, last access: 9 April 2021). The suites required for running JULES are available here: (last access: 9 April 2021). JULES-CN is u-ah896, JULES-C is u-ah932, and JULES-CNlayer is u-ai571.

Data availability

JULES output data are available at (Wiltshire et al.2020).

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.


Arora, V. K., Katavouta, A., Williams, R. G., Jones, C. D., Brovkin, V., Friedlingstein, P., Schwinger, J., Bopp, L., Boucher, O., Cadule, P., Chamberlain, M. A., Christian, J. R., Delire, C., Fisher, R. A., Hajima, T., Ilyina, T., Joetzjer, E., Kawamiya, M., Koven, C. D., Krasting, J. P., Law, R. M., Lawrence, D. M., Lenton, A., Lindsay, K., Pongratz, J., Raddatz, T., Séférian, R., Tachiiri, K., Tjiputra, J. F., Wiltshire, A., Wu, T., and Ziehn, T.: Carbon–concentration and carbon–climate feedbacks in CMIP6 models and their comparison to CMIP5 models, Biogeosciences, 17, 4173–4222,, 2020. a, b, c

Avitabile, V., Herold, M., Heuvelink, G. B., Lewis, S. L., Phillips, O. L., Asner, G. P., Armston, J., Ashton, P. S., Banin, L., Bayol, N., Berry, N. J., Boeckx, P., de Jong, B. H. J., DeVries, B., Girardin, C. A. J., Kearsley, E., Lindsell, J. A., Lopez-Gonzalez, G., Lucas, R., Malhi, Y., Morel, A., Mitchard, E. T. A., Nagy, L., Qie, L., Quinones, M. J., Ryan, C. M., Ferry, S. J. W., Sunderland, T., Laurin, G. V., Gatti, R. C., Valentini, R., Verbeeck, H., Wijaya, A., and Willcock, S.: An integrated pan-tropical biomass map using multiple reference datasets, Global Change Biol., 22, 1406–1420, 2016. a

Ballantyne, A., Smith, W., Anderegg, W., Kauppi, P., Sarmiento, J., Tans, P., Shevliakova, E., Pan, Y., Poulter, B., Anav, A., Friedlingstein, P., Houghton, R., and Running, R.: Accelerating net terrestrial carbon uptake during the warming hiatus due to reduced respiration, Nat. Clim. Change, 7, 148–152, 2017. a

Batjes, N.: Total carbon and nitrogen in the soils of the world, European J. Soil Sci., 65, 10–21, 2014. a, b

Batjes, N. H.: Harmonized soil property values for broad-scale modelling (WISE30sec) with estimates of global soil carbon stocks, Geoderma, 269, 61–68, 2016. a, b

Best, M. J., Pryor, M., Clark, D. B., Rooney, G. G., Essery, R. L. H., Ménard, C. B., Edwards, J. M., Hendry, M. A., Porson, A., Gedney, N., Mercado, L. M., Sitch, S., Blyth, E., Boucher, O., Cox, P. M., Grimmond, C. S. B., and Harding, R. J.: The Joint UK Land Environment Simulator (JULES), model description – Part 1: Energy and water fluxes, Geosci. Model Dev., 4, 677–699,, 2011. a, b

Boyer, E. W., Howarth, R. W., Galloway, J. N., Dentener, F. J., Green, P. A., and Vörösmarty, C. J.: Riverine nitrogen export from the continents to the coasts, Global Biogeochem. Cy., 20, GB1S91,, 2006. a, b

Brix, H. and Ebell, L.: Effects of nitrogen fertilization on growth, leaf area, and photosynthesis rate in Douglas-fir, Forest Sci., 15, 189–196, 1969. a, b

Burke, E. J., Chadburn, S. E., and Ekici, A.: A vertical representation of soil carbon in the JULES land surface scheme (vn4.3_permafrost) with a focus on permafrost regions, Geosci. Model Dev., 10, 959–975,, 2017. a, b, c, d, e, f, g, h

Burton, C., Betts, R., Cardoso, M., Feldpausch, T. R., Harper, A., Jones, C. D., Kelley, D. I., Robertson, E., and Wiltshire, A.: Representation of fire, land-use change and vegetation dynamics in the Joint UK Land Environment Simulator vn4.9 (JULES), Geosci. Model Dev., 12, 179–193,, 2019. a, b

Carswell, F., Meir, P., Wandelli, E., Bonates, L., Kruijt, B., Barbosa, E., Nobre, A., Grace, J., and Jarvis, P.: Photosynthetic capacity in a central Amazonian rain forest, Tree Physiol., 20, 179–186, 2000. a

Carvalhais, N., Forkel, M., Khomik, M., Bellarby, J., Jung, M., Migliavacca, M., Mu, M., Saatchi, S., Santoro, M., Thurner, M., Weber, U., Ahrens, B., Beer, C., Cescatti, A., Randerson, J. T., and Reichstein, M.: Global covariation of carbon turnover times with climate in terrestrial ecosystems, Nature, 514, 213–217, 2014. a, b, c, d, e

Chadburn, S., Burke, E., Essery, R., Boike, J., Langer, M., Heikenfeld, M., Cox, P., and Friedlingstein, P.: An improved representation of physical permafrost dynamics in the JULES land-surface model, Geosci. Model Dev., 8, 1493–1508,, 2015. a, b

Ciais, P., Sabine, C., Bala, G., Bopp, L., Brovkin, V., Canadell, J., Chhabra, A., DeFries, R., Galloway, J., Heimann, M., Jones, C., Le Quéré, C., Myneni, R. B., Piao, S., and Thornton, P.: Carbon and other biogeochemical cycles, in: Climate change 2013: the physical science basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, Cambridge University Press, 465–570, 2014. a

Clark, D. B., Mercado, L. M., Sitch, S., Jones, C. D., Gedney, N., Best, M. J., Pryor, M., Rooney, G. G., Essery, R. L. H., Blyth, E., Boucher, O., Harding, R. J., Huntingford, C., and Cox, P. M.: The Joint UK Land Environment Simulator (JULES), model description – Part 2: Carbon fluxes and vegetation dynamics, Geosci. Model Dev., 4, 701–722,, 2011. a, b, c, d, e, f, g, h, i, j, k, l, m, n

Cleveland, C. C., Townsend, A. R., Schimel, D. S., Fisher, H., Howarth, R. W., Hedin, L. O., Perakis, S. S., Latty, E. F., Von Fischer, J. C., Elseroad, A., and Wasson, M.: Global patterns of terrestrial biological nitrogen (N2) fixation in natural ecosystems, Global Biogeochem. Cy., 13, 623–645, 1999. a, b

Collalti, A. and Prentice, I.: Is NPP proportional to GPP? Waring’s hypothesis 20 years on, Tree Physiol., 39, 1473–1483, 2019. a

Cox, P., Huntingford, C., and Harding, R.: A canopy conductance and photosynthesis model for use in a GCM land surface scheme, J. Hydrol., 212, 79–94, 1998. a

Cox, P., Betts, R., Bunton, C., Essery, R., Rowntree, P., and Smith, J.: The impact of new land surface physics on the GCM simulation of climate and climate sensitivity, Clim. Dynam., 15, 183–203, 1999. a

Cox, P. M.: Description of the TRIFFID dynamic global vegetation model, Tech. rep., Met Office Hadley Centre, Exeter, UK, available at: (last access: 9 April 2021), 2001. a, b, c, d, e

Cox, P. M., Pearson, D., Booth, B. B., Friedlingstein, P., Huntingford, C., Jones, C. D., and Luke, C. M.: Sensitivity of tropical carbon to climate change constrained by carbon dioxide variability, Nature, 494, 341–344, 2013. a, b

Davies-Barnard, T. and Friedlingstein, P.: The Global Distribution of Biological Nitrogen Fixation in Terrestrial Natural Ecosystems, Global Biogeochem. Cy., 34, e2019GB006387,, 2020. a, b, c, d

Davies-Barnard, T., Meyerholt, J., Zaehle, S., Friedlingstein, P., Brovkin, V., Fan, Y., Fisher, R. A., Jones, C. D., Lee, H., Peano, D., Smith, B., Wårlind, D., and Wiltshire, A. J.: Nitrogen cycling in CMIP6 land surface models: progress and limitations, Biogeosciences, 17, 5129–5148,, 2020. a, b, c, d, e, f, g

Dise, N. B., Rothwell, J. J., Gauci, V., van der Salm, C., and de Vries, W.: Predicting dissolved inorganic nitrogen leaching in European forests using two independent databases, Sci. Total Environ., 407, 1798–1808,, 2009. a

Enquist, B. J., Brown, J. H., and West, G. B.: Allometric scaling of plant energetics and population density, Nature, 395, 163–165, 1998. a

Field, C. and Mooney, H.: The photosynthesis-nitrogen relationship in wild plants, Cambridge University Press, 25–55, 1986. a

Fisher, J. B., Sitch, S., Malhi, Y., Fisher, R. A., Huntingford, C., and Tan, S.-Y.: Carbon cost of plant nitrogen acquisition: A mechanistic, globally applicable model of plant nitrogen uptake, retranslocation, and fixation, Global Biogeochem. Cy., 24, GB1014,, 2010. a

Friedlingstein, P., O'Sullivan, M., Jones, M. W., Andrew, R. M., Hauck, J., Olsen, A., Peters, G. P., Peters, W., Pongratz, J., Sitch, S., Le Quéré, C., Canadell, J. G., Ciais, P., Jackson, R. B., Alin, S., Aragão, L. E. O. C., Arneth, A., Arora, V., Bates, N. R., Becker, M., Benoit-Cattin, A., Bittig, H. C., Bopp, L., Bultan, S., Chandra, N., Chevallier, F., Chini, L. P., Evans, W., Florentie, L., Forster, P. M., Gasser, T., Gehlen, M., Gilfillan, D., Gkritzalis, T., Gregor, L., Gruber, N., Harris, I., Hartung, K., Haverd, V., Houghton, R. A., Ilyina, T., Jain, A. K., Joetzjer, E., Kadono, K., Kato, E., Kitidis, V., Korsbakken, J. I., Landschützer, P., Lefèvre, N., Lenton, A., Lienert, S., Liu, Z., Lombardozzi, D., Marland, G., Metzl, N., Munro, D. R., Nabel, J. E. M. S., Nakaoka, S.-I., Niwa, Y., O'Brien, K., Ono, T., Palmer, P. I., Pierrot, D., Poulter, B., Resplandy, L., Robertson, E., Rödenbeck, C., Schwinger, J., Séférian, R., Skjelvan, I., Smith, A. J. P., Sutton, A. J., Tanhua, T., Tans, P. P., Tian, H., Tilbrook, B., van der Werf, G., Vuichard, N., Walker, A. P., Wanninkhof, R., Watson, A. J., Willis, D., Wiltshire, A. J., Yuan, W., Yue, X., and Zaehle, S.: Global Carbon Budget 2020, Earth Syst. Sci. Data, 12, 3269–3340,, 2020. a, b, c, d, e

Galloway, J. N., Dentener, F. J., Capone, D. G., Boyer, E. W., Howarth, R. W., Seitzinger, S. P., Asner, G. P., Cleveland, C. C., Green, P., Holland, E. A., Karl, D. M., Michaels, A. F., Porter, J. H., Townsend, A. R., and Vöosmarty, C. J.: Nitrogen cycles: past, present, and future, Biogeochemistry, 70, 153–226, 2004. a, b

Global Soil Data Task Group: Global Gridded Surfaces of Selected Soil Characteristics (IGBP-DIS), ORNL DAAC, Oak Ridge, Tennessee, USA,, 2000. a, b, c

Gruber, N. and Galloway, J. N.: An Earth-system perspective of the global nitrogen cycle, Nature, 451, 293–296, 2008. a

Harper, A. B., Cox, P. M., Friedlingstein, P., Wiltshire, A. J., Jones, C. D., Sitch, S., Mercado, L. M., Groenendijk, M., Robertson, E., Kattge, J., Bönisch, G., Atkin, O. K., Bahn, M., Cornelissen, J., Niinemets, Ü., Onipchenko, V., Peñuelas, J., Poorter, L., Reich, P. B., Soudzilovskaia, N. A., and Bodegom, P. V.: Improved representation of plant functional types and physiology in the Joint UK Land Environment Simulator (JULES v4.2) using plant trait information, Geosci. Model Dev., 9, 2415–2440,, 2016. a, b

Harper, A. B., Wiltshire, A. J., Cox, P. M., Friedlingstein, P., Jones, C. D., Mercado, L. M., Sitch, S., Williams, K., and Duran-Rojas, C.: Vegetation distribution and terrestrial carbon cycle in a carbon cycle configuration of JULES4.6 with new plant functional types, Geosci. Model Dev., 11, 2857–2873,, 2018. a, b, c, d

Hartley, A., MacBean, N., Georgievski, G., and Bontemps, S.: Uncertainty in plant functional type distributions and its impact on land surface models, Remote Sens. Environ., 203, 71–89, 2017. a, b

Hashimoto, S., Carvalhais, N., Ito, A., Migliavacca, M., Nishina, K., and Reichstein, M.: Global spatiotemporal distribution of soil respiration modeled using a global database, Biogeosciences, 12, 4121–4132,, 2015. a, b

Houlton, B. Z., Wang, Y.-P., Vitousek, P. M., and Field, C. B.: A unifying framework for dinitrogen fixation in the terrestrial biosphere, Nature, 454, 327–330, 2008. a

Hugelius, G., Loisel, J., Chadburn, S., Jackson, R. B., Jones, M., MacDonald, G., Marushchak, M., Olefeldt, D., Packalen, M., Siewert, M. B., Treat, C., Turetsky, M., Voigt, C., and Yu, Z.: Large stocks of peatland carbon and nitrogen are vulnerable to permafrost thaw, P. Natl. Acad. Sci. USA, 117, 20438–20446, 2020. a

Hurtt, G. C., Chini, L. P., Frolking, S., Betts, R., Feddema, J., Fischer, G., Fisk, J., Hibbard, K., Houghton, R., Janetos, A., Jones, C. D., Kindermann, G., Kinoshita, T., Klein Goldewijk, K., Riahi, K., Shevliakova, E., Smith, S., Stehfest, E., Thomson, A., Thornton, P., van Vuuren, D. P., and Wang, Y. P.: Harmonization of land-use scenarios for the period 1500–2100: 600 years of global gridded annual land-use transitions, wood harvest, and resulting secondary lands, Clim. Change, 109, 117–161, 2011. a

Jenkinson, D. and Coleman, K.: A model for the turnover of carbon in soil—Model description and windows user guide, Rothamsted Research, Harpenden, UK, 1999. a

Jenkinson, D. S., Andrew, S. P. S., Lynch, J. M., Goss, M. J., and Tinker, P. B.: The Turnover of Organic Carbon and Nitrogen in Soil [and Discussion], Philos. T. R. Soc. B, 329, 361–368, (last access: 9 April 2021), 1990. a

Jones, C. D. and Friedlingstein, P.: Quantifying process-level uncertainty contributions to TCRE and Carbon Budgets for meeting Paris Agreement climate targets, Environ. Res. Lett., 15, 074019,, 2020. a

Jung, M., Reichstein, M., Margolis, H. A., Cescatti, A., Richardson, A. D., Arain, M. A., Arneth, A., Bernhofer, C., Bonal, D., Chen, J., Gianelle, D., Gobron, N., Kiely, G., Kutsch, W., Lasslop, G., Law, B. E., Lindroth, A., Merbold, L., Montagnani, L., Moors, E. J., Papale, D., Sottocornola, M., Vaccari, F., and Williams, C.: Global patterns of land-atmosphere fluxes of carbon dioxide, latent heat, and sensible heat derived from eddy covariance, satellite, and meteorological observations, J. Geophys. Res.-Biogeo., 116, G00J07,, g00J07, 2011. a, b

Kim, D., Lee, M.-I., Jeong, S.-J., Im, J., Cha, D. H., and Lee, S.: Intercomparison of Terrestrial Carbon Fluxes and Carbon Use Efficiency Simulated by CMIP5 Earth System Models, Asia-Pac. J. Atmos. Sci., 54, 145–163, 2018. a, b

Koven, C. D., Riley, W. J., Subin, Z. M., Tang, J. Y., Torn, M. S., Collins, W. D., Bonan, G. B., Lawrence, D. M., and Swenson, S. C.: The effect of vertically resolved soil biogeochemistry and alternate soil C and N models on C dynamics of CLM4, Biogeosciences, 10, 7109–7131,, 2013. a

Lamarque, J.-F., Shindell, D. T., Josse, B., Young, P. J., Cionni, I., Eyring, V., Bergmann, D., Cameron-Smith, P., Collins, W. J., Doherty, R., Dalsoren, S., Faluvegi, G., Folberth, G., Ghan, S. J., Horowitz, L. W., Lee, Y. H., MacKenzie, I. A., Nagashima, T., Naik, V., Plummer, D., Righi, M., Rumbold, S. T., Schulz, M., Skeie, R. B., Stevenson, D. S., Strode, S., Sudo, K., Szopa, S., Voulgarakis, A., and Zeng, G.: The Atmospheric Chemistry and Climate Model Intercomparison Project (ACCMIP): overview and description of models, simulations and climate diagnostics, Geosci. Model Dev., 6, 179–206,, 2013. a

Le Quéré, C., Andrew, R. M., Friedlingstein, P., Sitch, S., Hauck, J., Pongratz, J., Pickers, P. A., Korsbakken, J. I., Peters, G. P., Canadell, J. G., Arneth, A., Arora, V. K., Barbero, L., Bastos, A., Bopp, L., Chevallier, F., Chini, L. P., Ciais, P., Doney, S. C., Gkritzalis, T., Goll, D. S., Harris, I., Haverd, V., Hoffman, F. M., Hoppema, M., Houghton, R. A., Hurtt, G., Ilyina, T., Jain, A. K., Johannessen, T., Jones, C. D., Kato, E., Keeling, R. F., Goldewijk, K. K., Landschützer, P., Lefèvre, N., Lienert, S., Liu, Z., Lombardozzi, D., Metzl, N., Munro, D. R., Nabel, J. E. M. S., Nakaoka, S., Neill, C., Olsen, A., Ono, T., Patra, P., Peregon, A., Peters, W., Peylin, P., Pfeil, B., Pierrot, D., Poulter, B., Rehder, G., Resplandy, L., Robertson, E., Rocher, M., Rödenbeck, C., Schuster, U., Schwinger, J., Séférian, R., Skjelvan, I., Steinhoff, T., Sutton, A., Tans, P. P., Tian, H., Tilbrook, B., Tubiello, F. N., van der Laan-Luijkx, I. T., van der Werf, G. R., Viovy, N., Walker, A. P., Wiltshire, A. J., Wright, R., Zaehle, S., and Zheng, B.: Global Carbon Budget 2018, Earth Syst. Sci. Data, 10, 2141–2194,, 2018. a, b, c

LeBauer, D. S. and Treseder, K. K.: Nitrogen limitation of net primary productivity in terrestrial ecosystems is globally distributed, Ecology, 89, 371–379,, 2008. a, b

Li, W., Ciais, P., Wang, Y., Yin, Y., Peng, S., Zhu, Z., Bastos, A., Yue, C., Ballantyne, A. P., Broquet, G., Canadell, J. G., Cescatti, A., Chen, C., Cooper, L., Friedlingstein, P., Le Quéré, C., Myneni, R. B., and Piao, S.: Recent changes in global photosynthesis and terrestrial ecosystem respiration constrained from multiple observations, Geophys. Res. Lett., 45, 1058—1068, 2018. a

Liang, J., Qi, X., Souza, L., and Luo, Y.: Processes regulating progressive nitrogen limitation under elevated carbon dioxide: a meta-analysis, Biogeosciences, 13, 2689–2699,, 2016. a

Mao, J., Mao, Q., Zheng, M., and Mo, J.: Responses of Foliar Nutrient Status and Stoichiometry to Nitrogen Addition in Different Ecosystems: A Meta-analysis, J. Geophys. Res.-Biogeo., 125, e2019JG005347,, 2020. a, b

McGuire, A. D., Melillo, J. M., and Joyce, L. A.: The role of nitrogen in the response of forest net primary production to elevated atmospheric carbon dioxide, Ann. Rev. Ecol. Sys., 26, 473–503, 1995. a

Mercado, L. M., Huntingford, C., Gash, J. H., Cox, P. M., and Jogireddy, V.: Improving the representation of radiation interception and photosynthesis for climate model applications, Tellus B, 59, 553–565, 2007. a, b

Mercado, L. M., Bellouin, N., Sitch, S., Boucher, O., Huntingford, C., Wild, M., and Cox, P. M.: Impact of changes in diffuse radiation on the global land carbon sink, Nature, 458, 1014–1017, 2009. a

Meyerholt, J. and Zaehle, S.: The role of stoichiometric flexibility in modelling forest ecosystem responses to nitrogen fertilization, New Phytol., 208, 1042–1055, 2015. a, b, c

Meyerholt, J., Zaehle, S., and Smith, M. J.: Variability of projected terrestrial biosphere responses to elevated levels of atmospheric CO2 due to uncertainty in biological nitrogen fixation, Biogeosciences, 13, 1491–1518,, 2016. a

Meyerholt, J., Sickel, K., and Zaehle, S.: Ensemble projections elucidate effects of uncertainty in terrestrial nitrogen limitation on future carbon uptake, Glob. Change Biol., 26, 3978–3996,, 2020a. a, b

Meyerholt, J., Sickel, K., and Zaehle, S.: Ensemble projections elucidate effects of uncertainty in terrestrial nitrogen limitation on future carbon uptake, Glob. Change Biol., 26, 3978–3996, 2020b. a

Ochoa-Hueso, R., Maestre, F. T., de los Ríos, A., Valea, S., Theobald, M. R., Vivanco, M. G., Manrique, E., and Bowker, M. A.: Nitrogen deposition alters nitrogen cycling and reduces soil carbon content in low-productivity semiarid Mediterranean ecosystems, Environ. Pollu., 179, 185–193,, 2013. a

Olson, D. M., Dinerstein, E., Wikramanayake, E. D., Burgess, N. D., Powell, G. V., Underwood, E. C., D'amico, J. A., Itoua, I., Strand, H. E., Morrison, J. C., Loucks, C. J., Allnutt, T. F., Ricketts, T. H., Kura, Y., Lamoreux, J. F., Wettengel, W. W., Hedao, P., and Kassem, K. R.: Terrestrial Ecoregions of the World: A New Map of Life on Earth A new global map of terrestrial ecoregions provides an innovative tool for conserving biodiversity, BioScience, 51, 933–938, 2001. a, b

Potter, P., Ramankutty, N., Bennett, E. M., and Donner, S. D.: Characterizing the spatial patterns of global fertilizer application and manure production, Earth Interactions, 14, 1–22, 2010. a

Poulter, B., MacBean, N., Hartley, A., Khlystova, I., Arino, O., Betts, R., Bontemps, S., Boettcher, M., Brockmann, C., Defourny, P., Hagemann, S., Herold, M., Kirches, G., Lamarche, C., Lederer, D., Ottlé, C., Peters, M., and Peylin, P.: Plant functional type classification for earth system models: results from the European Space Agency's Land Cover Climate Change Initiative, Geosci. Model Dev., 8, 2315–2328,, 2015. a

Ruesch, A. and Gibb, H. K.: New IPCC Tier-1 Global Biomass Carbon Map For the Year 2000, available: (last access: 9 April 2021), 2008. a

Saatchi, S. S., Harris, N. L., Brown, S., Lefsky, M., Mitchard, E. T., Salas, W., Zutta, B. R., Buermann, W., Lewis, S. L., Hagen, S., Petrova, S., White, L., Silman, M., and Morel, A.: Benchmark map of forest carbon stocks in tropical regions across three continents, P. Natl. Acad. Sci. USA, 108, 9899–9904, 2011. a

Schlesinger, W. H.: Biogeochemistry: An Analysis of Global Change, Academic Press, San Diego, 2nd Edn., 1997. a

Sellar, A. A., Jones, C. G., Mulcahy, J. P., Tang, Y., Yool, A., Wiltshire, A., O'Connor, F. M., Stringer, M., Hill, R., Palmieri, J., Woodward, S., de Mora, L., Kuhlbrodt, T., Rumbold, S. T., Kelley, D. I., Ellis, R., Johnson, C. E., Walton, J., Abraham, N. L., Andrews, M. B., Andrews, T., Archibald, A. T., Berthou, S., Burke, E., Blockley, E., Carslaw, K., Dalvi, M., Edwards, J., Folberth, G. A., Gedney, N., Griffiths, P. T., Harper, A. B., Hendry, M. A., Hewitt, A. J., Johnson, B., Jones, A., Jones, C. D., Keeble, J., Liddicoat, S., Morgenstern, O., Parker, R. J., Predoi, V., Robertson, E., Siahaan, A., Smith, R. S., Swaminathan, R., Woodhouse, M. T., Zeng, G., and Zerroukat, M.: UKESM1: Description and Evaluation of the U.K. Earth System Model, J. Adv. Model. Ea. Sy., 11, 4513–4558,, 2019. a, b, c, d, e, f, g

Shinozaki, K., Yoda, K., Hozumi, K., and Kira, T.: A quantitative analysis of plant form-the pipe model theory: I. Basic analyses, Jpn. J. Ecol., 14, 97–105, 1964a. a

Shinozaki, K., Yoda, K., Hozumi, K., and Kira, T.: A quantitative analysis of plant form-the pipe model theory: II. Further evidence of the theory and its application in forest ecology, Jpn. J. Ecol., 14, 133–139, 1964b. a

Sitch, S., Friedlingstein, P., Gruber, N., Jones, S. D., Murray-Tortarolo, G., Ahlström, A., Doney, S. C., Graven, H., Heinze, C., Huntingford, C., Levis, S., Levy, P. E., Lomas, M., Poulter, B., Viovy, N., Zaehle, S., Zeng, N., Arneth, A., Bonan, G., Bopp, L., Canadell, J. G., Chevallier, F., Ciais, P., Ellis, R., Gloor, M., Peylin, P., Piao, S. L., Le Quéré, C., Smith, B., Zhu, Z., and Myneni, R.: Recent trends and drivers of regional sources and sinks of carbon dioxide, Biogeosciences, 12, 653–679,, 2015. a

Smith, B., Warlind, D., Arneth, A., Hickler, T., Leadley, P., Siltberg, J., and Zaehle, S.: Implications of incorporating N cycling and N limitations on primary production in an individual-based dynamic vegetation model, Biogeosciences, 11, 2027–2054, 2014. a, b

Sterner, R. W. and Elser, J. J.: Ecological stoichiometry: the biology of elements from molecules to the biosphere, Princeton university press, 2002. a

Sullivan, B. W., Smith, W. K., Townsend, A. R., Nasto, M. K., Reed, S. C., Chazdon, R. L., and Cleveland, C. C.: Spatially robust estimates of biological nitrogen (N) fixation imply substantial human alteration of the tropical N cycle, P. Natl. Acad. Sci., 111, 8101–8106,, 2014. a

Thomas, R. Q., Bonan, G. B., and Goodale, C. L.: Insights into mechanisms governing forest carbon response to nitrogen deposition: a model–data comparison using observed responses to nitrogen addition, Biogeosciences, 10, 3869–3887,, 2013a. a

Thomas, R. Q., Zaehle, S., Templer, P. H., and Goodale, C. L.: Global patterns of nitrogen limitation: confronting two global biogeochemical models with observations, Glob. Change Biol., 19, 2986–2998,, 2013b. a, b

Viovy, N.: CRUNCEP Version 7 – Atmospheric Forcing Data for the Community Land Model, Research Data Archive at the National Center for Atmospheric Research, Computational and Information Systems Laboratory., 2018. a

Vitousek, P. M., Menge, D. N., Reed, S. C., and Cleveland, C. C.: Biological nitrogen fixation: rates, patterns and ecological controls in terrestrial ecosystems, Philos. T. Roy. Soc. B, 368, 20130119,, 2013. a, b, c

von Bloh, W., Schaphoff, S., Müller, C., Rolinski, S., Waha, K., and Zaehle, S.: Implementing the nitrogen cycle into the dynamic global vegetation, hydrology, and crop growth model LPJmL (version 5.0), Geosci. Model Dev., 11, 2789–2812,, 2018. a, b

Wang, D., Maughan, M. W., Sun, J., Feng, X., Miguez, F., Lee, D., and Dietze, M. C.: Impact of nitrogen allocation on growth and photosynthesis of Miscanthus (Miscanthus × giganteus), Gcb Bioenergy, 4, 688–697, 2012. a, b

Wania, R., Meissner, K. J., Eby, M., Arora, V. K., Ross, I., and Weaver, A. J.: Carbon-nitrogen feedbacks in the UVic ESCM, Geosci. Model Dev., 5, 1137–1160,, 2012. a, b, c

Weintraub, M. N. and Schimel, J. P.: Nitrogen Cycling and the Spread of Shrubs Control Changes in the Carbon Balance of Arctic Tundra Ecosystems, BioScience, 55, 408–415,[0408:NCATSO]2.0.CO;2, 2005. a

Wenzel, S., Cox, P. M., Eyring, V., and Friedlingstein, P.: Projected land photosynthesis constrained by changes in the seasonal cycle of atmospheric CO2, Nature, 538, 499–501, 2016. a, b

Wieder, W. R., Cleveland, C. C., Lawrence, D. M., and Bonan, G. B.: Effects of model structural uncertainty on carbon cycle projections: biological nitrogen fixation as a case study, Environ. Res. Lett., 10, 044016,, 2015a. a

Wieder, W. R., Cleveland, C. C., Smith, W. K., and Todd-Brown, K.: Future productivity and carbon storage limited by terrestrial nutrient availability, Nat. Geosci., 8, 441–444, 2015b. a

Wiltshire, A. J., Duran Rojas, M. C., Edwards, J. M., Gedney, N., Harper, A. B., Hartley, A. J., Hendry, M. A., Robertson, E., and Smout-Day, K.: JULES-GL7: the Global Land configuration of the Joint UK Land Environment Simulator version 7.0 and 7.2, Geosci. Model Dev., 13, 483–505,, 2020. a, b

Wiltshire, A., Burke, E., Chadburn, S., Jones, C., Cox, P., Davies-Barnard, T., Friedlingstein, P., Harper, A., Liddicoat, S., Sitch, S., and Zaehle, S.: JULES-CN: a coupled terrestrial Carbon-Nitrogen Scheme (JULES vn5.1) [Data set], Zenodo,, 2021. a

Xu-Ri and Prentice, I. C.: Terrestrial nitrogen cycle simulation with a dynamic global vegetation model, Glob. Change Biol., 14, 1745–1764,, 2008. a, b

Yang, X., Thornton, P. E., Ricciuto, D. M., and Post, W. M.: The role of phosphorus dynamics in tropical forests – a modeling study using CLM-CNP, Biogeosciences, 11, 1667–1681,, 2014. a

Zaehle, S.: Terrestrial nitrogen–carbon cycle interactions at the global scale, Philos. T. Roy. Soc. B, 368, 20130125,, 2013. a, b, c

Zaehle, S. and Dalmonech, D.: Carbon–nitrogen interactions on land at global scales: current understanding in modelling climate biosphere feedbacks, Curr. Opin. Env. Sust., 3, 311–320, 2011. a, b

Zaehle, S. and Friend, A.: Carbon and nitrogen cycle dynamics in the O-CN land surface model: 1. Model description, site-scale evaluation, and sensitivity to parameter estimates, Global Biogeochem. Cy., 24, GB1005,, 2010. a, b

Zaehle, S., Friend, A., Friedlingstein, P., Dentener, F., Peylin, P., and Schulz, M.: Carbon and nitrogen cycle dynamics in the O-CN land surface model: 2. Role of the nitrogen cycle in the historical terrestrial carbon balance, Global Biogeochem. Cy., 24, GB1006,, 2010. a

Zaehle, S., Medlyn, B. E., De Kauwe, M. G., Walker, A. P., Dietze, M. C., Hickler, T., Luo, Y., Wang, Y.-P., El-Masri, B., Thornton, P., Jain, A., Wang, S., Warlind, D., Weng, E., Parton, W., Iversen, C. M., Gallet-Budynek, A., McCarthy, H., Finzi, A., Hanson, P. J., Prentice, I. C., Oren, R., and Norby, R. J.: Evaluation of 11 terrestrial carbon–nitrogen cycle models against observations from two temperate Free-Air CO2 Enrichment studies, New Phytol., 202, 803–822,, 2013-16358, 2014.  a

Zaehle, S., Jones, C. D., Houlton, B., Lamarque, J.-F., and Robertson, E.: Nitrogen availability reduces CMIP5 projections of twenty-first-century land carbon uptake, J. Climate, 28, 2494–2511, 2015. a

Zhao, M. and Running, S. W.: Drought-induced reduction in global terrestrial net primary production from 2000 through 2009, Science, 329, 940–943, 2010. a, b

Zheng, M., Zhou, Z., Zhao, P., Luo, Y., Ye, Q., Zhang, K., Song, L., and Mo, J.: Effects of human disturbance activities and environmental change factors on terrestrial nitrogen fixation, Glob. Change Biol., 26, 6203–6217,, 2020. a

Short summary
Limited nitrogen availbility can restrict the growth of plants and their ability to assimilate carbon. It is important to include the impact of this process on the global land carbon cycle. This paper presents a model of the coupled land carbon and nitrogen cycle, which is included within the UK Earth System model to improve projections of climate change and impacts on ecosystems.