2.1 N₂O budget and production mechanisms

Prather et al. (2015) estimated that the preindustrial atmospheric lifetime of N₂O was 123 years. Together with the atmospheric N₂O concentration from ice cores, this figure constrains the total net preindustrial N₂O source to 10.5±1 Tg N yr⁻¹. Marine N₂O emissions were recently estimated by an observation-constrained approach, using water-column and surface N₂O observations as targets, to be 4.6 (±1 standard deviation range: 3.1 to 6.1) Tg N yr⁻¹. This calculation implies a natural terrestrial N₂O source of 5.9 (4.1 to 7.7) Tg N yr⁻¹ (Battaglia and Joos, 2018), in line with IPCC AR5 estimates of 6.6 (3.3 to 9.0) Tg N yr⁻¹.

N₂O is produced through a variety of pathways both in the ocean and on land and N₂O production is closely linked to the flows of C and N (Wrage et al., 2001; Chapuis-Lardy et al., 2006; Kato et al., 2013; Battaglia and Joos, 2018; Trimmer et al., 2016; Babbin et al., 2015; Gilly et al., 2013; Bange, 2008; Butterbach-Bahl et al., 2013; Firestone and Davidson, 1989). The dominant pathways of net N₂O production are thought to be respiratory denitrification (${\mathrm{NO}}_{\mathrm{3}}^{-}\to {\mathrm{NO}}_{\mathrm{2}}^{-}\to \mathrm{NO}\to {\mathrm{N}}_{\mathrm{2}}\mathrm{O}\to {\mathrm{N}}_{\mathrm{2}}$), under low oxygen conditions, and autotrophic nitrification (${\mathrm{NH}}_{\mathrm{3}}\to {\mathrm{NO}}_{\mathrm{2}}^{-}\to {\mathrm{NO}}_{\mathrm{3}}^{-}$), mediated by archaea, bacteria, and fungi. In addition, heterotrophic nitrification, which is the oxidation of organic N to nitrite (${\mathrm{NO}}_{\mathrm{2}}^{-}$) and subsequent reduction to N₂O by incomplete denitrification, was found to be the dominant path in a grassland ecosystem after fertilizer application (Moser et al., 2018). N₂O is produced as a byproduct of nitrification and as an intermediate product during denitrification. N₂O is also produced from nitrite through nitrification, often termed nitrifier denitrification, through anaerobic ammonium oxidation, chemautotrophic denitrification, abiotic processes, or photoautotrophic organisms in cryptogamic covers (Lenhart et al., 2015; Butterbach-Bahl et al., 2013).

Terrestrial N₂O production and emissions depend sensitively on environmental factors, including precipitation, soil temperature, soil moisture, soil texture, soil oxygen concentration or pH, topography, as well as substrate and nutrient availability and nutrient addition by deposition or fertilizer application (Zhuang et al., 2012; Stehfest and Bouwman, 2006; Wang et al., 2017). Field data- and model-based emission estimates show the highest emissions of N₂O in moist tropical areas and lower emissions at high latitudes (Stehfest and Bouwman, 2006; Zhuang et al., 2012; Werner et al., 2007; Potter et al., 1996; Xu et al., 2017; Xu-Ri et al., 2012; Wells et al., 2018). Moist soils typically show relatively high N₂O emission (Butterbach-Bahl et al., 2013; Zhuang et al., 2012), but the environmental dependencies of N₂O emission are often complex (Diem et al., 2017; Müller et al., 2015; Schmid et al., 2001a; Matson et al., 2017) and dependent on the production pathway (Kool et al., 2011). Higher N₂O fluxes at high soil water contents have been reported from laboratory and field studies and linked to increasing denitrification activity in response to reduced oxygen diffusion into the soil (Arias-Navarro et al., 2017). The sensitivity of denitrification to temperature is found to be higher than for CO₂ emissions from soil organic matter decomposition, and a positive feedback of soil warming on N₂O emissions is expected (Butterbach-Bahl et al., 2013). Warming treatment increased measured N₂O emissions in boreal peatlands (Cui et al., 2018). However, responses to warming treatment are found to be highly variable across a range of conditions and ecosystems, and it remains unclear whether warming will increase or reduce regional to global N₂O emissions (Dijkstra et al., 2012). Positive or neutral responses in N₂O emissions have been found in field experiments under elevated CO₂ in temperate and boreal forests and grasslands (van Groenigen et al., 2011; Dijkstra et al., 2012; Regan et al., 2011; Moser et al., 2018; Zhong et al., 2018). Nitrogen addition by mineral and organic fertilizer causes enhanced N₂O emissions. Emission factors are reported to depend sensitively on soil pH and are typically estimated to be around 0.5 % to 2 % of added N (Charles et al., 2017; Wang et al., 2017), but may vary by more than an order of magnitude.

The amount of N₂O produced per unit of N converted varies with environmental conditions and production pathway. For nitrification and denitrification this yield (or emission) factor depends on substrate availability linked to soil organic matter decomposition and C:N stoichiometry, on oxygen level influenced by soil moisture status and respiration, and on acidity and temperature (Diem et al., 2017; Davidson et al., 2000; Firestone and Davidson, 1989; Smith, 1997; Phillips et al., 2015; Saggar et al., 2013). In addition, different N₂O production pathways and N loss processes and their relative importance may evolve through time and influence the N₂O yield on local to global scales.

2.2 C–N–N₂O coupling

The main flow paths of reactive N and its link to N₂O production and C flows on land (Gruber and Galloway, 2008; Butterbach-Bahl et al., 2013; Vitousek et al., 2013; Zaehle, 2013; Firestone and Davidson, 1989) are schematically sketched in Fig. 1. These ecosystem flows can be assigned to an internal and an external N cycle. Within an ecosystem (green arrows in Fig. 1), N is primarily taken up in the form of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ by plant roots to support growth, while a large fraction of reactive N is taken up by soil microbes and fungi and is thus immobilized. Organic N is converted back to inorganic N during the mineralization of litter and soil organic matter. Gross N mineralization is thereby modified by the decomposers' carbon-use efficiency (Manzoni et al., 2008). Net N mineralization (as in Fig. 1) thus reflects the modified gross N mineralization minus N immobilization. The ${\mathrm{NO}}_{\mathrm{3}}^{-}$ pool is replenished by nitrification, the conversion of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ to ${\mathrm{NO}}_{\mathrm{3}}^{-}$. We note the acid–base equilibrium between ${\mathrm{NH}}_{\mathrm{4}}^{+}$ and NH₃ in soil water; for simplicity, we generally refer to ${\mathrm{NH}}_{\mathrm{4}}^{+}$ only.

Figure 1Terrestrial C and N cycles. Simple schematic of global annual C fluxes (black), recycling (green) and input/output (blue) N fluxes, and N₂O emissions (magenta) from nitrification and denitrification processes. Fluxes are illustrated with a quantitative budget from the LPX-Bern(v1.4N) model for preindustrial conditions (1450 AD). Plant net primary productivity (NPP) sets the N demand that is satisfied by N uptake from ammonium (${\mathrm{NH}}_{\mathrm{4}}^{+}$) and nitrate (${\mathrm{NO}}_{\mathrm{3}}^{-}$). N enters the ecosystem mainly through biological N fixation (BNF) by plant symbiotic and free-living microorganisms in the soil.

Download

Turning to the external cycle, reactive N enters land ecosystems (blue arrows in Fig. 1) through the conversion of dinitrogen (N₂) to organic N and eventually to ${\mathrm{NH}}_{\mathrm{4}}^{+}$ by BNF (Cleveland et al., 1999; Vitousek et al., 2013; Zaehle, 2013; Sullivan et al., 2014; Xu-Ri and Prentice, 2017), through rock weathering (Houlton et al., 2018), and the deposition (Lamarque et al., 2011; Vet et al., 2014; Dentener et al., 2006) of NH_x and NO_y (including sources by lightning). Reactive N is lost from land ecosystems through gaseous losses (including N₂O), leaching of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and dissolved organic matter by runoff, and emissions of N compounds by fire (Bouwman et al., 2013; Ciais et al., 2013; Hedin et al., 1995, 2003). On a larger scale, N lost by fires will be deposited again, and a large part of this N flux is therefore fed again to land ecosystems. In this sense, the fraction of the fire flux not lost to the ocean may be viewed as belonging to the internal global land N cycle. N leached as dissolved organic N is typically remineralized downstream and may undergo nitrification and denitrification or be taken up by aquatic organisms. N may also be absorbed by minerals and become unavailable for plants and microbiological assemblages. Following mass balance, losses of reactive N from ecosystems and changes in N stocks have to be compensated by N inputs, mainly by BNF.

