Thawing Permafrost as a Nitrogen Fertiliser: Implications for Climate Feedbacks

: Studies for the northern high latitudes suggest that, in the near term, increased vegetation uptake may offset permafrost carbon losses, but over longer time periods, permafrost carbon decomposition causes a net loss of carbon. Here, we assess the impact of a coupled carbon and nitrogen cycle on the simulations of these carbon ﬂuxes. We present results from JULES-IMOGEN—a global land surface model coupled to an intermediate complexity climate model with vertically resolved soil biogeochemistry. We quantify the impact of nitrogen fertilisation from thawing permafrost on the carbon cycle and compare it with the loss of permafrost carbon. Projections show that the additional fertilisation reduces the high latitude vegetation nitrogen limitation and causes an overall increase in vegetation carbon uptake. This is a few Petagrams of carbon (Pg C) by year 2100, increasing to up to 40 Pg C by year 2300 for the RCP8.5 concentration scenario and adds around 50% to the projected overall increase in vegetation carbon in that region. This nitrogen fertilisation results in a negative (stabilising) feedback on the global mean temperature, which could be equivalent in magnitude to the positive (destabilising) temperature feedback from the loss of permafrost carbon. This balance depends on the future scenario and initial permafrost carbon. JULES-IMOGEN describes one representation of the changes in Arctic carbon and nitrogen cycling in response to climate change. However there are uncertainties in the modelling framework, model parameterisation and missing processes which, when assessed, will provide a more complete picture of the balance between stabilising and destabilising feedbacks.


Introduction
Permafrost is defined as ground that has been frozen for a period of two years or more.It underlies about 22% of the northern hemisphere land surface [1,2].There is evidence that permafrost temperatures have been increasing in recent years [3,4] by as much as 0.39 ± 0.15 • C in the continuous permafrost zone (2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016).This is projected to continue to increase in the future, in response to increases in atmospheric greenhouse gas concentrations caused by the human burning of fossil fuels.For example, Ref. [5] project a loss of permafrost area of 4.0 +1.0 −1.1 × 10 6 km 2 / • C in conditions where the soil temperatures are in equilibrium with the air temperatures.This unprecedented thawing will have implications for both hydrological and biogeochemical cycles in the permafrost region, as noted by many recent research studies [6][7][8][9][10][11][12][13].
We present new model simulations showing the impact of permafrost thaw on the northern high latitude carbon terrestrial carbon (C) and nitrogen (N) cycle and how this feeds back onto global climate change.This work expands on the estimate of the positive (destabilising) permafrost carbon feedback by [6] who used the JULES land surface model with an interactive representation of the global C cycle coupled with the IMOGEN intermediate climate model.They projected that the additional warming caused by thawing of permafrost C contributed between 0.2 and 12% of the change in the global mean temperature (GSAT) by year 2100 and 0.5 and 17% by 2300.The latest modelling results presented here additionally quantify the impact on future climate change of including an explicit representation of an interactive C and N cycle within JULES [14].Here, we focus specifically on the increase in vegetation in response to additional nitrogen fertilisation from thawing permafrost-a negative (stabilising) feedback.Additional processes which may also feedback onto the climate are not discussed here but are highlighted in the Conclusions at the end of the paper.
Permafrost systems are typically nitrogen limited [15].Under future climate change, substantial permafrost thaw is expected (e.g., [5]), and this is expected to increase mineral N available to plants and potentially reduce their nitrogen limitation [16].This extra N availability could happen through two routes.The first is that newly thawed mineral N, which is currently locked within the permafrost [17][18][19], may become available.The second route is potentially via increased net mobilisation of the organic N in the deeper soils [16,20].Vegetation is expected to respond to this increase in nitrogen availability and thaw via deeper root profiles, increased growth and changes in species composition [16,21,22].This response will enhance the vegetation C uptake in the permafrost region and modify the overall permafrost C feedback [6,13].One modelling study by [23] used the Community Land Model (CLM) with a coupled carbon-nitrogen model that includes permafrost processes to show that under future climate change, increased soil N mineralization in the surface soil reduces plant nutrient limitation.CLM simulations also show that N deeper in the soil only has a small impact on the overall C budget.However, ref. [23] did not include the response of the mineral N pool currently locked within the permafrost to climate change [24].
Earth System Models (ESMs) are designed to project the future evolution of the Earth's climate in response, amongst others, to changes in atmospheric greenhouse gas concentrations.In addition to modelling the physical climate, ESMs also include climate interactions with biogeochemical cycles.Although an interactive N cycle is now included within many of the ESMs participating in the Sixth Phase of the Coupled Model Intercomparison Project (CMIP6), very few include permafrost C. Hence, ESMs are unable to capture the interactions between permafrost thaw and the nitrogen cycle.This paper adopts an intermediate complexity climate model tuned to represent a range of climate sensitivities and regional climate change features.This climate model is coupled with the JULES land surface model to quantify the impact of an interactive C and N cycle in the permafrost region.Specifically we explore the vegetation C response under a warming climate to increased available mineral N in the soil (minNfert) and its feedback onto global climate (minNfert-feed).We compare this response with the permafrost C loss (pfCloss) and its feedback onto global climate (pfCloss-feed) in the presence of an interactive C and N cycle, to determine the balance of their impact on the global carbon cycle.

