N 2 O changes from the Last Glacial Maximum to the preindustrial – Part 2: terrestrial N 2 O emissions and carbon–nitrogen cycle interactions

. Carbon–nitrogen (C–N) interactions regulate N availability for plant growth and for emissions of nitrous oxide (N 2 O) and the uptake of carbon dioxide. Future projections of these terrestrial greenhouse gas ﬂuxes are strikingly divergent, leading to major uncertainties in

Abstract. Carbon-nitrogen (C-N) interactions regulate N availability for plant growth and for emissions of nitrous oxide (N 2 O) and the uptake of carbon dioxide. Future projections of these terrestrial greenhouse gas fluxes are strikingly divergent, leading to major uncertainties in projected global warming. Here we analyse the large increase in terrestrial N 2 O emissions over the past 21 000 years as reconstructed from ice-core isotopic data and presented in part 1 of this study. Remarkably, the increase occurred in two steps, each realized over decades and within a maximum of 2 centuries, at the onsets of the major deglacial Northern Hemisphere warming events. The data suggest a highly dynamic and responsive global N cycle. The increase may be explained by an increase in the flux of reactive N entering and leaving ecosystems or by an increase in N 2 O yield per unit N converted. We applied the LPX-Bern dynamic global vegetation model in deglacial simulations forced with Earth system model climate data to investigate N 2 O emission patterns, mechanisms, and C-N coupling. The N 2 O emission changes are mainly attributed to changes in temperature and precipitation and the loss of land due to sea-level rise. LPX-Bern simulates a deglacial increase in N 2 O emissions but underestimates the reconstructed increase by 47 %. Assuming timeindependent N sources in the model to mimic progressive N limitation of plant growth results in a decrease in N 2 O emissions in contrast to the reconstruction. Our results appear consistent with suggestions of (a) biological controls on ecosystem N acquisition and (b) flexibility in the coupling of the C and N cycles during periods of rapid environmental change. A dominant uncertainty in the explanation of the reconstructed N 2 O emissions is the poorly known N 2 O yield per N lost through gaseous pathways and its sensitivity to soil conditions. The deglacial N 2 O record provides a constraint for future studies.