The external and internal N cycles are coupled. Reactive N in mineral forms serves as a substrate for the N loss fluxes (external cycle) as well as a reservoir for plant N uptake (internal cycle). Following mass balance, mineral N concentrations change until the balance of N input by BNF and other sources matches N loss and net ecosystem N uptake (reactive N uptake minus net N mineralization). For example, in a growing ecosystem an increasing amount of N is taken up by plants and converted from inorganic to organic forms. This net ecosystem N uptake tends to deplete reactive N in mineral forms. Correspondingly, N loss fluxes (including N₂O) would decrease and N limitation of plant growth would increase if N sources, such as BNF, do not adjust to the increasing ecosystem N demand. Whether N limitation increases or not in a growing ecosystem therefore depends critically on the flexibility of N input, hence BNF.

Empirical evidence from N-addition experiments, synthesized by Lu et al. (2011) and Niu et al. (2016), shows that the uptake of N by plants, net primary productivity, and biomass as well as NH_x and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ pools in soils, nitrification, nitrate leaching, denitrification, and N₂O emissions all increase simultaneously with N input. This finding points to a tight coupling between the availability of reactive N for plant growth, nitrification, and denitrification fluxes and N₂O. N₂O production on land is predominantly associated with denitrification and to a smaller extent with nitrification. Large gaseous losses and N₂O emissions are thus indicative of a large N throughput (input and loss), where ecosystem functioning may have adjusted by exploiting energetically costly but evolutionarily advantageous BNF to compensate for the losses (Batterman et al., 2013; Pons et al., 2007; Vitousek and Hobbie, 2000) or where losses are compensated by N input from weathering or deposition.

In summary of Sect. 2.1 and 2.2, changes in terrestrial N₂O emissions may be linked (i) to changes in the magnitude of reactive N entering and leaving ecosystems and (ii) to changes in the N₂O yield per unit reactive N converted in land ecosystems.

3.1 The LPX-Bern(v1.4N) dynamic global vegetation model

The dynamic global vegetation and land surface process model LPX-Bern (“Land surface Processes and eXchanges” model as implemented at the University of Bern, version 1.4N) (Lienert and Joos, 2018b) is applied here in transient mode over the last 21 000 years (Ruosch et al., 2016). The LPX-Bern model describes dynamical vegetation and terrestrial biogeochemical processes, integrates representations of non-peatland (Gerten et al., 2004; Joos et al., 2004; Sitch et al., 2003) and peatland (Spahni et al., 2013; Wania et al., 2009) ecosystems and their C and N dynamics (Stocker et al., 2013; Xu-Ri and Prentice, 2008; Xu-Ri et al., 2012), and describes the dynamic evolution of wetland and peatland extent (Stocker et al., 2014). The model calculates the release and uptake of the trace gases CO₂, N₂O (Stocker et al., 2013; Xu-Ri and Prentice, 2008; Xu-Ri et al., 2012), and CH₄ (Spahni et al., 2011; Wania et al., 2010; Zürcher et al., 2013).

Vegetation is represented by plant functional types (PFTs) that are in competition for resources (water, light, N) on each grid cell. Here a version with 15 PFTs is used. Eight generic tree PFTs and PFTs for C3- and C4-type grasses grow on natural land (excluding peat) and former peat. Two PFTs representing peat mosses and flood-tolerant C3 graminoids as well as three flood-tolerant tropical PFTs grow on peat and wetlands. The model accounts for the dynamic coupling of C and water cycles through photosynthesis and evapotranspiration, which also defines plant water use efficiency (Saurer et al., 2014; Keller et al., 2017). Seven C and N pools per PFT represent leaves, sapwood, heartwood, fine roots, aboveground leaf litter, aboveground woody litter, and belowground litter. Separate soil organic C and N pools receive input from litter of all PFTs. LPX-Bern uses a vertically resolved soil hydrology, heat diffusion, and an interactive thawing–freezing scheme (Wania et al., 2009).

The LPX-Bern vegetation and soil components interact with a dynamic N-cycle module (Xu-Ri and Prentice, 2008; Xu-Ri et al., 2012), accounting for the uptake of mineral N by N immobilization in soils (Bengtsson et al., 2003; Li et al., 2017; Gütlein et al., 2017) as in Xu-Ri and Prentice (2017). The module describes the relevant N and N₂O fluxes and pools for plants and soils as schematically depicted in Fig. 1 and summarized below.

In LPX-Bern, the source of reactive N is implied by maintaining prescribed soil N:C ratios associated with each of the plant functional types, reflecting their different litter chemistries and decomposer assemblages. Due to lower N:C ratios of litter than soil pools, the transfer of mass from litter to soil pools during litter decomposition therefore implies a given amount of N required to satisfy the soil N:C ratio. The required N is partly satisfied by a flux representing immobilization of mineral N. The remainder determines the total input of reactive N into the ecosystem, implicitly subsuming symbiotic and asymbiotic BNF, and any other potential N sources that may support plant growth, in addition to prescribed N deposition. The amount of N input required to close the N balance of soils and to maintain the soil pools at their high N:C ratios depends on the flux representing N immobilization. Constant fractions (frac_soil_immob, frac_litter_immob) of the N flux released by soil or litter remineralization are immediately returned to its pool of origin. Hence, the choice of these parameters simultaneously co-determines reactive N input rates and net N mineralization.

We assumed the parameters frac_soil_immob and frac_litter_immob to be invariant over time and space and calibrated their values (Table A1) to match a total preindustrial reactive N source of 128 Tg N yr⁻¹. This value implicitly includes contributions from symbiotic and asymbiotic BNF as well as other inputs of reactive N not included in the prescribed N deposition. Estimates of the global BNF for non-agricultural ecosystems are uncertain. They range from 40 to 470 Tg N yr⁻¹. Cleveland et al. (1999) used 100 plot-scale estimates of BNF to estimate global BNF on natural ecosystems to 195 Tg N yr⁻¹ (range: 100 to 290 Tg N yr⁻¹). Vitousek et al. (2013) suggest a plausible range for preindustrial BNF of 40 to 100 Tg N yr⁻¹ by computing BNF as the difference from all other global sources and sink fluxes of N. Cleveland et al. (2013) estimate symbiotic BNF to 105 Tg N yr⁻¹, based on cost–benefit modelling for N fixation. The same authors estimate asymbiotic N fixation to 22 Tg N yr⁻¹ by upscaling measurements reported in Cleveland et al. (1999). In contrast, Elbert et al. (2012) estimate asymbiotic N fixation by cryptogamic covers alone to 49 Tg N yr⁻¹ (27–99 Tg N yr⁻¹) by integrating experimental data from 200 studies. In addition, rock weathering is estimated to add 15 (10 to 20) Tg N yr⁻¹ to land ecosystems (Houlton et al., 2018). Xu-Ri and Prentice (2017) estimate global N sources to 340 (230–470) Tg N yr⁻¹ for the parameter settings they adopted within the LPJ-DyN model; this estimate includes contributions from rock weathering and other inputs that are not explicitly prescribed or simulated by LPJ-DyN. Meyerholt et al. (2016) implemented six different BNF formulations in their model to predict modern BNF ranging from 108 to 148 Tg N yr⁻¹.

Plant net primary productivity (NPP) and a prescribed constant N:C ratio of new production in different tissues sets the N demand that is satisfied by N uptake from ${\mathrm{NH}}_{\mathrm{4}}^{+}$ and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ pools which in turn depend on net N mineralization fluxes from litter and soil pools and loss fluxes of reactive N (e.g. denitrification, leaching, volatilization) (Fig. 1). In case that available inorganic N (sum of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ and ${\mathrm{NO}}_{\mathrm{3}}^{-}$) is insufficient to meet the demand, NPP is down-regulated, thereby inducing an effect of N limitation. The implied N source tends to re-establish a balance between the input and the loss of reactive N.