JULES Land Surface Scheme
The Joint UK Land Environment Simulator (JULES) forms the land-surface component of the UK Earth System Model (UKESM) [25][26][27].JULES simulates fluxes of energy, water, carbon and more recently nitrogen [14].It has a dynamic vegetation model, namely, the Top-down Representation of Interactive Foliage and Flora Including Dynamics (TRIFFID) where different plant functional types (PFTs) compete for space based on their height [28].Recent developments within JULES include the addition of a terrestrial nitrogen cycle which has a representation of all of the key terrestrial N processes and that integrate tightly with the structure of the C cycle [14].Inputs of N to the land surface are via biological fixation, fertilisation and N deposition, while losses of N from the land surface occur via leaching and gas loss.JULES simulates a N-limited ecosystem by reducing the net primary productivity if there is insufficient available N to accommodate the demand from the vegetation.
Ref. [14] describe in detail the vertically resolved soil biogeochemistry which is the focus of this paper.Below we briefly describe the model and highlight the relevant processes when assessing the permafrost-nitrogen feedbacks.A key requirement is to quantify the amount of inorganic nitrogen in the soil.In addition, the model needs to be able to divide the soil inorganic nitrogen pool into the fraction available for the plants to take up and the fraction not available either because there are no roots present or because the soil is frozen.
The soil biogeochemistry is represented using a vertically resolved version of the standard four-pool RothC soil carbon model [14,29].In each soil layer there is a soil carbon pool, that in turn disaggregates in to decomposable plant material (DPM), resistant plant material (RPM), microbial biomass (BIO), and humus (HUM), plus an equivalent soil organic nitrogen pool.There is a vertical mixing term representing either bioturbation (i.e., the soil mixing by, for example, animals and plant roots), or, in permafrost regions, cryoturbation (soil mixing is from frost heave and freeze-thaw processes).Each soil layer also has a mineral (or inorganic) N pool.N transfers between the organic and mineral N pools depending on respiration losses and the C-to-N ratio of the organic pool.Changes to the mineral N pool and at a given depth, z (N in (z) in kg [N]m −2 ) are modelled as follows: where N in is the inorganic N in kg [N] m −2 , z is the depth in m and t is the time in seconds.
The first right-hand term in Equation ( 1) is the biological N fixation, where v i is the PFT fraction for PFT i.
is the overall availability of biological N fixation for each PFT.The inputs from biological N fixation for each PFT are distributed vertically according to the normalised fraction of roots in each soil layer ( f R,i (z), dimensionless), that modulates the availability via F i .The second term in Equation ( 1) is the plant N uptake where Φ i (kg [N] m −2 s −1 ) is the plant N uptake for PFT i, and f U,i (z) is the fraction of available inorganic N in each soil layer for each PFT (described in more detail in the next paragraph).This is dependent on both the rooting depth and the freeze/thaw status of the soil.Specifically, if the soil layer is frozen or there are no roots present, there is no plant available nitrogen in that layer.The further terms are ) represents external inputs via surface deposition and (when applied) fertilisation; ) is the transfer of dissolved inorganic N between layers and N gas (z) (kg [N]m −2 s −1 ) is the total gaseous emission.
In order to calculate the plant uptake we assume that plants cannot access all the inorganic N (N in (z))-it can only access the plant available nitrogen (N avail (z) in kg [N] m −2 ) which is dependent on soil depth (z): If a soil layer is frozen then there is no plant-available nitrogen in the layer and plants cannot uptake any of the N in that layer and T(z) is zero, otherwise T(z) is one.Secondly, we assume that the plants only have direct access to a certain fraction of the soil, according to their root fraction, f root,i (z) (which reduces with depth).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., where z max is the maximum depth of the soil in m.JULES calculates the N demand based on the N required for spreading and growth [14].If the plant N demand is less than the available nitrogen the plant is not N limited and all of the carbon can be used for growth and spreading.However, if the plant N demand is greater than the available nitrogen there is excess carbon that cannot be assimilated by the plants and is considered to be an additional plant respiration term.This is used to reduce the Net Primary Productivity from its non N-limited value to its actual value.
We assess the impact of this vertically resolved soil biogeochemistry on the terrestrial carbon stores and their response to climate change.We also compare the impact on the terrestrial carbon cycle with the permafrost carbon loss from thawing permafrost and its feedback onto global mean temperature (pfCloss and pfCloss-feed).In the paragraph below, we briefly summarise the methodology used to quantify the permafrost carbon loss and feedback and describe it in full in the Appendix A.
In JULES we separate the permafrost carbon, defined as carbon below the maximum summer thaw depth, from the non-permafrost carbon using a numerical tracer which labels the permafrost carbon [6,29].However, JULES simulates an on-going small exchange between the carbon above and that below the permafrost table caused both by mixing processes and slow simulated decomposition in the permafrost in response to global warming.Hence in the simulations and over very many thousands of modelled years, all the labelled permafrost carbon will eventually turn over and be replaced by non-labelled carbon from above the permafrost table.This is an acknowledged artefact of the modelling framework.In this paper, unlike previous assessments with JULES-IMOGEN [6,30], this is now quantified using a stable climate (Figure A1).The permafrost carbon lost (pfCloss) is the difference between this drift and the loss of permafrost carbon under the future scenarios (Figure A2).The permafrost carbon feedback temperature is calculated by comparing the GSAT for model simulations with and without the permafrost carbon loss included in the land carbon exchange with the atmosphere.