Introduction
Nitrous oxide (N 2 O) is a sensitive proxy of biogeochemical and ecosystem processes on land and in the ocean, and its past atmospheric variations are recorded in ice cores. N 2 O is an important greenhouse gas and contributes to ongoing global warming . It is also involved in the destruction of stratospheric ozone (Myhre et al., 2013). N 2 O is produced primarily by nitrification and denitrification both on land and in the ocean and photochemically decomposed in the stratosphere (Ciais et al., 2013). Atmospheric N 2 O increased from 270 ppb (MacFarling Meure et al., 2006) to around 330 ppb (https://www.esrl.noaa.gov/gmd/, last access: 30 June 2020) over the industrial period due to human activities, including fertilizer application, fossil fuel use, and biomass burning Ciais et al., 2013). Atmospheric N 2 O varied naturally between around 180 and 300 ppb over glacial-interglacial cycles (Sowers et al., 2003;Spahni et al., 2005;Schilt et al., 2010). A quantitative explanation of these variations is lacking, and this knowledge gap renders projections of the feedbacks between N 2 O and climate change uncertain Battaglia and Joos, 2018;Kracher et al., 2016).
Variations in N 2 O emission, and thus in tropospheric N 2 O content, are closely linked to ecosystem processes governing the cycling of nitrogen (N) and carbon (C) on land and in the ocean (Gruber and Galloway, 2008). N availability to support land plant and phytoplankton growth and terrestrial and marine C storage is governed by the balance of N input and loss fluxes. Reactive N is added by biological nitrogen fixation (BNF) and, on land, by deposition and weathering. Reactive N is lost from ecosystems through nitrification followed by denitrification, as well as leakage and mineral adsorption. Terrestrial N 2 O emissions are a sensitive indicator of the flow of reactive N entering and leaving land ecosystems (Firestone and Davidson, 1989), and reactive N needs to be available for the production of N 2 O as well as to support carbon uptake and storage by plants and the land biosphere.
The land biosphere sequesters about a quarter of anthropogenic CO 2 emissions and is the largest natural N 2 O source (Ciais et al., 2013). However, key features of the C-N cycle are poorly understood, leading to major uncertainties in global warming projections (Joos et al., 2001;Arora et al., 2013;Friedlingstein et al., 2006Friedlingstein et al., , 2013Plattner et al., 2008;Zickfeld et al., 2013). The question to what extent N availability will limit future land C uptake and N 2 O emissions is unresolved. In particular, large uncertainties remain as to what extent N sources, net N mineralization, N uptake, and N loss processes will adjust and influence plant growth and N 2 O emissions under climatic and environmental change (Niu et al., 2016). It is debated how global warming, increasing CO 2 , and increased N deposition affect the current land carbon sink (Körner, 2015;Fatichi et al., 2014;Terrer et al., 2016) and how the land carbon sink and greenhouse gas emissions will evolve. Apparently conflicting observa-tions (Davidson et al., 2007;Luo et al., 2004;Reich et al., 2014;Vitousek et al., 2013;Xu-Ri and Prentice, 2008), theories (Zhu et al., 2017;Menge et al., 2017) and model projections (Hungate et al., 2003;Todd-Brown et al., 2013;Wieder et al., 2015b) of the role of N limitation in plant growth and the land C sink (Meyerholt et al., 2016;Zaehle et al., 2014;Walker et al., 2015) represent a major uncertainty in future projections of atmospheric CO 2 and climate Joos et al., 2013;Todd-Brown et al., 2013). Results from free-air CO 2 enrichment (FACE) and open-top chamber experiments show that N limitation -at least in some ecosystems -reduces the response of aboveground biomass growth to elevated CO 2 (McMurtrie et al., 2008;Norby et al., 2010;Terrer et al., 2016Terrer et al., , 2019 and that soil N and P availabilities are important controls on this magnitude in ectomycorrhizal and arbuscular mycorrhizal systems, respectively. Increased N immobilization in plant litter, biomass and soil (Luo et al., 2004), and multiple nutrient limitation (Körner, 2015;Vitousek et al., 2010;Elser et al., 2007) have been put forward as explanations for reduced growth responses.
Field and laboratory experiments provide important insights, but extrapolation of such results on the short-term response in C and N fluxes to the multi-decadal to century timescales, relevant for global warming projections, is uncertain. Additional hindrances to the improvement of our quantitative understanding of C-N cycle coupling are related to large variations in fluxes and inventories on small spatial and temporal scales (Arias-Navarro et al., 2017;Barton et al., 2015) and the diversity of responses across different organisms and ecosystems. Largely missing are long-term observational constraints on the global coupled C and N cycle during periods of rapid climate change and increasing atmospheric CO 2 .
High-resolution data on the isotopic composition of N 2 O from Antarctic ice cores have the potential to provide precise information on past variations in terrestrial and marine N 2 O emissions and thus on C-N coupling on timescales from decades to many centuries. A recent ice-core study on the stable isotope composition of N 2 O demonstrated the power of this approach . The isotopic data showed that both land and oceanic sources increased during the interval from 16 000 to 10 000 years before present (ka), when ocean circulation and climatic changes strongly affected the global cycling of both CO 2 and N 2 O . Schilt et al. concluded that natural N 2 O emissions will probably increase in response to global warming.
In part 1 of this study , this earlier work was extended to reconstruct the evolution in terrestrial versus oceanic emissions of N 2 O from the Last Glacial Maximum (LGM; 21 ka) to the late preindustrial Holocene using novel, high-precision N 2 O stable isotope data. The age scales of the records were aligned with the absolutely counted GICC05 age scale (Rasmussen et al., 2006) using high-resolution methane measurements leading to particularly small age-scale uncertainties around past rapid warming events, namely, at the onset of the Bølling-Allerød Northern Hemisphere warm period and of the Holocene interglacial. The ice-core data of N 2 O concentration and δ 15 N(N 2 O) were used in a model deconvolution to determine changes in global marine and terrestrial N 2 O emissions. The results indicate that N 2 O emissions from land and ocean increased over the deglaciation largely in parallel by 1.8±0.3 and 0.7± 0.3 Tg N yr −1 , respectively, relative to the Last Glacial Maximum level. However, during the abrupt Northern Hemisphere warming events, terrestrial emissions changed more rapidly than marine emissions. The ice-core data reveal large steplike changes in terrestrial emissions by up to 1 Tg N yr −1 at the onset of the warming events, realized over decades and, given the proxy data resolution, last 200 years at maximum. The multi-decadal timescales of these events are also relevant for 21st century climate projections and much longer than accessible in typical laboratory or field experiments. Further, the abrupt N 2 O emission changes indicate a very high sensitivity of the terrestrial C-N cycle to environmental change. But a detailed process-based investigation of terrestrial N 2 O emission changes over the deglaciation, rapid past warming events, and the Holocene warm period, and their links to the flow of N and C in land ecosystems, has not been attempted before.
The aim of part 2 of this study is to improve our understanding of the cycles of N 2 O, C, and N. The unique, new terrestrial N 2 O emission record of the past 21 000 years is used to explore mechanisms of the functioning of the C-N cycle on land and to quantify terrestrial drivers for atmospheric N 2 O concentration changes. Controls on variations in terrestrial emissions are elucidated in the framework of a dynamic global vegetation model. Part 1 of this study  presents the ice-core N 2 O concentration and isotope data and the reconstruction of global terrestrial and marine N 2 O emissions for the past 21 000 years and provides a discussion of the reconstructed marine emissions in the context of past climate and oceanographic changes. A model-based interpretation of the reconstructed marine N 2 O emission changes during past abrupt climate events is given by Joos et al. (2019). Here in part 2, we focus on the interpretation of the terrestrial N 2 O emission record and discuss terrestrial N 2 O emissions and C-N responses in transient deglacial simulations with a dynamical vegetation/terrestrial biogeochemistry model.  (Battaglia and Joos, 2018), in line with IPCC AR5 estimates of 6.6 (3.3 to 9.0) Tg N yr −1 . N 2 O is produced through a variety of pathways both in the ocean and on land and N 2 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 2 O production are thought to be respiratory denitrification ( , under low oxygen conditions, and autotrophic nitrification (NH 3 → NO − 2 → NO − 3 ), mediated by archaea, bacteria, and fungi. In addition, heterotrophic nitrification, which is the oxidation of organic N to nitrite (NO − 2 ) and subsequent reduction to N 2 O by incomplete denitrification, was found to be the dominant path in a grassland ecosystem after fertilizer application (Moser et al., 2018). N 2 O is produced as a byproduct of nitrification and as an intermediate product during denitrification. N 2 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 2 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 2 O in moist tropical areas and lower emissions at high latitudes Zhuang et al., 2012;Werner et al., 2007;Potter et al., 1996;Xu-Ri et al., 2012;Wells et al., 2018). Moist soils typically show relatively high N 2 O emission (Butterbach-Bahl et al., 2013;Zhuang et al., 2012), but the environmental dependencies of N 2 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 2 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 2 emissions from soil organic matter decomposition, and a positive feedback of soil warming on N 2 O emissions is expected (Butterbach-Bahl et al., 2013). Warming treatment increased measured N 2 O emissions in boreal peatlands (Cui F. Joos et al.: N 2 O changes from the Last Glacial Maximum to the preindustrial -Part 2 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 2 O emissions (Dijkstra et al., 2012). Positive or neutral responses in N 2 O emissions have been found in field experiments under elevated CO 2 in temperate and boreal forests and grasslands 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 2 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 2 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 2 O production pathways and N loss processes and their relative importance may evolve through time and influence the N 2 O yield on local to global scales.

C-N-N 2 O coupling
The main flow paths of reactive N and its link to N 2 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 NH + 4 and NO − 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 NO − 3 pool is replenished by nitrification, the conversion of NH + 4 to NO − 3 . We note the acid-base equilibrium between NH + 4 and NH 3 in soil water; for simplicity, we generally refer to NH + 4 only. Turning to the external cycle, reactive N enters land ecosystems (blue arrows in Fig. 1) through the conversion of dinitrogen (N 2 ) to organic N and eventually to NH + 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 2 O), leaching of NO − 3 and dissolved organic matter by runoff, and emissions of N compounds by fire Ciais et al., 2013;Hedin et al., 1995Hedin et al., , 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 2 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) andNiu et al. (2016), shows that the uptake of N by plants, net primary productivity, and biomass as well as NH x and NO − 3 pools in soils, nitrification, nitrate leaching, denitrification, and N 2 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 2 O. N 2 O production on land is predominantly associated with denitrification and to a smaller extent with nitrification. Large gaseous losses and N 2 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 2 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 2 O yield per unit reactive N converted in land ecosystems.  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).
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 2 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 codetermines 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 −1 . 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 −1 . Cleveland et al. (1999) used 100 plot-scale estimates of BNF to estimate global BNF on natural ecosystems to 195 Tg N yr −1 (range: 100 to 290 Tg N yr −1 ). Vitousek et al. (2013) suggest a plausible range for preindustrial BNF of 40 to 100 Tg N yr −1 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 −1 , based on cost-benefit modelling for N fixation. The same authors estimate asymbiotic N fixation to 22 Tg N yr −1 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 −1 (27-99 Tg N yr −1 ) by integrating experimental data from 200 studies. In addition, rock weathering is estimated to add 15 (10 to 20) Tg N yr −1 to land ecosystems (Houlton et al., 2018). Xu-Ri and Prentice (2017) estimate global N sources to 340 (230-470) Tg N yr −1 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 −1 .
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 NH + 4 and NO − 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 NH + 4 and NO − 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 NH + 4 soil pool with a temperature-dependent rate coefficient. Denitrification is modelled as a two-step process whereby NO − 3 is converted to nitrite (NO − 2 ) and NO − 2 is further converted to N 2 , following Michaelis-Menten kinetics with dependence on the substrates NO − 3 and NO − 2 as well as on labile C availability. Reaction rates are temperature dependent. NO − 3 leaching depends on soil NO − 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 3 volatilization.
N 2 O emissions are computed by assuming that fractions of the N fluxes associated with denitrification, nitrification, and N leaching are released as N 2 O. The calibration of the denitrification and nitrification yield factors is done to satisfy the constraints from the global atmospheric N 2 O budget, which suggests global terrestrial emissions to be 5.9 Tg N yr −1 at preindustrial (Battaglia and Joos, 2018). Since global terrestrial N 2 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 codetermine these rates, and then calibrated the yield parameters, given the N loss rates. For denitrification, the globally dominant N 2 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 2 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 2 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 2 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 2 O per N 2 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 2 O yield, and preindustrial N 2 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 −1 . 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 2 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 −1 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 2 O emissions is only 0.2 Tg N yr −1 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 2 O emissions are not sensitive to the parameter settings for the yield and the flux representing immobilization.