Nitrification is assumed to be proportional to the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ soil pool with a temperature-dependent rate coefficient. Denitrification is modelled as a two-step process whereby ${\mathrm{NO}}_{\mathrm{3}}^{-}$ is converted to nitrite (${\mathrm{NO}}_{\mathrm{2}}^{-}$) and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ is further converted to N₂, following Michaelis–Menten kinetics with dependence on the substrates ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ as well as on labile C availability. Reaction rates are temperature dependent. ${\mathrm{NO}}_{\mathrm{3}}^{-}$ leaching depends on soil ${\mathrm{NO}}_{\mathrm{3}}^{-}$, available water holding capacity, and daily runoff. Reactive N is also lost from a grid cell by fire, assuming complete loss of N from burned vegetation, and NH₃ volatilization.

N₂O emissions are computed by assuming that fractions of the N fluxes associated with denitrification, nitrification, and N leaching are released as N₂O. The calibration of the denitrification and nitrification yield factors is done to satisfy the constraints from the global atmospheric N₂O budget, which suggests global terrestrial emissions to be 5.9 Tg N yr⁻¹ at preindustrial (Battaglia and Joos, 2018). Since global terrestrial N₂O emissions depend also on the N loss rates (denitrification and nitrification), we first calibrated the parameters frac_soil_immob and frac_litter_immob, which co-determine these rates, and then calibrated the yield parameters, given the N loss rates. For denitrification, the globally dominant N₂O source path, this fraction is the product of a constant (RN2ODN) and a temperature-dependent factor f(T). f is unity at 22 ^∘C and its value roughly doubles for a temperature increase of 10 ^∘C. The amount of N₂O released per unit N denitrified is therefore higher in warm regions than in cold regions. The constant RN2ODN is set to 5.418 % of the denitrification rate. The corresponding yield fraction for N₂O emissions from nitrification is nominally set to 0.231 % at 22 ^∘C and the same temperature dependence (Smith, 1997) as for denitrification is assumed. The N₂O yield for the leaching N flux is assumed constant and, for simplicity, set to the nominal value of nitrification (0.231 %).

The two-step calibration described above resulted in yield factors that are higher than the range of published estimates. The global mean yield for denitrification, expressed as N₂O per N₂ produced, is 5.6 % and thus higher than the range of estimates (0.2 %–4.7 %) summarized by Xu-Ri and Prentice (2008). Similarly, the global yield for nitrification (0.26 %) is higher than observation-based estimates (0.01 %–0.2 %). The mismatches in these estimates for yield may suggest that current best estimates for the N source, N₂O yield, and preindustrial N₂O emissions are not fully consistent. We carried out a sensitivity simulation to explore uncertainties: the immobilization fraction is set to 0 % for litter and to 26.39 % for soil mineralization, leading to a preindustrial N source of 523 Tg N yr⁻¹. This is about a factor of 4 higher than in the standard setup. Correspondingly, yield factors are about a factor of 4 lower in this sensitivity run than in the standard simulation.

The sensitivity of simulated N₂O emission changes over the deglacial period to these parameter choices is relatively small, while the absolute magnitude of the N source has some implications for N stress and thus NPP. The increase in NPP over the deglaciation is larger in simulations with a high compared to a low N source (10.1 versus 5.9 GtC yr⁻¹ in the standard). Importantly, the difference in relative changes in modelled global N sources is small (16 % versus 10 % in the standard), and the deglacial increase in N₂O emissions is only 0.2 Tg N yr⁻¹ higher in the sensitivity run than in the standard run (see Sects. 3.2 and 5), despite the large difference in the implied N source. Thus, related model-based conclusions for N₂O emissions are not sensitive to the parameter settings for the yield and the flux representing immobilization.

3.2 Setup for transient glacial–interglacial simulations

A previously described LPX-Bern model setup for glacial–interglacial simulations is applied (Spahni et al., 2013; Ruosch et al., 2016) and input data are shown in Fig. 2. The evolution of monthly temperature, precipitation, cloud cover, and number of wet days, annual atmospheric CO₂ (Joos and Spahni, 2008), orbital insolation changes (Berger, 1978) modulating plant available light, and topography changes through ice-sheet and sea-level changes imposed by ICE-5G (Peltier, 2004) are prescribed. The monthly climate data are obtained by combining monthly values from the observation-based, modern climatology compiled by the Climate Research Unit (CRU) (Mitchell and Jones, 2005) with monthly anomalies for the past 21 000 years from a transient climate simulation over this period (TraCE-21ka) (Liu et al., 2009; Otto-Bliesner et al., 2014) with the Community Climate System Model, version 3 (CCSM3) maintained by the National Centre for Atmospheric Research (NCAR). Here the LPX-Bern model was run with a spatial resolution of 3.75^∘ longitude × 2.5^∘ latitude and a daily time step was applied in the photosynthesis, water, and N-cycle modules. Simulations started from an equilibrated spin-up at 21 ka. Annual N deposition from the atmosphere, distributed across days according to precipitation within a year, is prescribed at preindustrial values (Lamarque et al., 2011).

Figure 2Aggregated LPX-Bern(v1.4N) input data. (a) Evolution of atmospheric CO₂, (b) of temperature and precipitation, and (c) of changes in tropical and global land as available for plant growth relative to 21 ka (bottom). Temperature and precipitation represent mean values over global land areas that were not flooded and not ice-covered at 21 ka PB and at 0.5 ka.

Download

The same model parameter values as determined by Lienert and Joos (2018b) are used, except for the immobilization fractions and N₂O yield factors (Table A1). Regarding the N module, the parameters are the same as in previous studies (Stocker et al., 2013; Schilt et al., 2014) addressing N₂O emissions over the past and under future global warming. An exception is an adjustment in the upper limit of the fraction of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ nitrified per day from 0.1 to 0.09096 d⁻¹ at 20 ^∘C and in the yield factors for denitrification, nitrification, and leaching in response to new observation-constrained estimates of marine N₂O emissions and to a slightly revised estimate of the atmospheric N₂O lifetime from 120 to 123 years (Prather et al., 2015) as well as the simple representation of N uptake by N immobilization as discussed above.

Changes in N₂O emissions and other model outcomes are attributed to individual driving factors (temperature, precipitation, CO₂, orbital insolation, and land mask). One driver is kept at its preindustrial state in factorial simulations. BNF was kept at LGM values in an additional factorial run for each land use class and grid cell. In a further simulation, the yield factors for denitrification and nitrification were set to be constant in space and time for each land use class. The difference in results between the standard model setup (baseline) and a factorial run is attributed to the relevant driver. An interaction or synergy term, called “other drivers”, is quantified by the difference between the change in N₂O emissions (ΔeN₂O) simulated in the baseline run and the sum of the emission changes attributed to individual drivers: ${\mathrm{\Delta }{\mathrm{eN}}_{\mathrm{2}}\mathrm{O}}_{\mathrm{other}-\mathrm{drivers}}={\mathrm{\Delta }{\mathrm{eN}}_{\mathrm{2}}\mathrm{O}}_{\mathrm{baseline}}-{\mathrm{\Delta }{\mathrm{eN}}_{\mathrm{2}}\mathrm{O}}_{\mathrm{temperature}}-{\mathrm{\Delta }{\mathrm{eN}}_{\mathrm{2}}\mathrm{O}}_{\mathrm{precipitation}}-{\mathrm{\Delta }{\mathrm{eN}}_{\mathrm{2}}\mathrm{O}}_{{\mathrm{CO}}_{\mathrm{2}}}$. The dominant driver is identified as having the largest contribution to ΔeN₂O in the baseline run with the same sign of change. Grid cells that submerge under seawater or emerge from waning ice sheets during the period considered and grid cells with insignificant changes ($\left|\mathrm{\Delta }{\mathrm{eN}}_{\mathrm{2}}\mathrm{O}\right|\mathit{<}\mathrm{1}$ mg N m⁻² yr⁻¹) are excluded from the spatially resolved attribution.

The response timescales of LPX are investigated in a further “step-change” sensitivity simulation. Starting from the equilibrated spin-up at 21 ka, climatic conditions and atmospheric CO₂ are abruptly changed to conditions for 2.5 ka. The run is continued for another 1500 years with climate and CO₂ forcing for the period from 2.5 to 1 ka. The land mask is kept constant at the maximum extent possible for both LGM and late Holocene conditions. A corresponding reference simulation without step change was performed. The simulations, together with the above factorial runs, permit us to address how quickly N₂O emissions in the LPX model are able to respond to a sudden warming event, similar to the onset of Bølling–Allerød (B/A) and the end of the Younger Dryas (YD).