IMOGEN
The Integrated Model Of Global Effects of climatic aNomalies (IMOGEN) is an intermediate complexity climate model that can be fully coupled with any land surface model.In this case, IMOGEN supplies the meteorological data to drive JULES which simulates changes in the carbon stores of the land surface that feedback onto changes in the IMOGEN-simulated atmospheric CO 2 concentration [31].We denote this model JULES-IMOGEN.IMOGEN consists of a simple energy balance model (EBM; [32]) that relates changes in concentrations of different atmospheric greenhouse gases to the global mean land temperature response, via changes in their combined radiative forcing.IMOGEN uses a "pattern-scaling" approach [32] to then relate, linearly, the amount of annual global average warming over land (via the EBM) to monthly changes in local meteorology.This local meteorology is disaggregated to hourly timesteps, and then used to drive JULES.At the end of each year, the net annual global land-atmosphere carbon flux calculated by JULES is passed back to modulate calculated atmospheric CO 2 .The land-ocean CO 2 flux is calculated using a single "box" model, and is a function of global temperature increase and atmospheric CO 2 level [31].Hence IMOGEN characterises both land-and ocean-atmosphere climate-carbon cycle feedbacks.With these feedbacks switched on in the JULES-IMOGEN simulation framework, the global carbon cycle is therefore interactive and the model system is forced with anthropogenic CO 2 emissions.From the combination of emissions and feedbacks, IMOGEN determines evolving atmospheric CO 2 concentrations, which in turn then influence the EBM by altering atmospheric radiative forcing.Both the EBM and the spatial patterns of IMOGEN are calibrated to emulate the climate change patterns for 34 of the ESMs from CMIP5 described by [33,34].Here, we assess results from three of those ensemble members which were selected to span the range of climate.Hence, this modelling framework spans the full range of uncertainty in both overall climate sensitivity (i.e., thermal sensitivity to rising atmospheric GHGs), and a set of potential regional and seasonal distributions of climate change in response to background global warming.

Experimental Design
The standard version of JULES-IMOGEN was first "spun-up" using pre-industrial atmospheric CO 2 and repeating climatology of the Water and Global Change forcing data for the period of 1961-1990 [35].Although this period is not for pre-industrial times, it does contain three decades where there are a substantial number of weather stations available and with a large geographical spread.Hence, this is assumed to be representative of preindustrial conditions (for details see [6]).The initial simulations are designed to generate a stable time-invariant starting state, to initialise subsequent calculations with a changing climate.In particular, for the analysis here, there is an emphasis on ensuring that both the soil carbon and vegetation distributions are stable at all locations.
The global spun-up states were used to initialise an ensemble of transient simulations assumed to start in 1850 with pre-industrial atmospheric CO 2 , onto which anomalies calculated by IMOGEN (based on the climate change from the CMIP5 ESMs) were added.This ensemble describes the effects of historical climate change, followed by three future scenarios with changing atmospheric greenhouse gases.The historical ensemble was forced with known historical fossil fuel burning and cement production CO 2 emissions.These were then followed using the CO 2 emissions that correspond to three of the Representative Concentration Pathways (RCP) used in the fifth assessment report of the Intergovernmental Panel on Climate Change [6,36].These scenarios are the RCP2.6 high mitigation scenario, the RCP4.5 intermediate scenario and the RCP8.5 high emissions scenario.Simulations were carried out until year 2300 using the RCP extensions [37] to examine the long-term relationship between permafrost and climate.Neither non-CO 2 greenhouse gases, aerosols nor land use change emissions were included in this set of simulations.However, the impact of these extra emissions will be relatively minor for our analysis.Instead, we focus on the feedbacks caused by including the interactive carbon and nitrogen cycle within JULES, and in particular for regions of thawing permafrost.
The labelled permafrost C was initialised at the start of the simulation as the soil C that was below the top of the simulated permafrost table at that time.The standard carbon and nitrogen (CN) configuration (CN spun : Table 1) uses the amount and distribution of soil carbon that was calculated by JULES-IMOGEN at the end of the spin-up.There are two additional CN configurations (CN 1.5spun and CN 2spun ) where the labelled permafrost carbon was prescribed to be equivalent to 1.5 or 2 times (respectively) that in the spun up state of CN spun .These three configurations span the uncertainty in the observed permafrost carbon.In particular, they have different initial conditions for assessing the impact of uncertainty in the initial permafrost carbon content on the projected permafrost carbon loss under climate change (pfCloss).Table 1 summarises the initial permafrost C contents of the named CN simulations by JULES-IMOGEN.
Table 1.JULES-IMOGEN initialisation of the labelled permafrost carbon that is permanently frozen in the soil at the start of the simulation.Uncertainties in the initial labelled permafrost C pool are also represented, via 1.5 × and 2× the spun up labelled permafrost carbon (second and third rows).
Here the initial and labelled permafrost C pool excludes any carbon that is above the top of the permafrost table.

Simulation Name
Permafrost An ensemble of transient simulations was run for each of the named configurations in Table 1.In addition, a set of diagnostic experiments were developed to assess the impact of the interactive carbon and nitrogen cycle with the layered soil biogeochemistry on both the northern high latitudes and on global climate feedbacks (Table 2).For each diagnostic experiment an ensemble of simulations was carried out where each member is equivalent to a member of the standard transient simulations.• pfCloss is the difference in permafrost C (Gt C) between transient-ens and control-ens.

pfCloss-feed-ens Diagnostic experiment
Quantifying the positive permafrost C feedback onto GSAT.The diagnosed permafrost C emissions from control-ens are included in the CO 2 emissions file and only the non-permafrost C is visible to IMOGEN.
The following sections use the ensembles defined in Table 2 to quantify the impact of climate change on the northern high latitude C cycle in the presence of an interactive C and N cycle, which includes the mineral N pool locked within the present-day permafrost soils.