Setup for transient glacial-interglacial simulations
A previously described LPX-Bern model setup for glacialinterglacial simulations is applied 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 2 (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 observationbased, 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).
The same model parameter values as determined by Lienert and Joos (2018b) are used, except for the immobilization fractions and N 2 O yield factors (Table A1). Regarding the N module, the parameters are the same as in previous studies Schilt et al., 2014) addressing N 2 O emissions over the past and under future global warming. An exception is an adjustment in the upper limit of the fraction of NH + 4 nitrified per day from 0.1 to 0.09096 d −1 at 20 • C and in the yield factors for denitrification, nitrification, and leaching in response to new observation-constrained es-timates of marine N 2 O emissions and to a slightly revised estimate of the atmospheric N 2 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 2 O emissions and other model outcomes are attributed to individual driving factors (temperature, precipitation, CO 2 , 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 2 O emissions ( eN 2 O) simulated in the baseline run and the sum of the emission changes attributed to individual drivers: The dominant driver is identified as having the largest contribution to eN 2 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 (| eN 2 O|<1 mg N m −2 yr −1 ) 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 2 are abruptly changed to conditions for 2.5 ka. The run is continued for another 1500 years with climate and CO 2 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 2 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).

Reconstructed terrestrial N 2 O emissions and implications
We start the presentation of results by summarizing the main feature of the terrestrial N 2 O emission record (Fig. 3, green line) presented in part 1 of this study . In part 1, the global N 2 O emissions from land and from the ocean are jointly reconstructed by deconvolving novel icecore data of N 2 O and of its isotopic signature, δ 15 N(N 2 O), using an established method and relying on differences in the isotopic signature of land versus marine N 2 O emissions. Terrestrial emissions increased between LGM (21 ka) and PI (1500 CE) by about 1.8 Tg N yr −1 . Terrestrial emissions re- mained approximately invariant during the Heinrich Stadial I Northern Hemisphere (NH) cold phase (HS1; 17.4 to 14.6 ka, Rasmussen et al., 2014) until the start of the B/A NH warm period (14.6 to 12.8 ka). Then, land emissions increased at the start of the B/A, declined again during the YD NH cold period (12.8 to 11.7 ka), and peaked at the start of the preboreal period (PB), followed by a modest increase during the Holocene.
Remarkably, the overall deglacial increase in terrestrial N 2 O emissions was mainly realized in two large steps at the onset of the B/A and at the end of the YD, two major Northern Hemisphere warming events. The detailed analysis of the ice-core N 2 O concentration and isotope data (see  reveals that global terrestrial N 2 O emissions started to rise at the beginning of the warming events. The enclosure process of atmospheric air into firn and ice acts like a low-pass filter, smoothing any fast variations in atmospheric N 2 O, in its isotopic signature, and, correspondingly, in inferred emissions. Fischer et al. conclude that global terrestrial N 2 O emissions reacted within a maximum of 200 years to the large-scale climate reorganizations associated with the two major deglacial Northern Hemisphere warming events. Overall, the ice-core data show that land ecosystem N 2 O emissions responded sensitively to climatic and environmental changes over the deglaciation. The rapid increase in ter-restrial N 2 O emissions at the onset of the B/A and at the end of the YD and the overall increase in emissions over the past 21 000 years either point (i) to an increase in N 2 O yield per unit N converted for emissions to the atmosphere, averaged globally and across all N 2 O production pathways, or/and (ii) to an increase in the global flux of converted N. Reactive N was available in sufficient amounts to support an increase in global terrestrial N 2 O emissions during periods where environmental conditions became, on a global scale, more favourable for plant growth and C sequestration (Ciais et al., 2012;Bird et al., 1994;Joos et al., 2004;.  deglacial increase in N 2 O emissions is only 0.96 Tg N yr −1 and about 47 % smaller than the reconstructed increase of 1.8 ± 0.3 Tg N yr −1 (±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 2 O emissions over the last 11 000 years, the Holocene period. Simulated terrestrial N 2 O emissions decreased by 0.06 Tg N yr −1 between 10 and 0.5 ka, whereas reconstructed emissions increased by 0.41 Tg N yr −1 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.
The changes in global terrestrial N 2 O emissions are the result of spatially differentiated responses. LPX-Bern simulates high natural emissions of N 2 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 2 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 −2 yr −1 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 2 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 2 O emissions between 15 and 8 ka. The net loss of terrestrial N 2 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 2 O emissions. Changes in the land extent caused by changes in sea level represent an important factor for past global terrestrial N 2 O emissions.
The patterns of change in terrestrial N 2 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 2 O emissions during the Holocene, LPX-Bern simulated large re- gional shifts in source strength, linked to changes in temperature and precipitation. This includes for example a decrease in N 2 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.

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 2 O emission change in the model: (i) changes in the source of reactive N, fueling nitrification and denitrification, and (ii) changes in the N 2 O yield per unit N converted.
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 −2 yr −1 at northern mid and high latitudes and around 1 g N m −2 yr −1 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 −2 yr −1 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 −2 in tropical forests and around 1 g N m −2 in temperate and boreal forests and tundra regions, while higher values of up to and more than 10 g N m −2 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 2 O yield factors, i.e. the N 2 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 2 O yield in the standard compared to this sensitivity run, and 0.44 Tg N yr −1 of the deglacial increase in land N 2 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 2 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 −1 at the LGM to 128 Tg N yr −1 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 2 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 2 O emissions (Fig. 3, black line). The decrease in global land N 2 O emissions is related to the loss of land in response to sea-level rise (−0.11 Tg N yr −1 ) and a small decrease in N 2 O emissions on remaining land (−0.088 Tg N yr −1 ) 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 −1 ) instead of an increase (5.9 Pg C yr −1 ) 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 nonflooded 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 −1 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 −1 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 2 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 extratropics. 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 2 O emissions from the land biosphere. The increases in N 2 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 2 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 2 O emissions. If the model is al-lowed 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 2 O emissions increase substantially.

Attribution of simulated terrestrial N 2 O emissions to climatic and environmental drivers
The complex spatio-temporal changes in land N 2 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 −1 to the overall emission increase, followed by changes in precipitation (0.7 Tg N yr −1 , Fig. 7). Changes due to increasing atmospheric CO 2 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 2 O emissions (−0.5 Tg N yr −1 ). 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). 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 2 O emissions, which occurred earlier in the simulation than the ice-core reconstruction, is not driven by CO 2 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 2 , 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.
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 2 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 2 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 2 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 2 in this region (Fig. 8e, f) most likely reflect non-linear interactions between precipitation, temperature, and CO 2 . 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 2 generally play a secondary role.

The timescales of response to a step change in climate and CO 2
Next, we investigate potential causes of the relatively slow N 2 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 2 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 2 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 2 ) 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 2 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 2 O of more than 100 years into account, such a fast multi- 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 NO − 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 NO − 3 , leading to the decrease in N leakage, while warming accelerates soil decomposition and the release of NH + 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 2 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 2 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 2 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. . 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 2 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 Figure 8. Attribution of simulated global terrestrial N 2 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 2 (e, f), and non-linear interactions among the drivers plus the small effect due to changes in incident solar radiation (g, h). of the YD, strongly promoted N 2 O emissions globally. This suggests, as noted earlier  and as projected , that terrestrial N 2 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 2 O emissions imply either an increase in N 2 O yield per unit N converted, averaged globally and across all N 2 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 2 O emission record are linked to uncertainties in the ratio between N 2 O produced to N converted. This yield factor is known to vary with environmental conditions and across N 2 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 2 O/NO − 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 2 O/N 2 yield increases with increasing NO − 3 availability, increasing oxygen concentration, decreasing decomposable carbon, and decreasing pH.  Changes in precipitation may have altered the N 2 O yield over the deglaciation. Higher soil water content and associated anoxic conditions generally favour the conversion of N 2 O to N 2 and result in a low yield of N 2 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 . Yield is low under low NO − 3 availability (Diem et al., 2017) and generally found to increase when NO − Figure 11. Response to a step-like change in (a) environmental conditions for globally averaged terrestrial (b) N 2 O emissions and (c) NPP. Atmospheric CO 2 (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. when considered in isolation, implies increasing N 2 O emissions to be indicative of increasing NO − 3 availability. The ratio of N 2 O to N 2 emitted has been thought to decrease with increasing temperature (Firestone and Davidson, 1989). However, the observational challenges to determine N 2 O and N 2 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 2 O and N 2 produced via denitrification. Phillips et al. (2015) show for New Zealand pasture soil samples that denitrified N shifts towards N 2 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 decompos-able 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 2 to N 2 O under elevated CO 2 compared to ambient CO 2 (Zhong et al., 2018;Baggs et al., 2003). In addition, different N 2 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 2 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 2 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 2 O emissions. Such a scenario appears in conflict with the reconstruction, given that reconstructed N 2 O emissions remained high during the early B/A and even slightly increased throughout the Holocene.
Terrestrial N 2 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 2 O emissions (Davidson et al., 2000). In addition, cryptogamic covers fix large amounts of N and contribute substantially to global N 2 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 2 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 2 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 2 O remained scarce and that marine emissions dominated past atmospheric N 2 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 2 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 2 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 costeffective under more severe N limitation at mid and high latitudes . 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 . Regional changes may have contributed to the deglacial N 2 O emission increase.

Simulating deglacial terrestrial N 2 O emission changes with LPX-Bern
We applied the LPX-Bern dynamic global vegetation model to investigate regional N 2 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 2 O release by lichens and bryophytes is also not modelled (Porada et al., 2017). Only net N 2 O production during denitrification is considered, while gross production and consumption of N 2 O during denitrification (Schmid et al., 2001b) are not explicitly modelled. In future work, N 2 O may be treated as an intermediate product of denitrification to explicitly simulate N 2 O consumption as well as the related ratio of N 2 O to N 2 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 NO − 3 availability, oxygen status, decomposable carbon, or pH are not considered (see the discussion in the previous section). Uncertainties in modelled N 2 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 2 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 2 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 2 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 2 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 3 , 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 NH + 4 available for aerobic nitrification and of NO − 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 2 .
LPX-Bern was applied in an earlier study to simulate climate-N 2 O feedbacks under global warming , 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 2 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 2 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 2 O emission changes remain similar between versions v1.4 (change 1.13 Tg N yr −1 ) and v1.4N (0.96 Tg N yr −1 ). 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 2 O yield fac-3530 F. Joos et al.: N 2 O changes from the Last Glacial Maximum to the preindustrial -Part 2 tors are calibrated in both versions to the same preindustrial N 2 O emissions of 5.9 Tg N yr −1 . 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 2 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 2 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 2 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 2 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 2 O emissions. Regarding the timescales of response, results of a sensitivity simulation demonstrate that LPX-Bern is able to represent abrupt, multidecadal-scale change in global terrestrial N 2 O emissions as reconstructed for past rapid warming events. Further, factorial simulations demonstrate a complex interplay of different driving factors for deglacial N 2 O emissions with considerable synergistic or antagonistic interactions related to changes in temperature, precipitation, and CO 2 . Remarkably, about a third of the simulated global increase in terrestrial N 2 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 2 O emissions over the Holocene period, whereas the reconstruction shows an increase by about 0.4 Tg N yr −1 over this period. The exact reason for this discrepancy between modelled and reconstructed change is unclear. The model results suggest that the global N 2 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 increas-ing 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 2 O production from terrestrial ecosystems as well as increased carbon storage on land. The increase in N 2 O emissions, linked to this increase in N source, is further amplified by an increase in the yield of N 2 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 2 O emissions, in clear conflict with the ice-core data. The increase in N 2 O emissions attributable to the increase in the model's N 2 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 2 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; 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 2 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 2 O emissions (Zhu et al., 2016;Saikawa et al., 2013;Tian et al., 2015;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 2 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 2 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).

Conclusions
Parts 1  and 2 of this study present three novel elements to gain insights into the functioning of the global carbon-nitrogen cycle: first, a record of the N 2 O isotopic history over the past 21 000 years, complementing the existing records over parts of the termination ; second, the first reconstructions of marine and terrestrial N 2 O emissions from the Last Glacial Maximum (LGM) F. Joos et al.: N 2 O changes from the Last Glacial Maximum to the preindustrial -Part 2 3531 to the preindustrial period as obtained by deconvolving the N 2 O and isotope ice-core records using a robust, established method; third, the application of a dynamic global vegetation model to simulate deglacial terrestrial N 2 O emissions and to attribute deglacial N 2 O emissions changes to different regions and governing processes. Our model results provide insight into the multi-decadal-to-millennial dynamics of the terrestrial C-N cycling by showing that the ice-core terrestrial N 2 O emission record could potentially be explained with a rapid adjustment of N cycling to the climate-and CO 2driven acceleration of the C cycle, but they do not exclude the possibility that alternative explanations linked to changes in N 2 O yield could be important.
The records provide the opportunity to evaluate models that simulate biological nitrogen fixation and nitrogen cycling (e.g. Meyerholt et al., 2016;Nevison et al., 2016;Thornton et al., 2007;Xu-Ri and Prentice, 2017;Wieder et al., 2015a) and feedbacks between climate change, land carbon, and terrestrial greenhouse gas emissions (Arneth et al., 2010;Stocker et al., 2013;Tian et al., 2018) for improved projections. Here, we evaluated LPX-Bern in transient simulations over the past 21 000 years. The model simulates an increase in N 2 O emission over the deglaciation but underestimates the magnitude of the deglacial increase in comparison with the ice-core data. This may add confidence in the sign of the climate-N 2 O feedback projected with this model for the next centuries , while suggesting that the magnitude of the projected feedback may be underestimated. Model results for regional changes and specific processes and forcings await confirmation by other studies.
The ice-core data presented in part 1 reveal a highly dynamic terrestrial nitrogen cycle. Remarkably, substantial changes in global terrestrial N 2 O emissions were realized within maximum 200 years, and possibly faster, in response to past Northern Hemisphere warming events. Reconstructed N 2 O emissions changed by up to 1 Tg N yr −1 and within less than 2 centuries at the onset of the NH warming events around 14.6 and 11.7 ka. These timescales are relevant for 21st century climate projections and much longer than accessible in typical laboratory or field experiments. Taken at face value, these results suggest that the nitrogen cycle and N 2 O emissions will also adjust on the global scale in the coming centuries.
The ice-core data and the N 2 O emission reconstruction in combination with models provide a new window to integrate across systems and relevant timescales to improve our understanding of the coupled C-N cycle in the Earth system. The results may help to put emerging results from observational field and laboratory studies into the context of decadal-to century-scale environmental and climate change. A few parameters of the nitrogen module of LPX-Bern are updated for this study compared to the previous version (v1.4) (Lienert and Joos, 2018a). These updated parameter values are set to obtain a preindustrial implied N source of 128 Tg N yr −1 and land biosphere N 2 O emissions of 5.9 Tg N yr −1 . The difference in parameter values and in key model outcomes between version v1.4N, used here, and version (v1.4) are summarized in Table A1. The parameter frac_soil_immob governing the fraction of the C and N flux from gross soil organic matter remineralization which is transferred back to the soil organic matter pool is adjusted to a higher value. This yields a smaller simulated BNF and smaller simulated N loss fluxes in v1.4N compared to v1.4. In turn, the parameters governing the yield in N 2 O emissions per unit N converted by nitrification (RN2ON), denitrification (RN2ODN), and leaching (RN2ONL) are re-calibrated in order to get preindustrial N 2 O emissions from the land biosphere of 5.9 Tg N yr −1 as in v1.4 and as suggested by observational estimates. In v1.4, the yield factor from nitrification is assumed constant. In v1.4N the same temperature dependency as for the yield for denitrification is assumed. This is a relatively minor modification as nitrification contributes only 12 % to the modelled preindustrial N 2 O emissions. We note that corresponding results and figures for version v1.4 are given in the paper published in Biogeosciences Discussion (https://www.biogeosciences-discuss.net/bg-2019-118/, last access: 30 June 2020).    Supplement. The supplement related to this article is available online at: https://doi.org/10.5194/bg-17-3511-2020-supplement.
Author contributions. FJ, RS, and BDS designed the study. RS, SL, and JM performed the simulations with LPX-Bern. SL, JM, RS, and FJ prepared the figures. HF and JS provided the ice-core data and BOB and ZL climate data. All the authors contributed to the writing and preparation of the manuscript.