5.1 Simulated changes in global terrestrial N₂O emissions and spatial patterns of change

We next explore governing mechanisms of the deglacial terrestrial N₂O emissions and potential implications for the C–N cycles in the spatially resolved LPX-Bern model. LPX-Bern v1.4N simulates a general increase in global land N₂O emissions over the deglacial period (Fig. 3). The modelled deglacial increase in N₂O emissions is only 0.96 Tg N yr⁻¹ and about 47 % smaller than the reconstructed increase of 1.8±0.3 Tg N yr⁻¹ (±1 standard deviation; Fig. 3). Modelled changes are typically less abrupt than reconstructed variations. The model shows emission variations during the B/A and YD periods with peaks in emissions at the onset of the B/A (14.6 ka) and the preboreal (11.7 ka) and an emission peak around 13.5 ka and corresponding minima at 14 ka and during the YD (12.8 to 11.7 ka). Reconstruction and models show moderate and small changes in global terrestrial N₂O emissions over the last 11 000 years, the Holocene period. Simulated terrestrial N₂O emissions decreased by 0.06 Tg N yr⁻¹ between 10 and 0.5 ka, whereas reconstructed emissions increased by 0.41 Tg N yr⁻¹ over the Holocene, leading to a discrepancy between simulated and reconstructed anomalies. The reconstruction suggests relatively constant emissions from the land biosphere during HS1 and a rapid rise in emissions at the start of the B/A, whereas the model simulates steadily increasing emissions over the second half of the HS1 interval, reflecting the gradual climate warming in the TraCE-21ka climate input data during HS1 (see discussion below). The model also shows a much slower and a smaller emission increase at the YD–PB transition than reconstructed.

Figure 4Preindustrial patterns (a) and changes in terrestrial N₂O emissions for the last glacial termination (b; 11 ka minus 21 ka) and the Holocene (c; 0.5 ka minus 11 ka) as simulated with LPX-Bern(v1.4N). Numbers indicate integrated values for three latitudinal bands (90–20^∘ S; 20^∘ S–20^∘ N; 20–90^∘ N).

The changes in global terrestrial N₂O emissions are the result of spatially differentiated responses. LPX-Bern simulates high natural emissions of N₂O in the tropics and low emissions at high latitudes (Fig. 4a). The spatial pattern and the magnitude of emissions are consistent with data-based estimates of natural N₂O emissions from soils (Zhuang et al., 2012; Stehfest and Bouwman, 2006; Potter et al., 1996). At 1450 AD, emissions in the tropics can be as high as 250 mg N m⁻² yr⁻¹ and the integrated flux between 20^∘ S and 20^∘ N amounts to 54 % of the global emissions, while emissions per unit area are low at high latitudes and northern and southern extra-tropics contribute only a share of 29 % and 17 % to the global terrestrial emissions. In contrast, emissions increased strongly in the northern extra-tropics over the glacial termination (Fig. 4b), while the integrated change in emissions is negligible in the tropical belt. Large increases per unit area are simulated over the termination at mid and low latitudes on the North and South American continents, in the southern boreal zone in Eurasia, and in parts of eastern Asia, India, Indonesia, and Africa. N₂O emissions decreased in a few regions, namely in Africa around 15^∘ S and in northern Australia (Fig. 4b). Loss of tropical land due to rising sea level and the addition of land emerging from waning ice sheets (Peltier, 2004) are important drivers of modelled terrestrial N₂O emissions between 15 and 8 ka. The net loss of terrestrial N₂O emissions by land area changes was dominated by the submergence of the high-productivity Sunda and Sahul Shelf regions. This loss offsets about one-third of the G-IG increase in global terrestrial N₂O emissions. Changes in the land extent caused by changes in sea level represent an important factor for past global terrestrial N₂O emissions.

The patterns of change in terrestrial N₂O emissions (Fig. 4c), as well as spatial and seasonal patterns in precipitation and temperature, are complex for the Holocene period. Despite small changes in global terrestrial N₂O emissions during the Holocene, LPX-Bern simulated large regional shifts in source strength, linked to changes in temperature and precipitation. This includes for example a decrease in N₂O emissions from 11 to 0.5 ka in boreal Asia, in the sub-Sahara region in Africa, and in parts of the conterminous United States of America (USA) and an increase in emissions in tropical Africa, parts of Australia and Latin America, as well as in western Europe.

5.2 Biological nitrogen fixation and carbon–nitrogen coupling

In this section, we address C–N coupling in LPX-Bern and analyse the spatial patterns for the source of reactive N, soil mineral N, net primary productivity (NPP), and C stocks (Fig. 5) and their changes over the deglaciation (Fig. 6). We quantify two decisive factors for N₂O emission change in the model: (i) changes in the source of reactive N, fueling nitrification and denitrification, and (ii) changes in the N₂O yield per unit N converted.

Figure 5Preindustrial patterns of the input flux of new nitrogen (N) into ecosystems (a), the bioavailable concentration of N in the form of ammonium and nitrate in soils (b), net primary productivity (c), and total carbon in the land biosphere (d) as simulated by LPX-Bern(v1.4N). The N source in (a) does not include the prescribed N deposition. Deserts with high mineral nitrogen concentrations are masked for clarity in (b). Numbers indicate integrated values for three latitudinal bands (90–20^∘ S; 20^∘ S–20^∘ N; 20–90^∘ N).

Figure 6Biological nitrogen fixation allowed to adjust to demand (left) versus time-invariant N source (right). The figure shows simulated deglacial changes (0.5 ka minus 21 ka) for the standard setup of LPX-Bern(v1.4N) (left column) and a setup where the nitrogen source flux is kept constant at values simulated for the Last Glacial Maximum (right column) for the same variables shown in Fig. 5. Note that a decreasing N source is shown in panel (b) for land areas lost due to deglacial sea-level rise and in some tropical cells where mineral soils with a relatively high N source are replaced by peat with a relatively low N source.

We first address C and N fluxes for the preindustrial mean state (Fig. 5), before turning to deglacial change. The N source (Fig. 5a) is typically smaller than 0.5 g N m⁻² yr⁻¹ at northern mid and high latitudes and around 1 g N m⁻² yr⁻¹ in tropical rainforest and in central America, comparable to observational estimates (Sullivan et al., 2014; Wurzburger and Hedin, 2016a, b). An N source in the range of 2 to 6 g N m⁻² yr⁻¹ is simulated in many semi-arid regions, including southern Africa, the sub-Sahara region, India, northern Australia, and the southern parts of North America. Soil mineral N (Fig. 5b) is typically below 0.5 g N m⁻² in tropical forests and around 1 g N m⁻² in temperate and boreal forests and tundra regions, while higher values of up to and more than 10 g N m⁻² are simulated in extra-tropical semi-arid and arid regions. Annually integrated NPP (Fig. 5c) is largest in the tropics, while large carbon stocks (Fig. 5d) are simulated in tropical forests and in the northern boreal zone. We note that the patterns of these C and N fluxes and stocks are all different, pointing to spatially distinct relationships between these four variables.

The N₂O yield factors, i.e. the N₂O produced per unit N converted by denitrification and nitrification, are assumed to vary with temperature and thus in space and time in LPX-Bern v1.4N. In a sensitivity run, these yield factors are set constant with all other settings as in the standard. The deglacial warming leads to a higher N₂O yield in the standard compared to this sensitivity run, and 0.44 Tg N yr⁻¹ of the deglacial increase in land N₂O emissions are attributed to this change in yield (Fig. 7, black line). In other words, changes in the yield factors further amplify the increase in N₂O emissions as driven by the increase in the flow of reactive N in LPX-Bern.