Model Evaluation
The region over which the CN spun model simulates permafrost, for the period 2000-2009, is shown in Figure 1a.JULES-IMOGEN simulations are at a relatively coarse resolution and have no representation of sub-grid scale variability, so they only show either presence or absence of permafrost in any particular grid cell.The observations ( [1]; our Figure 1b) are provided at a much a higher resolution and give the proportion of permafrost in each grid cell.In the observations, most of the regions where the annual mean air temperature is less than 0 • C (MAAT area < 0 • C) have a non-zero probability of permafrost.We would prefer the model to represent permafrost in all of the cells where the observations have a probability of permafrost presence greater than 50%.Figure 1 suggests that the model over-estimates the permafrost extent in eastern Russia and underestimates it in eastern Canada.Table 3 compares the modelled permafrost extent of 15.2 million km 2 with the observed area of permafrost of 14.3 and 16.2 million km 2 [1,38].(These observed areas are derived from the area of each grid cell weighted by the proportion of the grid cell containing permafrost.)Another relevant metric shown in Table 3 is the permafrost extent as a percentage of the land surface area where MAAT area < 0 • C.This is again comparable between the model and the observations, with the modelled value (67%) falling between the two observed values of 62% to 71%.Table 3. Evaluation metrics for CN spun using the transient-ens, based on model simulations of permafrost extent, and vegetation and soil attributes.The carbon and nitrogen cycle metrics ("Area of trees" and below) are derived only for the grid cells where the model simulates permafrost for the period 2000-2009.Soil organic C and inorganic N contents are defined over the top 3 m and the soil organic N is defined for the top 1 m of soil * .GPP stands for gross primary productivity, NPP for net primary productivity and NEE for net ecosystem exchange with positive reflecting net uptake at the land surface.Ecosystem turnover is defined as (vegetation C + soil C)/GPP and soil turnover is defined as soil C/NPP.The response ratio a is the ratio of the potential amount of C that can be allocated to growth and spreading of the vegetation when the nitrogen is freely available compared with the actual amount achieved in the natural state when N may be limited.It is calculated from the multi-annual means of all the grid cells in the respective biomes.The boreal and tundra biomes in the model were characterised by [40] based on the World Wildlife Fund terrestrial ecoregions [41].Additional to metrics of permafrost extent, a brief summary of the ability of CN spun to recreate the present day carbon and nitrogen cycle is also shown in Table 3.The dynamic vegetation model component of JULES-IMOGEN calculates the proportion of the different types of land cover in the permafrost region.In general, CN spun represents the correct proportions of vegetation cover and bare soil with only slightly lower tree cover than expected.However, the model represents too much grass cover and insufficient shrub cover compared with the observed data.The modelled carbon and nitrogen stocks and fluxes broadly match the selected observations.The plant productivity is slightly higher in the model than in the observations, causing an associated higher vegetation carbon content.However, given the uncertainties in the observations, we assume the model and observations are broadly comparable.Both the model and observation-based assessments have a slightly positive net ecosystem exchange, meaning that the northern high latitudes are expected to be a sink of CO 2 from the atmosphere and in the present day.

Metric
A major recent advance for the JULES land model is the addition of a nitrogen cycle ("JULES-CN"; [14]).In this CN model configuration, the N-cycle affects vegetation by acting to limit the simulated NPP.This limitation 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 when nitrogen is assumed to be freely available (NPP pot ) compared with the actual amount achieved in the natural state when the nitrogen may be limiting (NPP).A comparison with the observations [50,51] for the tundra and boreal biomes (both regions with permafrost) suggests CN spun has a response ratio which falls at or below the lower end of the observations.This suggests that CN spun has a relatively small nitrogen limitation compared with the plausible range.
Residence times quantify how long carbon remains in the land system and can provide an extra check of whether the combined simulated carbon stocks and fluxes are correct.Whilst residence times do provide a good constraint on land models, their observation remains highly uncertain.For soils in particular, residence times are derived using the observed soil carbon stocks, which are themselves highly uncertain.Here, ecosystem residence time was estimated from a multi-annual mean on a grid cell by grid cell basis and then averaged spatially for both observations [48] and model.It is noted that the soil carbon used in the [48] data set is lower than that estimated by [49].The ecosystem residence time for the permafrost region is 154 years in the model, but 97 years in the data set of [48].Encouragingly for this assessment of the permafrost carbon feedback, the soil carbon residence time shown in Table 3 is very similar between model and observations.Changes in the N fluxes and stocks are key to the response of the permafrost in biogeochemistry in a changing climate.Firstly, we examine the present-day behaviour.Figure 2 shows the seasonal cycle of the net mineralisation, plant demand and plant uptake under present-day conditions.The plant N demand and N uptake reaches a peak slightly earlier in the summer than the net mineralisation, suggesting that the plants can become more N-limited early in the summer if there is not enough available mineral N in the soil.
The mean modelled net mineralisation for April to October averaged over the permafrost region is 5.2 g N/m 2 /year and plant uptake is 6.2 g N/m 2 /year.This is less than the mean demand for April to October of 8.4 g N/m 2 /year, suggesting that plants are not able to acquire all the mineral N they could potentially utilise.These model-based estimates are difficult to verify in full, as observations of N fluxes in the permafrost regions are limited and highly localised.The net mineralisation for April to October is higher than the field-based observations of ~1 g N/m 2 /year by [24,53].Ref. [54] suggest that the maximum daily plant N uptake during the growing season is 8.3 mg N/m 2 /day compared with our estimate of around 20 mg N/m 2 /day.These limited results suggest that nitrogen cycling in JULES in the permafrost regions may be faster than observed.This may well be related to the low N limitation compared with the observations (Table 3).