Turning to changes in the sources of reactive N, the simulated N source increased by 12 % from 113 Tg N yr⁻¹ at the LGM to 128 Tg N yr⁻¹ at 0.5 ka in the standard setup. BNF, computed to match prescribed N:C ratios in soils (see methods), is implicitly assumed to adjust on the timescales of vegetation and soil carbon turnover to an increase in N demand by raising NPP and partly alleviating N limitation of plant growth in LPX-Bern. It is implicitly assumed that limitation by other nutrients does not affect the cycling of N and C through ecosystems on multi-decadal to century timescales and that nutrient input into ecosystems by deposition and weathering (plus BNF for N) is large enough to support plant growth. These assumptions are controversial (Hungate et al., 2003; Luo et al., 2004; Körner, 2015; Wieder et al., 2015b; Vitousek et al., 2010; Fatichi et al., 2014). Our novel ice-core isotope reconstruction of terrestrial N₂O emissions allows us to critically evaluate this assumption in a quantitative way on the multi-decadal to centennial timescales as relevant for the anthropogenic perturbation. To this end, we implement an N-cycle representation leading to strong long-term nutrient limitation in the model. This is achieved by fixing the rate of the variable N source to its glacial value in each grid cell and land use class and by keeping these rates of N source constant in a sensitivity simulation over the past 21 000 years. This more strongly N-limited simulation yields not an increase, as reconstructed, but a decrease in global terrestrial N₂O emissions (Fig. 3, black line). The decrease in global land N₂O emissions is related to the loss of land in response to sea-level rise (−0.11 Tg N yr⁻¹) and a small decrease in N₂O emissions on remaining land (−0.088 Tg N yr⁻¹) in this sensitivity simulation.

The comparison of changes in C and N fluxes and stocks between the simulations with variable and constant N sources also provides further insight into the N and C coupling (Fig. 6). The setup with a constant N source leads to smaller changes in the N cycle, a small deglacial decrease (−1.2 Pg C yr⁻¹) instead of an increase (5.9 Pg C yr⁻¹) in global NPP, and a decrease in terrestrial C stocks of 47 Pg C compared to an increase of 482 Pg C simulated in the standard setup.

The N source increases in many land regions in the standard setup but remains by design almost constant on non-flooded land in the simulation, with a constant N source per land use class (Fig. 6a, b). In the standard, the N source decreased over the deglaciation by 7 Tg N yr⁻¹ in the tropics, due to land loss and a decrease in parts of Africa, while the N source increased by 16 and 4 Tg N yr⁻¹ in the northern and southern extra-tropics, respectively. Corresponding changes are simulated for nitrification and denitrification in the standard setup, whereas changes in nitrification and denitrification remain small on non-flooded land in the N-limited run.

In response to the increased N input, the availability of reactive N remained roughly constant or increased in most land areas in the standard run, despite increased storage of N in plant and soil organic matter and accelerated nitrification and denitrification (Fig. 6c). A decrease in mineral N is simulated in the standard in regions where C stocks also decreased or changed little, including parts of mid-latitude Eurasia, Africa, Australia, and the southern US. In contrast, soil reactive N concentrations decreased in many extra-tropical regions in the simulation with a constant N source (Fig. 6d). Globally integrated, soil mineral N decreased by 19 % over the deglaciation in the run with a constant N source, compared to a negligible change in the standard setup.

These differences in BNF (including abiotic input) and the N cycle between the run with a constant N source and the standard case have profound impacts, not only for the emission of N₂O, but also for NPP (Fig. 6e, f) and C sequestration (Fig. 6g, h). Global NPP increased over the deglaciation by 13 % in the standard case, compared to a decrease by 3 % in the run with a constant N source. LPX-Bern predicts that most of the NPP increase is realized in the northern extra-tropics. The NPP increase in the extra-tropics is almost twice as large in the standard as in the sensitivity run. As a result of the generally higher NPP, the deglacial change in C storage is in the standard, about 402 Pg C larger in the extra-tropics and about 127 Pg C larger in the tropics than in the run with a constant N source.

In summary, the simulation with a constant N source completely fails to reproduce the reconstructed N₂O emissions from the land biosphere. The increases in N₂O yield, as well as in soil and litter C and N turnover rates, under deglacial warming are not sufficient to overcome the effect of N limitation on N₂O emissions in this sensitivity simulation. We note, however, that changes in yield due to processes not incorporated into LPX-Bern could potentially explain the reconstructed increase in N₂O emissions. If the model is allowed to satisfy the demand of N, and thereby implicitly of other elements to support the growth of N fixers, nitrifiers and denitrifiers, and plants, simulated terrestrial N₂O emissions increase substantially.

5.3 Attribution of simulated terrestrial N₂O emissions to climatic and environmental drivers

The complex spatio-temporal changes in land N₂O emissions are attributed to climatic and environmental drivers (Figs. 7 to 10). This attribution helps us to better understand the simulated changes and data–model mismatches. Globally, deglacial warming is the most important factor in the model, contributing 1.4 Tg N yr⁻¹ to the overall emission increase, followed by changes in precipitation (0.7 Tg N yr⁻¹, Fig. 7). Changes due to increasing atmospheric CO₂ and orbitally driven changes in photosynthetically active radiation slightly offset the emission increase. Non-linear interactions among the different drivers contribute negatively to the global increase. The change in land extent due to ice-sheet melting and sea-level rise leads to a reduction in terrestrial N₂O emissions (−0.5 Tg N yr⁻¹). As described in the previous section, a substantial part of the increase attributed to temperature is due to a temperature-driven change in the yield factor of nitrification and denitrification (Fig. 7, red and black lines).

Figure 7Attribution of simulated global terrestrial N₂O emission changes to individual drivers. Results are from LPX-Bern(v1.4N) for the standard and factorial simulations, where an individual driver (temperature (red), precipitation (blue), atmospheric CO₂ (green), orbital parameters affecting photosynthetic active radiation (light blue), land mask accounting for sea-level changes) is kept at the preindustrial level. The violet curve (other non-linear effect) is the difference between simulated emissions in the standard setup and the sum of all these attributed emissions; it reflects the effect of non-linear interactions between different drivers. In addition, the black curve indicates the change in N₂O emissions attributed to variations in the N₂O yield factor for nitrification and denitrification.

Download

Turning to HS1 and the transition to the B/A, the majority of the simulated emission increase is in response to changes in physical climate, i.e. changes in temperature and precipitation (Fig. 7). In other words, the simulated rise in land N₂O emissions, which occurred earlier in the simulation than the ice-core reconstruction, is not driven by CO₂ fertilization and an associated change in N ecosystem flows, but by the prescribed climate change from TraCE-21ka. It is not clear whether the model's failure to represent the reconstructed emission changes during the HS1 period and at the onset of the B/A in the standard simulation is due to deficiency in the response of LPX-Bern to early deglacial climate change, to deficiencies in the climate input data, or a combination of the two.

Individual drivers exert regionally distinct influences, and these may vary between different periods. Here, we attribute the spatial changes to changes in temperature, precipitation, CO₂, and their non-linear interactions for the glacial termination and the Holocene using factorial simulations (Fig. 8) and by examining the temperature and precipitation changes in the TraCE-21ka input data (Fig. 9). Generally, attributed changes in emissions find their counterpart in underlying changes in individual drivers, but sometimes non-linear interactions and non-additivity of individual responses hamper the attribution to individual drivers.

Figure 8Attribution of simulated global terrestrial N₂O emission changes to individual drivers for the last glacial termination (11 ka minus 21 ka; left) and the Holocene (0.5 ka minus 11 ka; right). Changes are attributed to changes in temperature (a, b), precipitation (c, d), atmospheric CO₂ (e, f), and non-linear interactions among the drivers plus the small effect due to changes in incident solar radiation (g, h).

Figure 9Changes in Northern Hemisphere summer temperature (June, July, August, panels a and b) and annual precipitation (c, d) as simulated by CCSM3 and prescribed to LPX-Bern(v1.4N) for the glacial termination (a, c) and the Holocene (b, d).

Changes in temperature over the termination caused emissions to rise substantially in Eurasia, in the conterminous USA as well as in Argentina and southern Brazil and in eastern Australia (Fig. 8a). On the other hand, regional cooling in parts of the Amazon region and in parts of Africa and India caused emissions to fall in this period. The pattern of change in N₂O emissions attributed to changes in temperature (Fig. 8b) is very different for the Holocene compared to the termination. The summer (June, July, August) cooling also found in climate data (Wanner et al., 2008) and simulated over the Holocene period in boreal Eurasia and the western part of North America results in a decrease in N₂O emissions over large parts of Eurasia and North America. The slight summer warming in eastern Canada and Scandinavia in the model has little impact on simulated emissions. An increase in terrestrial emissions is simulated in tropical Africa and Latin America, including the Amazon region, in response to changing temperatures.

The attribution of emissions to changes in precipitation (Fig. 8c, d) reveals some well-expressed dipole patterns, partly indicative of spatial shifts in precipitation regions (Fig. 9c). Attributed emission changes reflect corresponding changes in precipitation with generally increasing emissions under increasing precipitation. Increasing CO₂ exerts a generally positive influence on simulated emissions in tropical regions, except in semi-arid regions including eastern Africa and Brazil, which dominate the globally integrated signal (Fig. 8e, f). Non-linear interactions between the three drivers can be significant regionally and either enhance or reduce simulated emission changes (Fig. 8g, h). The temperature-attributed decrease in emissions over the termination at around 15 to 20^∘ S in Africa (Fig. 8a) is not linked to a corresponding change in seasonal or annual temperatures. This particular attribution as well as the attribution of negative emission changes to increasing CO₂ in this region (Fig. 8e, f) most likely reflect non-linear interactions between precipitation, temperature, and CO₂.

Figure 10 further illustrates which driver exerts the largest influence on simulated regional emission changes in a given grid cell. Over the termination, changing temperature is the dominant influence on the simulated emission changes in mid- and high-latitude Eurasia and North America, while changes in precipitation dominate emission changes south of the Sahara and in India and large parts of Australia, where temperatures were already rather high. Over the Holocene, temperature is the dominant driver in northern Eurasia, while changes in precipitation dominate the emission response in Africa. Changes in CO₂ generally play a secondary role.

Figure 10The dominant drivers for changes in terrestrial N₂O emissions over (a) the glacial termination and (b) the Holocene are temperature (red), precipitation (blue), CO₂ (green), and other non-linear effects in LPX-Bern(v1.4N) on these driver combinations (violet). The colour shading indicates the fraction of the total emission change attributed to the dominant driver typically explaining between 40 % and 90 %. Areas with minimal changes are masked.

Figure 11Response to a step-like change in (a) environmental conditions for globally averaged terrestrial (b) N₂O emissions and (c) NPP. Atmospheric CO₂ (dotted line in a) and climate are prescribed to suddenly change from Last Glacial Maximum conditions to late Holocene conditions in a sensitivity simulation with LPX-Bern(v1.4N) (“step run”, blue lines). Annual values for the step run are shown in dark blue and smoothed (31-year average) values in light blue for the step run and in dark green for a corresponding reference run. The dashed horizontal lines represent averages from the reference run over the model period 0 to 1200 years. Surface air temperature in (a) represents the average over the model's land area. Note the change in the time-axis scaling at year 100.

Download

5.4 The timescales of response to a step change in climate and CO₂

Next, we investigate potential causes of the relatively slow N₂O increase during past abrupt NH warming events in the standard simulation compared to the ice-core reconstruction. To this end, we assess on which timescales N₂O emissions adjust in LPX-Bern after a sudden change in environmental conditions. The question is whether LPX-Bern can simulate an equally fast response in N₂O emissions as reconstructed from the ice-core data for the abrupt warming events at the onset of the B/A and the end of the YD or whether the intrinsic response timescales of the model are too long to match reconstructed emissions. We construct an extreme case for rapid warming events by switching climate (and CO₂) instantaneously, step-like from LGM to late Holocene conditions in LPX-Bern (Fig. 11a).

LPX-Bern shows a fast response, followed by a relatively small century-scale adjustment (Fig. 11). About 80 % of the final response in global N₂O emissions to the step change is realized within 40 years and about 90 % within a century, while it takes about 700 years to reach a near equilibrium (Fig. 11b). Taking the atmospheric lifetime of N₂O of more than 100 years into account, such a fast multi-decadal increase in N₂O emissions, would cause a century-scale increase in atmospheric N₂O concentrations, similar to what is seen in the ice-core record. The adjustments in global nitrification and denitrification fluxes evolve similarly to those for N₂O emissions, with the main increase in these fluxes again realized within about 40 years (not shown). The global N source peaks 20 years after the step (not shown). Then, it decreases again to reach a new steady-state value, which is about 20 % higher than before the step. In contrast, the model's global N–leaching flux decreases within about a decade (by about 20 % relative to LGM) after the step, followed by a slight increasing trend over the next 1500 years. This transient evolution in modelled N leakage is indicative of a corresponding evolution of soil ${\mathrm{NO}}_{\mathrm{3}}^{-}$ availability. Modelled global NPP (Fig. 11c) increases immediately after the step, and about three-quarters of the NPP change are realized in the first year after the step. Then it takes again about 700 years to reach the new equilibrium. The fast initial response is explained as follows. NPP responds immediately to the change in environmental conditions. The associated enhanced plant N uptake depletes soil ${\mathrm{NO}}_{\mathrm{3}}^{-}$, leading to the decrease in N leakage, while warming accelerates soil decomposition and the release of ${\mathrm{NH}}_{\mathrm{4}}^{+}$. The newly assimilated C and N are allocated to leaves, sapwood and hardwood, and roots, before being released to litter and labile soil organic matter overturning on timescales from years to decades. These annual-to-decadal vegetation and litter turnover timescales govern the initial response time of N source, nitrification, denitrification, and N₂O emissions. The century-scale response is linked to the equilibration timescales of C and N in the slowly overturning soil pools as well as to the poleward expansion of vegetation. Finally, there is a small remaining offset between the reference run and the near equilibrium of the step simulation which is linked to remaining peat in the step simulation formed under LGM conditions but is absent in the reference simulation.

In conclusion, this sensitivity experiment demonstrates that N₂O emissions in the LPX-Bern model are indeed able to adjust within decades to less than a century to abrupt warming events, if prescribed forcing is changing quickly enough. We conclude that the exact timing of simulated N₂O emissions at the onset of the B/A and the end of the YD depends sensitively on the climate forcing data prescribed in the LPX-Bern model.

6.1 Reconstructed terrestrial N₂O emissions and implications for the C–N cycle

Terrestrial N₂O emissions show a 40 % increase from the Last Glacial Maximum to the late preindustrial period (Fischer et al., 2019). Most of the deglacial increase was realized in two large steps, linked to rapid, decadal-scale, and widespread Northern Hemisphere warming and to shifts in the Intertropical Convergence Zone (ITCZ) and precipitation patterns. Each step was realized within a maximum of 200 years and thus on a timescale similar to that of the ongoing anthropogenic climate perturbation.

It has remained somewhat unclear whether warming will increase or reduce regional- to global-scale N₂O emissions, as responses to warming treatment are found to be highly variable across a range of conditions and ecosystems (Dijkstra et al., 2012). The ice-core record shows that past warming events, such as those at the onset of the B/A and the end of the YD, strongly promoted N₂O emissions globally. This suggests, as noted earlier (Schilt et al., 2014) and as projected (Stocker et al., 2013), that terrestrial N₂O emissions from natural land will likely increase further as climate warms, implying the existence of a positive climate feedback linked to the terrestrial N cycle (Xu-Ri et al., 2012).

The increasing terrestrial N₂O emissions imply either an increase in N₂O yield per unit N converted, averaged globally and across all N₂O production pathways, an increase in the global flux of converted N, or a combination of these factors. Regarding yield, uncertainties in the interpretation of the terrestrial N₂O emission record are linked to uncertainties in the ratio between N₂O produced to N converted. This yield factor is known to vary with environmental conditions and across N₂O production pathways, though quantitative experimental evidence is often not unequivocal (Diem et al., 2017; Davidson et al., 2000; Firestone and Davidson, 1989; Saggar et al., 2013; Butterbach-Bahl et al., 2002). For nitrification, the ratio of N₂O/${\mathrm{NO}}_{\mathrm{3}}^{-}$ produced increases with increasing acidity and decreasing oxygen (Firestone and Davidson, 1989) and increasing temperature (Smith, 1997). For denitrification, the production ratio of N₂O∕N₂ yield increases with increasing ${\mathrm{NO}}_{\mathrm{3}}^{-}$ availability, increasing oxygen concentration, decreasing decomposable carbon, and decreasing pH. Changes in precipitation may have altered the N₂O yield over the deglaciation. Higher soil water content and associated anoxic conditions generally favour the conversion of N₂O to N₂ and result in a low yield of N₂O (Weier et al., 1993), though dependencies of yield on water-filled pore space are sometimes complex (Diem et al., 2017) and soil texture and drainage affect water-filled pore space and yield (Bouwman et al., 2002). Yield is low under low ${\mathrm{NO}}_{\mathrm{3}}^{-}$ availability (Diem et al., 2017) and generally found to increase when ${\mathrm{NO}}_{\mathrm{3}}^{-}$ becomes more available (Weier et al., 1993). This relation, when considered in isolation, implies increasing N₂O emissions to be indicative of increasing ${\mathrm{NO}}_{\mathrm{3}}^{-}$ availability. The ratio of N₂O to N₂ emitted has been thought to decrease with increasing temperature (Firestone and Davidson, 1989). However, the observational challenges to determine N₂O and N₂ emissions are substantial, and available results are equivocal. Butterbach-Bahl et al. (2002), analysing results for a series of soil cores at spruce and beech sites in Germany, find it difficult to derive a general trend in the effect of temperature on the ratio of N₂O and N₂ produced via denitrification. Phillips et al. (2015) show for New Zealand pasture soil samples that denitrified N shifts towards N₂O with increasing temperature under anaerobic conditions. The warming over the glacial termination is also expected to accelerate organic matter turnover and thus availability of decomposable C, which would tend to lower this ratio (Firestone and Davidson, 1989). Further, experimental studies on grasslands yield both higher and lower ratios of N₂ to N₂O under elevated CO₂ compared to ambient CO₂ (Zhong et al., 2018; Baggs et al., 2003). In addition, different N₂O production pathways and their relative importance may evolve through time, and changes in soil organic matter decomposition pathways and in the stoichiometry of plant-derived material may affect yield. Given this complexity, we are not in a position to evaluate whether the N₂O yield factor for the combined N conversion processes increased or decreased over the deglaciation on the global scale. Changes in globally averaged yield may possibly explain the reconstructed terrestrial N₂O emission increase.