Permafrost Soils
Figure 3 shows the spatial distribution of soil carbon in the top 3 m of soil, in simulation CN spun compared with the NCSCD observations [49].In general, the model has less soil carbon in the northern permafrost regions and more in the southern permafrost regions compared with the observations.At present, JULES is missing processes such as peat soil formation, the deposition of wind-blown soil carbon, and the warmer paleo-climate conditions that may have been more favourable for carbon accumulation, so it is not expected to obtain the initial distribution of permafrost carbon accurately.Figure 4a shows the present-day profile of soil C density in the permafrost region compared with the observations (the observations are discussed in more detail in the Supplement).Although the observations of organic C have a wide uncertainty (Figure 4; panel a, grey curves), the model falls within the spread of the observed values.At the large pan-Arctic scale, the profiles (model and data) all have generally higher values at the surface which decrease with increasing depth.It is noted that the model simulates very little soil organic carbon at depths of 3 m and greater.Moreover, shown in Figure 4 are the profiles of organic (b) and inorganic (c) N density.As with the organic C density the organic N density decreases with increasing depth.Although, Figure 4b does suggest the model organic N content decreases more slowly than the observations.The simulated organic N density is comparable with the observations-slightly lower near the surface and higher at depth.The simulated present-day profile of inorganic N is shown in yellow in Figure 4c and the plant available inorganic N is shown in black.The profile of total inorganic N has a different shape to the organic C or N densities.There are higher values at the surface, which reduce just below the surface then increase to a maximum at between 1 and 2 m depth (typically just below the top of the permafrost table).These then decrease again deeper in the profile.At the surface, the higher concentrations are caused by the nitrogen deposition and fixation.A net mineralisation flux which decreases with increasing depth adds to the inorganic N deeper in the soil.This inorganic N is reduced by a temperature dependent turnover.The inorganic nitrogen below the top of the permafrost table is relatively inert at these cold temperatures and not part of the active nitrogen cycling.In addition, the plants have very limited access to this deeper inorganic N because of a combination of soil freezing and lack of roots.Figure 4c shows the plant-accessible inorganic N with the black dotted line.This is a very small fraction of the total inorganic N and indicates that the model can simulate a supply of inorganic N that is locked within the permafrost under the current climate.As the temperature warms the available inorganic nitrogen will increase when some of the additional, but unavailable inorganic nitrogen becomes active in the nitrogen cycling.
Measurements of the inorganic N profile have been made at selected sites in the Arctic.For example, Ref. [19] measured ammonium content from cores in eastern Siberia and found that ammonium was significantly higher in the permafrost than in the active layer-between 2 and 40 times depending on location.Additionally, Ref. [53] found, for a site in Alaska, that extractable dissolved inorganic nitrogen was 35 times higher in the permafrost than in the active layer.This is larger than the seven-fold increase found by [24].Figure 4c suggests that the permanently frozen soils have about four times as much inorganic nitrogen than the active layer.This falls towards the lower limit of the range of available site observations.

Soil Nitrogen
In a changing climate, both the net mineralisation flux and the plant N uptake increase with increasing temperature.However, there is asymmetry in time between the increase in net N mineralisation which occurs later in the year and increase in plant N uptake which occurs earlier in the year leading to a build up of the inorganic N pool in the autumn and a depletion in the spring.This also results in more inorganic N nearer the surface and less deeper in the profile (Figure 5a).In the future, any thawing of permafrost or change in vegetation distribution will change the proportion of the inorganic N available for plant uptake.This response is shown in Figure 5 for the permafrost region at a range of global mean temperature changes.As the temperature increases, the soil thaws and more of the inorganic N becomes available for the plants to uptake.Figure 5a shows a reduction in the total inorganic N pool at depths greater than 30 cm with increased temperature alongside an increase in the inorganic N pool at the surface.It should be noted that the availability of mineral N to the plants is highly dependent on the root distribution relative to the frozen soil.In general, below-ground plant traits are modelled less comprehensively than their above-ground characteristics, yet in high latitudes, these features will impact simulated ecosystem response to permafrost thaw [22,55].

Nitrogen Limitation
The simulated change in N limitation over the permafrost region and over time is shown via the response ratio in Figure 6a.A higher response ratio reflects more N limitation and a lower response ratio means less N limitation with a response ratio of 1 meaning no N limitation.Figure 6a shows that during the 20th century when there is a relatively small amount of climate change, N limitation is relatively stable.Over this period, despite the increase in atmospheric CO 2 , the vegetation is unable to take up as much CO 2 as it could potentially use because it does not have enough N available in the soil.However, as the permafrost begins to thaw and there is an increase in plant available N from the soil which fertilises the vegetation and increases the amount of carbon it can take up.This reduces the response ratio.Figure 6a shows that reduction in nitrogen limitation in the Arctic is greatest in the high emissions scenario, RCP8.5.The spatial distribution of the response ratio is shown in Figure 6b,c.The model suggests that the vegetation is more N limited in the tundra regions than in the boreal regions in the present day (see also Table 3).By 2200 for the RCP8.5 scenario, only the very high Arctic tundra show notable amounts of N limitation.

Implications for the Carbon Cycle
Figure 7 shows the transient response to the applied climate change of the permafrost extent along with the vegetation and soil carbon and net ecosystem productivity for the different scenarios.There is a gradual loss of permafrost from ~15 million km 2 in the present day to between 12 and 14 million km 2 depending on scenario by 2100.The vegetation carbon increases in the northern high latitudes because conditions become more favourable for plant growth.This increase is 10-20 Gt C by 2100 and 10-120 Gt C by 2300.This large increase in vegetation carbon in the model is related to a poleward expansion of trees.
Table 3 shows the model has 2 million km 2 of trees at the beginning of the 21 century.This increases to 3.2 million km 2 by 2100 under RCP8.5 up to as much as 10 million km 2 by 2300 with an associated reduction in grass and bare soil.The soil carbon increases slightly between the start of the simulation and 2100 and then levels off and decreases slightly (or decreases more significantly in the case of RCP8.5).During the historical period, the net effect of changes in soil and vegetation carbon on the net ecosystem productivity is an increase in carbon uptake, i.e., the permafrost region is a net sink as in [56].Future projections using our model suggest that this net uptake reduces after ~2050 and, depending on scenario and model initialisation, the permafrost region may become a source of carbon.The permafrost region remains a sink for both RCP2.6 and RCP4.5 scenarios.For the RCP8.5 scenario CN spun shows the permafrost region becoming a small source of carbon of around 0.2 Gt C per year after 2100.

Mineral Nitrogen Vegetation Fertilisation (minNfert) and Feedback (minNfert-Feed)
Approximately half of the increase in vegetation C seen in the CN spun simulation is caused by the reduction in N limitation as the permafrost thaws (Figure 8a). Figure 8 shows this response of the vegetation in CN spun to the mineral N vegetation fertilisation.This extra uptake of carbon by the land surface can be translated to a negative feedback onto the global surface air temperature (Figure 8b). Figure 8b shows that the impact of increased mineral N reduces the global surface air temperature by up to 0.05 • C depending on scenario.It should be noted that this feedback is independent of the amount of permafrost soil carbon at the start of the simulation.The uncertainties shown as plumes in Figure 8 represent the difference in climate sensitivities and Arctic amplification between the driving climate models.Figure 9 shows that the mineral N vegetation fertilisation is independent of the initial permafrost carbon.This is in contrast to the permafrost carbon lost, which increases with increased initial permafrost carbon (See also the Appendix A).

Comparison with the Permafrost Carbon Feedback (pfCloss-Feed)
The permafrost carbon loss (pfCloss) and feedback (pfCloss-feed) for the different future scenarios were quantified using the method summarised above and discussed in detail in the Appendix.The three different versions of the CN model each with a different plausible initial labelled permafrost carbon content are shown.The top row of Figure 9 shows the loss of labelled permafrost carbon as a function of the initial labelled permafrost carbon content for 2100 (top left figure) and 2300 (top right figure).For comparison purposes, also included on these figures are the impacts of increasing available mineral nitrogen in the soil on the carbon emissions via an increase in vegetation carbon (minNfert).Moreover, pfCloss and minNfert have opposing signs and their combined interaction gives the net effect of the impact of thawing permafrost on the northern high latitudes biogeochemical cycling simulated within JULES.The bottom row shows the temperature feedbacks onto GSAT representing pfCloss-feed and minNfert-feed for the same time periods.2 for more details on how these are estimated.
As might be expected, pfCloss and pfCloss-feed are dependent on the initial labelled permafrost carbon content-the higher the initial labelled permafrost carbon the larger the C loss and the feedback temperature.This loss is between 2 and 8 Gt C by 2100.By 2300, this loss is between 5 and 75 Gt C with RCP8.5 releasing notably more labelled permafrost carbon than the other two scenarios.At the lower initial permafrost carbon contents the pfCloss is generally offset by minNfert.Whereas at the higher values of the initial carbon content the pfCloss is greater than minNfert.The net effect on global surface air temperature of pfCloss-feed (warming effect) and minNfert-feed (cooling) is relatively small by 2100-up to 1 hundredth of a °C.This increases slightly by 2300 up to about 0.1 °C.These two processes offset each other and the net effect of the two feedbacks (not shown) can vary in sign depending on the initial permafrost carbon content.