Regarding the flow of reactive N, the reactive N lost from ecosystems must be replaced if ecosystem N pools are not to be depleted. This mass balance consideration suggests that the total input flux of N changed hand in hand with related loss fluxes over the deglacial period and sufficiently to support changes in ecosystem N storage. An alternative to a balanced N input, N storage change, and output flux would be that N losses are fueled by existing reservoirs of reactive N. Large nitrate deposits (“caliche”) are currently found in the hyperarid Atacama desert (Ericksen, 1981; Tapia et al., 2018). But such large deposits are very unusual and require several hundred thousand years to accumulate (Michalski et al., 2004). Microorganisms require readily decomposable carbon substrate for the denitrification of nitrate; otherwise, nitrate may accumulate in the absence of substrate availability (Weier et al., 1993). Carbon substrate availability for denitrifiers might have increased at the onset of the B/A and at the end of the YD when climate warmed, precipitation patterns shifted, and organic matter remineralization accelerated in many regions. This in turn could have led to the conversion of nitrate that had potentially accumulated previously. This would lead to a depletion of this nitrate pool and, in turn, one would expect a decrease in N₂O emissions. Such a scenario appears in conflict with the reconstruction, given that reconstructed N₂O emissions remained high during the early B/A and even slightly increased throughout the Holocene.

Terrestrial N₂O emissions result primarily from nitrification and denitrification (Firestone and Davidson, 1989). A requirement for these processes to take place is the availability of ammonium and nitrate, as consumed by nitrifiers, denitrifiers, and plants alike. Variations in reactive N availability among diverse land classes are found to correspond to variations in N₂O emissions (Davidson et al., 2000). In addition, cryptogamic covers fix large amounts of N and contribute substantially to global N₂O emissions (Elbert et al., 2012). The ice-core data show that reactive N must have been available in sufficient quantities to support increasing terrestrial N₂O production over the deglacial period. Sources of reactive N on land, e.g. from BNF and weathering, may possibly have increased under warming climate and increasing CO₂ over the deglacial period and contributed to meet the N demand of plants, nitrifiers, denitrifiers, and cryptogamic covers under more favourable growth conditions. At the same time, the ice-core reconstructions do not suggest that reactive N for the production of N₂O remained scarce and that marine emissions dominated past atmospheric N₂O variations. The magnitude of possible adjustments in the N source is unclear. Wang and Houlton (2009) modelled an increase in N fixation under increasing CO₂ and temperature. However, their modelled increase in BNF was too small to meet projected N demand in global warming scenarios; the authors therefore proposed to downscale future C uptake in the warming projections of C4MIP models. However, this conclusion is based on an assumed optimum temperature for BNF of around 25 ^∘C, an assumption that has been challenged by Liao et al. (2017), who found the abundance of N-fixing trees to increase monotonically with temperature by analysing more than 125 000 forest plots in the USA and Mexico. More work is needed to improve projections of BNF under past and future environmental change, considering also the cost of BNF (Shi et al., 2018) and relying on data-based and model approaches (Fisher et al., 2012). Similarly, it remains unclear how other smaller sources of reactive N changed over the deglacial period and influenced N₂O emissions.

The ice-core data and emission estimates represent global values with century-scale resolution. They do not permit us to discriminate regional variations or interannual to decadal changes. In particular, N limitation is considered to be common in high-latitude regions. But this does not exclude a decadal- to century-scale increase in N sources under changing environmental conditions. Sufficient N and other nutrients were available to support the northward expansion of boreal forest during the deglaciation. Indeed, it has been hypothesized that an increase in BNF may become cost-effective under more severe N limitation at mid and high latitudes (Menge et al., 2017). Field data support a role of BNF in supporting forest growth, in particular for early successional forests. Experimental evidence comes from chronosequence studies in the Amazon region and in Panama (Batterman et al., 2013; Gehring et al., 2005), from forest clay and sand box studies in Brazil and in the USA (Davidson et al., 2018; Bormann et al., 2002), as well as from forest inventory data covering tropical-to-temperate climatic conditions (Liao et al., 2017). Regional changes may have contributed to the deglacial N₂O emission increase.

6.2 Simulating deglacial terrestrial N₂O emission changes with LPX-Bern

We applied the LPX-Bern dynamic global vegetation model to investigate regional N₂O emissions patterns, governing mechanisms, and C–N coupling in simulations over the past 21 000 years. The model represents the cycling of N and C in soils and plants in a simplified way. C and N are stored in PFT-specific plant and litter pools and a fast and a slow overturning soil pool in each grid cell. This chain of pools represents a spectrum of overturning timescales on each grid cell. LPX-Bern, being a spatially explicit model, reacts differently in different regions to changes in climatic and environmental drivers (Figs. 4, 8). However, microbial and fungal biomasses are – unlike in microbial-explicit models (Schimel and Weintraub, 2003; Zhu et al., 2017; Allison and Gessner, 2012) – not explicitly modelled, and organic matter decomposition does not depend on microbial mass and physiology. Instead, a mass balance approach is applied with C:N stoichiometry prescribed at observation-based PFT-dependent values for litter and soils. A constant fraction of remineralized N is returned immediately to its source soil pool and the N budget is closed by the implied N source flux (Fig. 1). There is also no distinction between different classes of organic matter according to their accessibility to microbial action (Averill and Waring, 2018). Potential N₂O release by lichens and bryophytes is also not modelled (Porada et al., 2017). Only net N₂O production during denitrification is considered, while gross production and consumption of N₂O during denitrification (Schmid et al., 2001b) are not explicitly modelled. In future work, N₂O may be treated as an intermediate product of denitrification to explicitly simulate N₂O consumption as well as the related ratio of N₂O to N₂ emissions, i.e. the yield factor. In LPX-Bern v1.4N, the yield factor varies as a function of temperature for nitrification and denitrification, whereas dependencies of the yield on ${\mathrm{NO}}_{\mathrm{3}}^{-}$ availability, oxygen status, decomposable carbon, or pH are not considered (see the discussion in the previous section).

Uncertainties in modelled N₂O emissions are also associated with uncertainties in the representation of N sources and N loss processes. The N source and its changes in LPX are implied by maintaining the soil C:N ratio at observed PFT-dependent values. This N source thereby accounts for any other source, except explicitly prescribed N deposition (Fig. 1). Carbon costs of N acquisition are not directly considered. However, a fraction of 6 % of NPP is directly transferred to a pool with a short overturning time to represent root exudates. The largest contribution by far to the implied N source is thought to come from BNF.

Different approaches to represent BNF are used in different models, and their effects on modelled BNF, NPP, carbon stocks, and N₂O emissions are compared in recent studies for modern conditions and future projections (Meyerholt et al., 2016; Wieder et al., 2015a). The simplest approach is a prescribed static map of BNF (Zaehle and Friend, 2010). Most frequently, empirical models are used, describing BNF as a linear function of evapotranspiration (Cleveland et al., 1999) or as an exponential function of NPP (Thornton et al., 2007). More process-oriented models heuristically account for the dependency of symbiotic BNF on N demand by vegetation, soil N status, and light limitations in extra-tropical regions (Gerber et al., 2010) or on the optimization of plant C investment in resource acquisition (Fisher et al., 2010). Meyerholt et al. (2016) describe asymbiotic BNF as a function of temperature, shading, and soil moisture.

The implementation of such different parameterizations in LPX-Bern would likely lead to different estimates for deglacial changes in BNF and N₂O emissions. We may expect that prescribing a constant modern BNF field would lead to approximately constant fluxes for nitrification, denitrification, and leaching. In turn, deglacial changes in N₂O emissions would only be due to changes in yield factors and be smaller than modelled in our standard simulation and smaller than reconstructed. A similar result is expected for BNF depending on evapotranspiration, because globally averaged evapotranspiration changes little over the deglaciation in our standard simulation (Table A1) and global BNF would remain at its modern value. In contrast, some parameterizations would likely yield similar or larger responses than simulated here. A strong increase in BNF is found in global warming simulations for the N demand (Meyerholt et al., 2016) and the NPP-based (Wieder et al., 2015a) BNF parameterizations. In this latter parameterization, BNF responds immediately to NPP and grows exponentially with NPP, whereas the deglacial increase in BNF of 11 % is smaller than in NPP (13 %) in our standard run (Table A1). Overall, one might expect a similar or even larger increase in BNF and N₂O emissions over the deglaciation when replacing our implicit N source approach with a demand, NPP, or cost-driven BNF parameterization. However, a corresponding quantitative analysis is beyond the scope of this study and is left to future efforts.

The N loss in LPX is predominantly driven by local gaseous loss from denitrification, with a much smaller role for fire, leaching, minor contributions from volatilization of NH₃, and gas release during nitrification (Fig. 1). Denitrification and nitrification are thought to occur at anaerobic and aerobic microsites in the soil, which are challenging to represent. In LPX-Bern the fraction of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ available for aerobic nitrification and of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ available for anaerobic denitrification within a grid cell is assumed to scale linearly with water-filled pore space in the top 50 cm (Xu-Ri and Prentice, 2008). Meyerholt and Zaehle (2018) investigated different algorithms for N loss processes and find variable responses in N loss and in the partitioning of N losses between gaseous and leaching losses under elevated CO₂.

LPX-Bern was applied in an earlier study to simulate climate–N₂O feedbacks under global warming (Stocker et al., 2013), and results for global emissions for the period from 16 to 10 ka are presented in Schilt et al. (2014). Recently, model parameters were updated using a set of modern observational constraints in a Bayesian approach (version v1.4) (Lienert and Joos, 2018a) and here further modified towards a lower N source (version v1.4N) by adjusting the fraction of N remineralized that is returned to the soil pool. Except for the matching of preindustrial N₂O emission estimates, the model has not been tuned in any way towards matching the ice-core reconstruction. These updated versions also account for the uptake of mineral N by N immobilization in soils. Results in terms of deglacial N₂O emission changes are similar between these different versions. Changes in parameter values and selected globally averaged model results for version v1.4 versus v1.4N are documented in Appendix Table A1. Appendix Figs. A1 and A2 show the spatial patterns for the implied N source, soil mineral N, NPP, total carbon, and their changes over the deglaciation. The most striking difference is that the implied N source is a factor of 4 higher in v1.4 compared to v1.4N. Nevertheless, results for deglacial N₂O emission changes remain similar between versions v1.4 (change 1.13 Tg N yr⁻¹) and v1.4N (0.96 Tg N yr⁻¹). The reason is that the relative, percentage changes in the N source flux and in N loss fluxes are similar between the two versions; however, the N₂O yield factors are calibrated in both versions to the same preindustrial N₂O emissions of 5.9 Tg N yr⁻¹. The lower N source leads, however, to smaller mineral N concentration. In turn, N limitation of NPP is larger and, hence, NPP and carbon storage are smaller, both in steady state and for deglacial change, in v1.4N than in v1.4.

LPX-Bern forced by TraCE-21ka output simulates a deglacial increase in N₂O emission but underestimates the reconstructed increase by about 50 %. This may be related to the model's C cycle and/or to the model formulation for denitrification and other shortcomings discussed above. LPX-Bern is known for a relatively low sensitivity to increasing CO₂ and simulates a modest increase in NPP and the terrestrial carbon sink over the industrial period (Lienert and Joos, 2018a). A larger deglacial increase in NPP would tend to increase the implied N source in the model and potentially increase nitrification, denitrification, and N₂O emissions. Further, denitrification processes are assumed to respond to the relevant substrate concentration following Michaelis–Menten kinetics (Xu-Ri and Prentice, 2008). This limits N conversion rates and N₂O production at high substrate concentrations, whereas a recent synthesis of N fertilization experiments (Niu et al., 2016) points towards an exponential relationship between N load and N₂O emissions. Regarding the timescales of response, results of a sensitivity simulation demonstrate that LPX-Bern is able to represent abrupt, multi-decadal-scale change in global terrestrial N₂O emissions as reconstructed for past rapid warming events. Further, factorial simulations demonstrate a complex interplay of different driving factors for deglacial N₂O emissions with considerable synergistic or antagonistic interactions related to changes in temperature, precipitation, and CO₂. Remarkably, about a third of the simulated global increase in terrestrial N₂O emissions over the past 21 000 years was counteracted by land loss due to flooding of previously exposed shelf regions.

LPX-Bern simulates a slight decrease in N₂O emissions over the Holocene period, whereas the reconstruction shows an increase by about 0.4 Tg N yr⁻¹ over this period. The exact reason for this discrepancy between modelled and reconstructed change is unclear. The model results suggest that the global N₂O emission change results from the difference of, in some regions, large positive minus, in other regions, large negative emission changes (Fig. 4c) and that these regional emission changes depend sensitively on regional changes in precipitation and temperature (Fig. 8). Thus, it may remain somewhat challenging to simulate the balance between positive and negative emission changes correctly, given uncertainties in land models and the climatic drivers as well as in the deconvolution of the ice-core concentration and isotope data.

We applied LPX-Bern to analyse the consequences of alternative C–N cycle hypotheses. In the model's standard setup, the N cycle responds dynamically to changes in environmental conditions; the N source increases under increasing ecosystem N demand. This additional N source contributed to meet the rising N demand of plants simulated under more favourable growth conditions and resulted in a higher flow of N through land ecosystems and increased N₂O production from terrestrial ecosystems as well as increased carbon storage on land. The increase in N₂O emissions, linked to this increase in N source, is further amplified by an increase in the yield of N₂O per N converted by nitrification and denitrification.

In an alternative deglacial simulation, the sources of reactive N were prescribed at the level simulated for the Last Glacial Maximum. In this case, the model yields a decrease in N₂O emissions, in clear conflict with the ice-core data. The increase in N₂O emissions attributable to the increase in the model's N₂O yield factors for denitrification and nitrification under deglacial warming as well as the positive response of soil and litter turnover rates to increasing temperatures are not sufficient to offset the decrease in N₂O emissions in this run with a constant N source. The simulated terrestrial C inventory decreased over the deglaciation in this simulation, again in conflict with data-based reconstructions (Lindgren et al., 2018; Jeltsch-Thömmes et al., 2019) that suggest a growing land biosphere carbon inventory as also predicted in the standard setup.

Within the current setup of LPX-Bern, acknowledging its limitations, a growing N source appears necessary to simulate increasing N₂O emissions over the last deglaciation. It would be interesting to see whether other models incorporating the cycling of N and C (e.g. Zaehle et al., 2014; Averill and Waring, 2018; Thornton et al., 2007; Fisher et al., 2010) and N₂O emissions (Zhu et al., 2016; Saikawa et al., 2013; Tian et al., 2015; Xu et al., 2017; Huang and Gerber, 2015; Zaehle and Friend, 2010; Olin et al., 2015; Goll et al., 2017) reveal similar or alternative mechanisms and whether these models are able to represent the reconstructed terrestrial N₂O emissions over the past 21 000 years. Future research efforts are needed to test the robustness of our results and to further evaluate the potential role of alternative or complementary explanations. In particular, the significantly smaller amplitudes of our modelled N₂O changes compared to the data reconstruction justify further research. The setup of our deglacial simulation could be used to compare and evaluate models, e.g. in the framework of the model intercomparison initiated by the Global Carbon Project (Tian et al., 2018).