Discussion and Conclusions
This paper uses the JULES-IMOGEN framework to assess the effect of including a mineral nitrogen store in permafrost regions on the vegetation carbon.Such nitrogen is preserved within the present-day permafrost, but is projected to act as a fertiliser causing increased vegetation carbon uptake as the soil thaws.We find that the uptake of this mineral nitrogen by the plants reduces the N limitation of vegetation resulting in a negative feedback onto the climate by slowing the rate of increase of atmospheric CO 2 .Using our JULES-IMOGEN framework, we find this extra N fertilisation causes a negative (stabilising) climate feedback potentially equal in magnitude to the positive (destabilising) permafrost carbon climate feedback.It should be noted that the permafrost carbon feedback is strongly dependent on estimates of the soil carbon initial conditions.That is, the more permafrost carbon present at the start of the simulation, the greater the positive feedback and the more likely the net effect of including both permafrost carbon and nitrogen biogeochemistry in JULES is a positive feedback onto the global surface air temperature.
The C and N cycle used within the JULES model [14] may well have less N limitation in the Arctic compared to the observations.In addition, the model may not have quite as much frozen mineral N as is observed.Moreover, JULES does not include processes such as invasive pests, pathogens, or altered herbivory regimes which could impact the relative importance of N availability.These factors mean that JULES-IMOGEN may well under-estimate the negative mineral N feedback on vegetation carbon.There are additional nitrogen cycling processes not discussed here which will also feedback onto the climate.For example, the response of the N cycle to changes in hydrology will impact the balance of nitrogen loss between leaching and denitrification.Any subsequent change in nitrous oxide emissions will feedback onto the global climate [57].
where the climate, CO 2 emissions and nitrogen deposition are set to that at the start of the simulation, i.e., there is no change in climate drivers.Therefore, in the simulations of the control-ens any loss or drift in the labelled permafrost carbon is not related to external forcing.Figure A1 shows the loss in labelled permafrost carbon in a stable climate for the three different simulations with varying initial permafrost carbon amounts shown in Table 1.This loss is dependent on both the initial permafrost carbon and the time after the initialisation of the simulation.In general, the more labelled permafrost carbon present at the start of the simulation the greater the loss with the annual loss decreasing slowly over time.At the beginning of the simulation, the loss is around 0.008% of the labelled permafrost carbon each year and falls to around 0.007% per year after 350 years.This gives a total loss of 6 Gt C by 2000, reaching around 16.5 Gt C by 2300 for CN spun (black line; Figure A2).  Figure A2 shows the control-ens for CN spun (black line) alongside the equivalent transient-ens.Ref [6] used a set of simulations equivalent to the transient-ens but used JULES-C and did not remove the drift.They ( [6]) found a loss of labelled permafrost carbon of ∼40 Gt C by 2100 (independent of scenario) and an ensemble mean of ∼70 Gt C for RCP2.6/∼200Gt C for RCP8.5 by 2300.In this latest version using the CN model ∼13 Gt C are lost by 2100 increasing to between 20 and 60 Gt C by 2300 depending on scenario.If the drifts were comparable in JULES-C [6], the impact of adding the interactive N cycle is to reduce the loss of permafrost carbon.This is likely related to differences in the soil carbon dynamics in JULES-C compared with the interactive CN configuration, in particular the depth distribution of the carbon.The loss of labelled permafrost carbon caused by the external forcing is the difference between the control-ens and transient-ens and is around ∼5 Gt C in 2100 and ∼8 Gt C for RCP2.6/∼44Gt C for RCP8.5 by 2300.
The permafrost carbon feedback temperature is quantified by comparing equivalent simulations where only the non-permafrost carbon is visible to IMOGEN (pfCloss-feed-ens) with those where all of the soil carbon is visible to IMOGEN (transient-ens).In this paper, the drift in the labelled permafrost carbon is accounted for in the prescribed CO 2 emissions in the pfCloss-feed-ens (see Table 2).Figure A3 shows the additional change in the global mean temperature caused by the loss of labelled permafrost carbon caused by the external forcing is 0.005 to 0.012 °C by 2100 and 0.01-0.04°C by 2300.As expected, this is lower than the additional increase of 0.01-0.1 °C (5th-95th percentile) by 2100 and 0.02-0.28°C by 2300 found by [6].Here, we assess the permafrost carbon and loss for the simulations highlighted in Table 1.This enables us to evaluate the sensitivity to initial conditions.The three different CN simulations-CN spun , CN 1.5spun and CN 2spun have increasing amounts of permafrost carbon in the pre-industrial which fall within the range of the observed permafrost carbon quantified in the Supplement.In CN spun , the permafrost carbon reaches its spun up state with a total of 484 Gt C (Figure A4).This is at the low end of the permafrost carbon found from the observation constraint (Supplement Figure ).In CN 1.5spun and CN 2spun the labelled permafrost carbon was prescribed at the start of the simulation to be equivalent to 1.5 and 2 times the spun up state of CN spun .These both fall within the range suggested by the observational constraint with CN 2spun at around the upper limit (Supplement Figure ). Figure A4 shows the spatial distribution of the total soil carbon in the top 3m for the three different versions of JULES CN.Additionally shown are the NCSCD soil carbon data for the top 3 m [49].There are some differences in the spatial distribution with more soil carbon in southern Asia in the model than in the observations and less in southern Canada but the overall amount is very similar.
Figure A5 shows the sensitivity of the loss of permafrost carbon to the initial conditions after the drift in the labelled carbon has been removed.In general, by 2100 less than 10 Gt C of permafrost carbon is lost from any of the simulations.This increases after 2100 and by 2300 to between 10 and 120 Gt C depending on initial condition.Overall, as might be expected, the scenarios with the greater increase in global mean temperature have a greater loss of permafrost carbon.The width of the plumes in Figure A5 give an indication of the impact of the differences in climate sensitivities between the driving models on the permafrost carbon lost.By 2300, this is between 15 and 30 Gt C.However, the difference between the permafrost carbon lost for CN spun and CN 2spun is around 40 Gt C. Therefore the uncertainties caused by differences in the initial permafrost carbon content are larger than the uncertainties caused by differences in the climate sensitivities of the driving climate models This suggests that, in order to reduce projection uncertainty, the initial permafrost carbon needs to be suitable defined.As the initial carbon content increases the total amount of labelled permafrost carbon lost increases.However, Figure A6 shows this loss as a percentage of the initial permafrost carbon is relatively independent of the initial conditions for the CN configurations.For example, CN 2spun loses more than twice as much permafrost carbon than CN spun (Figure A5) but CN 2spun also has approximately twice as much permafrost carbon present initially compared with CN spun .Therefore, as a percentage of initial permafrost carbon, both CN spun and CN 2spun lose approximately the same amount as is shown in (Figure A6).They lose about 1% of the initial permafrost carbon by 2100 and this is relatively independent of scenario.This percentage loss rises by 2300 to 2% for RCP4.5 and 7% for RCP8.5.All three CN versions have plausible amounts of initial carbon and hence plausible losses of labelled permafrost carbon.The positive permafrost carbon feedback onto GSAT is shown in Figure A7 as a function of time.As expected, it is dependent on the initial carbon content and ranges between 0.001 and 0.03 °C by 2100 and 0.002 and 0.09 °C by 2300, with the higher the initial carbon content the higher the feedback.In summary, Figures A5-A7 demonstrate that it is important to quantify the amount of permafrost carbon present at the start of the simulation in order to reduce uncertainties in the global temperature response to permafrost carbon loss (Figure A7) in a changing climate.

Figure 1 .
Figure 1.Permafrost extent (a) simulated by CN spun model configuration for the period 2000-2009 and (b) from the observations provided by [1].The model assumes either no permafrost or 100% permafrost, whereas the observations capture subgrid scale heterogeneity and therefore give the proportion of permafrost in each grid cell.Superimposed in yellow on each plot is the zero °C annual mean isotherm from the WFDEI 2000-2009 2 m air temperature [39].

Figure 2 .
Figure 2. Ensemble mean seasonal cycle of net N mineralisation, N demand and N uptake for the period 2000-2009 and spatially averaged over the simulated permafrost region shown in Figure 1a.

Figure 3 .
Figure 3. Distribution of the total soil carbon in top 3 m (kg m −2 ) for CN spun (a).Shown in (b) is the total soil carbon for the top 3 m of the NCSCD observations [49].

Figure 4 .
Figure 4. Profile of (a) organic C, (b) organic N and (c) inorganic (mineral) N within the permafrostfor the period 2000-2009 simulated using the CN spun configuration and the transient-ens.The observations are derived from a range of available data sets described in more detail in the Supplement.

Figure 5 .
Figure 5. Profile of (a) total inorganic N and (b) plant available inorganic N along with its change over the future for the permafrost region using the CN spun configuration and the transient-ens for a range of GSAT changes in °C.

Figure 6 .
Figure 6.(a) The time series of the annual and ensemble mean response ratio (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 the natural state calculated from the area-total fluxes) for the CN spun configuration and transient-ens simulation over the permafrost region.The spatial distribution of the response ratio for (b) the year 2000 and (c) the year 2200.A higher response ratio reflects more N limitation and a lower response ratio means less N limitation.A response ratio of 1 means no N limitation.

Figure 7 .
Figure 7. Output from the transient-ens simulations for the three RCP scenarios and CN spun .The changes in permafrost extent, vegetation carbon, soil carbon and net ecosystem productivity (NEP) are shown.

Figure 8 .
Figure 8.(a) The time series of the extra vegetation C uptake after fertilisation by mineral N released from thawing permafrost (minNfert).(b) The time series of the global surface air temperature (GSAT) response caused by the mineral nitrogen vegetation fertilisation feedback (minNfert-feed).See Table2for more details on how these are estimated.

Figure 9 .
Figure 9. Permafrost carbon loss (pfCloss) and permafrost carbon feedback (pfCloss-feed) for the different RCP scenarios.Additionally shown, for comparison purposes, are the impacts of increasing the available mineral nitrogen in the soil on the carbon emissions and the global temperature feedback (minNfert).

Figure A1 .
Figure A1.The annual loss of labelled permafrost carbon (in Gt C per year) caused by the on-going small exchange between the carbon above the permafrost table and permafrost carbon below the permafrost table (denoted labelled permafrost carbon in this paper) in a stable climate using the control-ens of the JULES-CN configuration.

Figure A2 .
Figure A2.The cumulative loss of labelled permafrost carbon caused by the on-going small exchange between the carbon above and the labelled carbon below the permafrost table in a stable climate for CN spun .The black line is for the control-ens and the shaded areas show the range of the transient-ens and the differences between the black line and shaded area represent the loss of permafrost carbon.

Figure A3 .
Figure A3.The permafrost carbon feedback temperature for the three different future scenarios for CN spun .

Figure A4 .
Figure A4.Spatial distribution of soil carbon in kg m −2 for the top 3 m of the soil and the three different versions of JULES CN spun , CN 1.5spun , and CN 2spun .Additionally shown are the NCSCD observed soil carbon data for the top 3 m [49].

Figure A5 .
Figure A5.Permafrost carbon lost for the different model initialisation shown in Table2.

Figure A6 .
Figure A6.Loss of permafrost carbon as a percent of the initial labelled permafrost carbon and plotted as a function of initial permafrost carbon.

Figure A7 .
Figure A7.Permafrost carbon feedback temperature onto GSAT for the different JULES configurations and scenarios.

Table 2 .
List of ensemble simulations assessed here-the transient simulation plus diagnostic experiments.These diagnostic experiments are required to quantify the mineral N fertilisation on vegetation C (minNfert) and its climate feedback (minNfert-feed) and the permafrost C loss (PfCloss) and permafrost climate feedback (PfCloss-feed).Quantifying the vegetation fertilisation by mineral N released from thawing permafrost.The vegetation is only allowed access to the mineral N that is in the soil layers that are not permafrost at the start of the simulation and cannot access mineral N in any subsequently thawed permafrost.
• minNfert is the difference in vegetation C (Gt C) between minNfert-ens and transient-ens.minNfert-feed-ensDiagnosticexperimentsQuantifying the negative climate feedback onto global surface air temperature (GSAT) from fertilisation of the vegetation by newly thawed mineral N. The prescribed CO 2 emissions are increased by the increase in vegetation C diagnosed by the minNfert-ens.•minNfert-feedisthe difference in GSAT (°C) between minNfert-feed-ens and transient-ens.control-ensStationarysimulationsQuantifying the non-climate related drift in permafrost C emissions.The climate, prescribed CO 2 emissions and N deposition are set to the values at the start of the simulation (see Section 2.1).