N2O changes from the Last Glacial Maximum to the preindustrial-part II: terrestrial N2O emissions constrain carbon-nitrogen interactions

Carbon-nitrogen (C-N) interactions regulate N availability for plant growth and for emissions of nitrous oxide (N2O) as well as 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 N2O emissions 20 over the past 21,000 years as reconstructed from ice-core isotopic data and presented in part I of this study. Remarkably, the increase occurred in two steps, each realized over decades and within maximum two 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 N2O yield per unit N converted. We applied the LPX-Bern dynamic global vegetation model in deglacial simulations forced with Earth 25 System Model climate data to investigate N2O emission patterns, mechanisms, and C-N coupling. The N2O emission changes are mainly attributed to changes in temperature and precipitation and the loss of land due to sea level rise. The reconstructed increase in terrestrial emissions is broadly reproduced by the model, given the assumption that sources of reactive N positively respond to increasing N demand by plants. In contrast, assuming timeand demand-independent N sources in the model to mimic progressive N limitation of plant growth results in N2O emissions that are incompatible with the reconstruction. Our 30 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. Alternative mechanisms to explain the reconstructed N2O emissions include changes in N2O yield per N lost through gaseous pathways. [297 words]


Introduction
Nitrous oxide (N2O) 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. N2O is an important greenhouse gas and contributes to ongoing global 5 warming . It is also involved in the destruction of stratospheric ozone (Myhre et al., 2013). N2O 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 N2O increased from 270 ppb (MacFarling Meure et al., 2006) to around 330 ppb (https://www.esrl.noaa.gov/gmd/) over the industrial period due to human activities including fertilizer application, fossil fuel use and biomass burning Ciais et al., 2013). Atmospheric N2O varied naturally between around 180 10 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 N2O and climate change uncertain Battaglia and Joos, 2018;Kracher et al., 2016).
Variations in N2O emission, and thus in tropospheric N2O content, are closely linked to ecosystem processes governing the 15 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 N2O emissions are a sensitive indicator of the flow of reactive N entering and leaving land ecosystems (Firestone and Davidson, 20 1989) and reactive N needs to be available for the production of N2O as well as to support carbon uptake and storage by plants and the land biosphere.
The land biosphere sequesters about a quarter of anthropogenic CO2 emissions and is the largest natural N2O source (Ciais et al., 2013). Yet key features of the C-N cycle are poorly understood, leading to major uncertainties in global warming 25 projections (Joos et al., 2001;Friedlingstein et al., 2006;Friedlingstein et al., 2013;Zickfeld et al., 2013). The question to what extent N availability will limit future land C uptake and N2O 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 N2O emissions under climatic and environmental change (Niu et al., 2016). It is debated how global warming, increasing CO2, and increased N deposition affect the current land carbon sink 30 (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 observations (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 for 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 CO2 and climate Todd-Brown et al., 2013). Results from free-air CO2 enrichment (FACE) and open-top chamber experiments show that N limitation at least in some ecosystems reduces the response of aboveground biomass growth to 5 elevated CO2 (McMurtrie et al., 2008;Norby et al., 2010;Terrer et al., 2016) and that soil N and P availability 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 of reduced growth responses. Recent syntheses of the results from CO2 enrichment experiments suggests that mycorrhizal association exerts an important control on the magnitude of the realized 10 CO2 fertilization effect (Terrer et al., 2016) and that soil N and P availability are important controls on this magnitude in ectomycorrhizal and arbuscular mycorrhizal systems, respectively (Terrer et al., 2019).
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 time scales, relevant for global warming projections, is uncertain. For example , 15 Reich et al. (2018) found an unexpected reversal of C3 versus C4 grass response to elevated CO2 and shifts in soil N mineralization rates during a 20-year field experiment. These authors concluded that even the best-supported short-term drivers of plant response to global change might not predict long-term results. 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 20 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 CO2.
High-resolution data on the isotopic composition of N2O from Antarctic ice cores have the potential to provide precise information on past variations in terrestrial and marine N2O emissions and thus on C-N coupling on time scales from decades 25 to many centuries. A recent ice-core study on the stable isotope composition of N2O demonstrated the power of this approach (Schilt et al., 2014). The isotopic data showed that both land and oceanic sources increased during the interval from 16 to 10 thousand years before present (ka BP), when ocean circulation and climatic changes strongly affected the global cycling of both CO2 and N2O (Schilt et al., 2014). Schilt et al. concluded that natural N2O emissions will probably increase in response to global warming. In part I of this study (Fischer et al., 2019), this earlier work was extended to reconstruct the evolution in 30 terrestrial versus oceanic emissions of N2O from the Last Glacial Maximum (LGM; 21 ka BP) to the late preindustrial Holocene using novel N2O stable isotope data. The ice core data reveal large step-like changes in terrestrial emissions at the onset of warming events, realized over decades and, given the proxy data resolution, last 200 years at maximum. But a detailed process-based investigation of terrestrial N2O 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 II of this study is to improve our understanding of the cycles of N2O, C and N. The unique, new terrestrial N2O emission record of the past 21 kyr is used to explore mechanisms of the functioning of the C-N cycle on land and to quantify 5 terrestrial drivers for atmospheric N2O concentration changes. Controls on variations in terrestrial emissions are elucidated in the framework of a dynamic global vegetation model. Part I of this study (Fischer et al., 2019) presents the ice core N2O concentration and isotope data and the reconstruction of global terrestrial and marine N2O emissions for the past 21,000 years and provides a discussion on the reconstructed marine emissions in the context of past climate and oceanographic changes. A model-based interpretation of the reconstructed marine N2O emission changes during past abrupt climate events is given by 10 Joos et al. (2019). Here in part II we focus on the interpretation of the terrestrial N2O emission record and discuss terrestrial N2O emissions and C-N responses in transient deglacial simulations with a dynamical vegetation/terrestrial biogeochemistry model.
2 Introduction to N2O, terrestrial nitrogen and carbon flows 2.1 N2O budget and production mechanisms 15 Prather et al. (2015) estimated that the pre-industrial atmospheric lifetime of N2O was 123 years. Together with the atmospheric N2O concentration from ice cores, this figure constrains the total net pre-industrial N2O source to 10.5 ± 1 Tg N yr -1 . Marine N2O emissions were recently estimated by an observation-constrained approach, using water-column and surface N2O observations as targets, to be 4.6 (± 1 standard deviation range: 3.1 to 6.1) Tg N yr -1 . This calculation implies a natural terrestrial N2O source of 5.9 (4.1 to 7.7) Tg-N yr -1 (Battaglia and , in line with IPCC AR5 estimates of 6.6 (3.3 to 9.0) Tg N 20 N2O is produced through a variety of pathways both in the ocean and on land and N2O production is closely linked to the flows of C and N (Wrage et al., 2001;Chapuis-Lardy et al., 2006;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 25 dominant pathways of net N2O production are thought to be respiratory denitrification ( under low oxygen conditions, and autotrophic nitrification (NH3 →NO2 − → NO3 − ), mediated by archaea, bacteria and fungi. In addition heterotrophic nitrification, which is the oxidation of organic N to nitrite (NO2 -) and subsequent reduction to N2O by incomplete denitrification, was found to be the dominant path in a grassland ecosystem after fertilizer application (Moser et al., 2018). N2O is produced as a byproduct of nitrification and as an intermediate product during denitrification. N2O is also 30 produced from nitrite through nitrification, often termed nitrifier-denitrification, through anaerobic ammonium oxidation, chemautotrophic denitrification, abiotic processes or from photoautotrophic organisms in cryptogamic covers (Lenhart et al., 2015;Butterbach-Bahl et al., 2013).
Terrestrial N2O 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 on substrate and nutrient 5 availability and on 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 N2O in moist tropical areas, and lower emissions in high latitudes (Stehfest and Bouwman, 2006;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 N2O emission (Butterbach-Bahl et al., 2013;Zhuang et al., 2012), but the environmental dependencies of N2O emission are often complex 10 (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 N2O 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 CO2 emissions from soil organic matter decomposition, and a positive feedback of soil warming on N2O emissions is expected (Butterbach-Bahl et al., 2013). Warming 15 treatment increased measured N2O 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 N2O emissions (Dijkstra et al., 2012). Positive or neutral responses in N2O emissions have been found in field experiments under elevated CO2 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 20 fertilizer causes enhanced N2O 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 N2O produced per unit of N converted varies with environmental conditions and production pathway. For 25 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, 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 N2O production pathways and N loss processes and their relative importance may evolve through time and influence the N2O yield on local to global scales. 30

C-N-N2O coupling
The main flow paths of reactive N and its link to N2O production and C flows on land (Gruber and Galloway, 2008;Butterbach-Bahl et al., 2013;Vitousek et al., 2013;Zähle, 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 NH4 + and NO3by plant roots to support growth, while a large fraction of reactive N is taken 5 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 NO3pool is replenished by nitrification, the conversion of NH4 + to NO3 -. We note the acid-base equilibrium between NH4 + and NH3 in soil water; for simplicity, we generally refer to NH4 + only. 10 Turning to the external cycle, reactive N enters land ecosystems (blue arrows in Fig. 1) through the conversion of dinitrogen (N2) to organic N and eventually to NH4 + by BNF (Cleveland et al., 1999;Vitousek et al., 2013;Zähle, 2013;Sullivan et al., 2014;Xu 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 NHx and NOy (including sources by lightning). Reactive N is lost from land ecosystems 15 through gaseous losses (including N2O), leaching of NO3and dissolved organic matter by runoff, and emissions of N compounds by fire Ciais et al., 2013;Hedin et al., 1995;Hedin et al., 2003). On 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 to belong 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. 20 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 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 25 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 N2O) 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 depends 30 therefore critical on the flexibility of N input, hence BNF.
Empirical evidence from N-addition experiments, synthesized by Lue et al. (2011) andNiu et al. (2016), shows that the uptake of N by plants, net primary productivity and biomass as well as NHx and NO3pools in soils, nitrification, nitrate leaching, denitrification, and N2O 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 N2O. N2O production on land is predominantly associated with denitrification, and to a smaller extent with nitrification. Large gaseous losses and N2O 5 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 section 2.1 and 2.2, changes in terrestrial N2O emissions may be linked (i) to changes in the magnitude of 10 reactive N entering and leaving ecosystems, and (ii) to changes in the N2O yield per unit reactive N converted in land ecosystems.

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 15 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 Wania et al., 2009) ecosystems and their C and N dynamics Xu-Ri and Prentice, 2008;Xu-Ri et al., 2012), and describes the dynamic evolution of wetland and peatland extent . The model 20 calculates the release and uptake of the trace gases CO2, N2O Xu-Ri and Prentice, 2008;Xu-Ri et al., 2012) and CH4 (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 fifteen PFTs is used. Eight generic tree PFTs, and PFTs for C3 and C4-type grasses grow on natural 25 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 uses a 30 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 N2O fluxes and pools for plants and soils as schematically depicted in Fig. 1 and summarized below. 5 In LPX, 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. 10 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 -15 determine 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 (Tab. A.1) to match a total preindustrial reactive N source of 128 TgN 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. 20 Estimates of the global BNF for non-agricultural ecosystems are uncertain. They range from 40 to 470 TgN yr -1 , with most published estimates around 100 to 150 TgN yr -1 . Cleveland et al. (1999) used 100 plot-scale estimates of BNF to estimate global BNF on natural ecosystems to 195 TgN 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 TgN 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 TgN yr -1 , based on cost-benefit modelling for N fixation. 25 The same authors estimate asymbiotic N fixation to 22 TgN yr -1 by upscaling measurements reported in Cleveland et al. (1999).
In contrast, Elbert et al. (2012)  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 NH4 + and NO3pools 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 NH4 + and NO3 -) 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. 5 Nitrification is assumed to be proportional to the NH4 + soil pool with a temperature dependent rate coefficient. Denitrification is modelled as a two-step process whereby NO3is converted to nitrite (NO2 -) and NO2is further converted to N2, following Michaelis-Menten kinetics with dependence on the substrates NO3 -, NO2 -, as well as on labile C availability. Reaction rates are temperature dependent. NO3leaching depends on soil NO3 -, available water holding capacity, and daily runoff. Reactive 10 N is also lost from a grid cell by fire, assuming complete loss of N from burned vegetation, and NH3 volatilization. N2O emissions are computed by assuming that fractions of the N fluxes associated with denitrification, nitrification and N leaching are released as N2O. The calibration of the denitrification and nitrification yield factors is done to satisfy the constraints from the global atmospheric N2O budget, which suggests global terrestrial emissions to be 5.9 TgN yr -1 at pre-15 industrial (Battaglia and . Since global terrestrial N2O 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 N2O source path, this fraction is the product of a constant (RN2ODN) and a temperature-dependent factor f(T). f is unity at 22 o C and its value roughly doubles for a temperature increase of 10 o C. The amount of N2O released per unit N denitrified is therefore higher in 20 warm than in cold regions. The constant RN2ODN is set to 5.418 % of the denitrification rate. The corresponding yield fraction for N2O emissions from nitrification is nominally set to 0.231 % at 22 o C and the same temperature dependence (Smith, 1997) as for denitrification is assumed. The N2O yield for the leaching N flux is assumed constant and, for simplicity, set to the nominal value of nitrification (0.231 %).

25
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 N2O per N2 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, N2O yield, and preindustrial N2O emissions are not fully consistent. We carried out a sensitivity simulation 30 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 TgN yr -1 . This is about a factor of four higher than in the standard setup. Correspondingly, yield factors are about a factor four lower in this sensitivity run than in the standard simulation.
The sensitivity of simulated N2O 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 source is small (16 % versus 10 % in the standard) and the deglacial increase in N2O emissions is only 0.2 TgN yr -1 higher in the sensitivity than in the standard run (see Sect. 3.2 and 5   5), despite the large difference in the implied N source. Thus, related model-based conclusions for N2O 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 glacial-interglacial simulations is applied Ruosch et 10 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 CO2 (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)  The same model parameter values as determined by Lienert and Joos (2018b) are used, except for the immobilization fractions and N2O yield factors (Tab. A.1). Regarding the N module, the parameters are the same as in a previous studies Schilt et al., 2014) addressing N2O emissions over the past and under future global warming. An exception is an 25 adjustment in the upper limit of the fraction of NH4 + nitrified per day from 0.1 day -1 to 0.09096 day -1 at 20 o C and in the yield factors for denitrification, nitrification, and leaching in response to new observation-constrained estimates of marine N2O emissions and to a slightly revised estimate of the atmospheric N2O life time from 120 yr to 123 yr (Prather et al., 2015) as well as the simple representation of N uptake by N immobilization as discussed above.

30
Changes in N2O emissions and other model outcomes are attributed to individual driving factors (temperature, precipitation, CO2, 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 N2O emissions (eN2O) simulated in the baseline run and the sum of the emission changes attributed to individual drivers: eN2Oother-drivers= eN2Obaseline -eN2Otemperature -eN2Oprecipitation -eN2OCO2. The dominant driver is identified as having the largest contribution to eN2O 5 in the baseline run with the same sign of change. Grid cells that submerge under sea water or emerge from waning ice sheets during the period considered and grid cells with insignificant changes (|eN2O| < 1 mg N m -2 yr -1 ) are excluded from the spatially-resolved attribution.
The response time scales of LPX are investigated in a further "step-change" sensitivity simulation. Starting from the 10 equilibrated spin-up at 21 ka BP, climatic conditions and atmospheric CO2 are abruptly changed to conditions for 2.5 ka BP.
The run is continued for another 1500 years with climate and CO2 forcing for the period from 2.5 ka BP to 1 ka BP. 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 fast N2O emissions in the LPX model are able to respond to a sudden warming event, similar to the onset of 15 the B/A and the end of the Younger Dryas.

Reconstructed terrestrial N2O emissions and implications
We start the presentation of results by summarizing the main feature of the terrestrial N2O emission record (Fig. 3, green line) presented in part I of this study (Fischer et al., 2019). In part I, the global N2O emissions from land and from the ocean are jointly reconstructed by deconvolving novel ice core data of N2O and of its isotopic signature, 15N(N2O), using an established 20 method and relying on differences in the isotopic signature of land versus marine N2O emissions. Terrestrial emissions increased between LGM (21 kyr BP) and PI (1500 CE) by about 1.7 TgN yr -1 . Terrestrial emissions remained approximately invariant during the Heinrich Stadial I Northern Hemisphere (NH) cold phase (HS1; 17.4 to 14.6 ka BP, (Rasmussen et al., 2014)) until the start of the Bølling/Allerød NH warm period (B/A; 14.6 to 12.8 ka BP). Then, land emissions increased at the start of the B/A, declined again during the Younger Dryas NH cold period (YD, 12.8 to 11.7 ka BP) and peaked at the start of 25 the Preboreal period (PB), followed by a modest increase during the Holocene.
Remarkably, the overall deglacial increase in terrestrial N2O 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 N2O concentration and isotope data (see Fischer et al. (2019)) reveals that global terrestrial N2O emissions started to rise at the 30 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 N2O, its isotopic signature, and, correspondingly, in inferred emissions. Fischer et al. conclude that global terrestrial N2O emissions reacted within maximum 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 N2O emissions responded sensitively to climatic and environmental changes over the deglaciation. The rapid increase in terrestrial N2O emissions at the onset of the B/A and at the end of the YD 5 and the overall increase in emissions over the past 21,000 years either point (i) to an increase in N2O yield per unit N converted for emissions to the atmosphere, averaged globally and across all N2O production pathways, or/and (ii) to an increase in the global flux of converted N. Reactive N was available in sufficient amount to support an increase in global terrestrial N2O emissions during periods where environmental conditions became, on a global scale, more favorable for plant growth and C sequestration (Ciais et al., 2012;Bird et al., 1994;Joos et al., 2004;. 10 5 Transient simulations of terrestrial N2O emissions and the C-N cycle over the past 21,000 years

Simulated changes in global terrestrial N2O emissions and spatial patterns of change
We next explore governing mechanisms of the deglacial terrestrial N2O emissions and potential implications for the C-N cycles in the spatially resolved, mechanistic LPX-Bern model. LPX-Bern v1.4N simulates a general increase in global land N2O emissions over the deglacial period (Fig. 3). The simulated evolution matches the reconstructed change in terrestrial N2O 15 emissions from the ice core isotopic records relatively well, although modelled changes are smaller and typically less abrupt than reconstructed variations. The model represents the emission variations during the Bølling/Allerød (B/A) and Younger Dryas (YD) periods with peaks in emissions at the onset of the BA (14.6 ka BP) and the preboreal (11.7 ka BP) and an emission peak around 13.5 ka BP and corresponding minima at 14 ka BP and during the YD (12.8 to 11.7 ka BP). Reconstruction and models both show small changes in global terrestrial N2O emissions over the last 11 ka, the Holocene period. Simulated 20 terrestrial N2O emissions decreased somewhat between 9 and 8 ka BP, whereas reconstructed emissions slightly increased over the Holocene, leading to a discrepancy between simulated and reconstructed anomalies. Overall, the agreement between proxy reconstruction and model results supports the plausibility of the LPX-Bern model as well as of the underlying TraCE-21kyr climate input data used to force LPX-Bern on these long timescales. On the other hand, the reconstruction suggests relatively constant emissions from the land biosphere during Heinrich Stadial 1 (HS1) and a rapid rise in emissions at the start 25 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-21kyr 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 N2O emissions are the result of spatially differentiated responses. LPX-Bern simulates high 30 natural emissions of N2O 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 N2O 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 mgN m -2 yr -1 and the integrated flux between 20 o S and 20 o N amounts to 54 % of the global emissions, while emissions per unit area are low in 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 in mid-and low-5 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. N2O emissions decreased in a few regions, namely in Africa around 15 o 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 N2O emissions between 15 and 8 ka BP. The net loss of terrestrial N2O emissions by land area changes was dominated by the submergence of the high-productivity Sunda and Sahul Shelf regions. This loss 10 offsets about one-third of the G-IG increase in global terrestrial N2O emissions (Fig. 3). Changes in the land extent caused by changes in sea level represent an important factor for past global terrestrial N2O emissions.
The patterns of change in terrestrial N2O 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 N2O emissions during the 15 Holocene, LPX-Bern simulated large regional shifts in source strength, linked to changes in temperature and precipitation.
This includes for example a decrease in N2O emissions from 11 to 0.5 ka BP 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 20
In this section, we address C-N coupling in LPX-Bern and analyze 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 N2O emission change in the model: (i) changes in the source of reactive N, fueling nitrification and denitrification, and (ii) changes in the N2O yield per unit N converted. 25 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 gN m -2 yr -1 in northern mid and high latitudes and around 1 gN 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 gN m -2 yr -1 is simulated in many semi-arid regions, including southern Africa, the sub-Sahara region, India, northern Australia, and in the southern parts of North America. Soil mineral N (Fig. 5B) is typically 30 below 0.5 gN m -2 in tropical forests and around 1 gN m -2 in temperate and boreal forests and tundra regions, while higher values of up to and more than 10 gN m -2 are simulated in exatratropical 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 N2O yield factors, i.e., the N2O 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 5 all other settings as in the standard. The deglacial warming leads to a higher N2O yield in the standard compared to this sensitivity run and 0.44 TgN yr -1 of the deglacial increase in land N2O 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 N2O emissions as driven by the increase in the flow of reactive N in LPX-Bern.

10
Turning to changes in the sources of reactive N, the simulated N source increased by 12 % from 113 TgN yr -1 at the LGM to 128 TgN yr -1 at 0.5 ka BP in the standard setup. Biological nitrogen fixation (BNF) is assumed to adjust dynamically to an increase in N demand 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 time scales and that nutrient input into ecosystems by deposition, and weathering (plus BNF for N) is large enough to support plant growth. These 15 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 N2O emissions allows us to critically evaluate this assumption in a quantitative way on the multi-decadal to centennial time scales 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 20 keeping these rates of N source constant in a sensitivity simulation over the past 21 kyr. This more strongly N-limited simulation yields not an increase, as reconstructed, but a decrease in global terrestrial N2O emissions (Fig. 3, black line). The decrease in global land N2O emissions is related to the loss of land in response to sea level rise (-0.11 TgN yr -1 , and a small decrease in N2O emissions on remaining land (-0.088 TgN yr -1 ) in this sensitivity simulation.

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

30
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 constant N source per land use class (Fig. 6A,B). In the standard, the N source decreased over the deglaciation by 7 TgN 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 TgN yr -1 in the northern and southern extratropics, 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 5 denitrification (Fig. 6C). A decrease in mineral N is simulated in the standard in regions where also C stocks decreased or changed little including parts of mid-latitude Eurasia, Africa, Australia and of the southern US. In contrast soil reactive N concentrations decreased in many extratropical regions in the simulation with constant N source (Fig. 6D). Globally integrated, soil mineral N decreased by 19 % over the deglaciation in the run with constant N source, compared to a negligible change in the standard setup. 10 These differences in BNF (including abiotic input) and the N cycle between the run with constant N source and the standard case have profound impacts, not only for the emission of N2O, 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 constant N source. LPX-Bern predicts that most of the NPP increase is realized in the northern extratropics. The NPP increase 15 in the extratropics is almost twice as large in the standard than in the sensitivity run. As a result of the generally higher NPP , the deglacial change in C storage is in the standard about 402 PgC larger in the extratropics and about 127 PgC larger in the tropics than in the run with constant N source.
In summary, the simulation with constant N source completely fails to reproduce the reconstructed N2O emissions from the 20 land biosphere. The increase in N2O 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 N2O emissions in this sensitivity simulation. We note, however, that changes in yield due to processes not incorporated in LPX-Bern could potentially explain the reconstructed increase in N2O 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, terrestrial N2O emissions increase as reconstructed. 25

Attribution of simulated terrestrial N2O emissions to climatic and environmental drivers
The complex spatio-temporal changes in land N2O emissions are attributed to climatic and environmental drivers (Fig. 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 TgN yr -1 to the overall emission increase, followed by changes in precipitation (0.7 TgN yr -1 , Fig. 7). Changes due to increasing atmospheric CO2 and orbitally driven changes in 30 photosynthetically active radiation slightly offset the emission increase, while non-linear interactions among the different drivers contribute to the global increase. The change in land extent due to ice sheet melting and sea level rise leads to a reduction in terrestrial N2O emissions (0.5 TgN 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 N2O 5 emissions, which occurred earlier in the simulation than the ice core reconstruction, is not driven by CO2 fertilization and an associated change in N ecosystem flows, but by the prescribed climate change from TraCE-21kyr. 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 BA in the standard simulation is due to deficiency in the response of the LPX-Bern to early deglacial climate change, or to deficiencies in the climate input data, or a combination of the two. 10 Individual drivers exert regionally distinct influences and these may vary between different periods. Here, we attribute the spatial changes to changes in temperature, precipitation, CO2 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 of the TraCE-21kyr input data (Fig. 9). Generally, attributed changes in emissions find their counterpart in underlying changes in individual 15

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 20 Amazon region and in parts of Africa and India caused emissions to fall in this period. The pattern of change in N2O 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 N2O 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. 25 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 30 precipitation with generally increasing emissions under increasing precipitation. Increasing CO2 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 o S to 20 o 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 CO2 in this region (Fig. 8E,F) most likely reflect non-linear interactions between precipitation, temperature and CO2 . Figure 10 further illustrates which driver exerts the largest influence on simulated regional emission changes in a given grid 5 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, 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 CO2 generally play a secondary role. 10

The time scales of response to a step change in climate and CO2
Next, we investigate potential causes of the relatively slow N2O increase during past abrupt NH warming events in the standard simulation compared to the ice core reconstruction. To this end, we assess on which time scales N2O 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 N2O emissions as reconstructed from the ice core data for the abrupt warming events at the onset of the B/A and 15 the end of the YD or whether the intrinsic response time scales of the model are too long to match reconstructed emissions.
We construct an extreme case for rapid warming events by switching climate (and CO2) 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 20 response in global N2O 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 N2O of more than 100 years into account, such a fast multi-decadal increase in N2O emissions would cause a century scale increase in atmospheric N2O concentrations, similar to what is seen in the ice core record. The adjustments in global nitrification and denitrification fluxes evolve similar as for N2O emissions with the main increase in these fluxes again realized within about 40 years (not shown). 25 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 NO3availability. Modelled global NPP ( Figure 11C) increases immediately after the step and about three quarters of the NPP change are realized in the first year 30 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 NO3leading to the decrease in N leakage, while warming accelerates soil decomposition and the release of NH4 + . The newly assimilated C and N is allocated to leaves, sapwood and hardwood and roots, before released to litter and labil soil organic matter overturning on time scales from years to decades. These annual-to-decadal vegetation and litter turnover time scales govern the initial response time of N source, nitrification, denitrification and N2O emissions. The century-scale response is linked to the equilibration time scales 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 5 which is linked to remaining peat in the step-simulation formed under LGM conditions, but absent in the reference simulation.
In conclusion, this sensitivity experiment demonstrates that N2O emissions in the LPX-Bern model are indeed able to adjust within decades to less than a century to abrupt warming events, given prescribed forcing is changing fast enough. We conclude that the exact timing of simulated N2O emissions at the onset of the B/A and the end of the YD depends sensitively on the 10 climate forcing data prescribed in the LPX-Bern model.

Reconstructed terrestrial N2O emissions and implications for the C-N cycle
Terrestrial N2O 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 15 hemisphere warming and to shifts in the Intertropical Convergence Zone (ITCZ) and precipitation patterns. Each step was realized within maximum 200 years. and thus on a time scale similar to that of the ongoing anthropogenic climate perturbation.
It has remained somewhat unclear whether warming will increase or reduce regional to global scale N2O emissions, as responses to warming treatment are found to be highly variable across a range of conditions and ecosystems (Dijkstra et al., 20 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 N2O emissions globally. This suggests, as noted earlier (Schilt et al., 2014) and as projected , that terrestrial N2O 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).

25
The increasing terrestrial N2O emissions imply either an increase in N2O yield per unit N converted, averaged globally and across all N2O production pathways, or an increase in the global flux of converted N or a combination of these factors.
Regarding yield, uncertainties in the interpretation of the terrestrial N2O emission record are linked to uncertainties in the ratio between N2O produced to N converted. This yield factor is known to vary with environmental conditions and across N2O production pathways, though quantitative experimental evidence is often not unequivocal (Diem et al., 2017;Davidson et al., 30 2000;Firestone and Davidson, 1989;Saggar et al., 2013;Butterbach-Bahl et al., 2002). For nitrification, the ratio of N2O/NO3produced increases with increasing acidity and decreasing oxygen (Firestone and Davidson, 1989) and increasing temperature (Smith, 1997). For denitrification, the production ratio of N2O/N2 yield increases with increasing NO3availability, increasing oxygen concentration, decreasing decomposable carbon, decreasing pH and decreasing temperature. Changes in precipitation may have altered the N2O yield over the deglaciation. Higher soil water content and associated anoxic conditions generally favor the conversion of N2O to N2 and results in a low yield of N2O (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 5 and yield (Bouwman et al., 2002). Yield is low under low NO3availability (Diem et al., 2017) and generally found to increase when NO3becomes more available (Weier et al., 1993). This relation, when considered in isolation, implies increasing N2O emissions to be indicative of increasing NO3availability. The ratio of N2O to N2 emitted generally decreases with increasing temperature (Firestone and Davidson, 1989). On the other hand, 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 10 Davidson, 1989). Further, experimental studies on grasslands yield both higher and lower ratio of N2 to N2O under elevated CO2 compared to ambient CO2 (Zhong et al., 2018;Baggs et al., 2003). In addition, different N2O 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 the position to evaluate whether the N2O yield factor for the combined N conversion processes increased or decreased over the deglaciation on the global scale. 15 Changes in globally-averaged yield may possibly explain the reconstructed terrestrial N2O 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 and 20 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 of years to accumulate (Michalski et al., 2004) . Microorganisms require readily decomposable carbon substrate for the denitrification of nitrate, otherwise nitrate may accumulate under 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 25 of the YD when climate warmed, precipitation pattern 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 N2O emissions. Such a scenario appears in conflict with the reconstruction, given that reconstructed N2O emissions remained high during the early B/A and even slightly increased throughout the Holocene. 30 Terrestrial N2O 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. Variation in reactive N availability among diverse land classes are found to correspond to variations in N2O emissions (Davidson et al., 2000). In addition, cryptogamic covers fix large amounts of N and contribute substantially to global N2O emissions (Elbert et al., 2012). The ice core data show that reactive N must have been available in sufficient quantity to support increasing terrestrial N2O 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 CO2 over the deglacial period and contributed to meet the N demand of plants, nitrifiers, denitrifiers and cryptogamic covers under more favorable growth conditions. At the same time, 5 the ice core reconstructions do not suggest that reactive N for the production of N2O remained scarce and that marine emissions dominated past atmospheric N2O variations. The magnitude of possible adjustments in the N source is unclear. Wang and Houlton (2009) modelled an increase in N fixation under increasing CO2 and temperature. Yet 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 10 for BNF of around 25 o C, an assumption that has been challenged by  who found the abundance of N-fixing trees to increase monotonically with temperature by analyzing 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  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 N2O emissions. 15 The ice-core data and emission estimates represent global values with century-scale resolution. They do not permit us to discriminate regional variations, nor 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 20 deglaciation. Indeed, it has been hypothesized that an increase in BNF may become cost-effective under more severe N limitation in 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 25 conditions . Regional changes may have contributed to the deglacial N2O emission increase.

Simulating deglacial terrestrial N2O emission changes with LPX-Bern
We applied the LPX-Bern dynamic global vegetation model to investigate regional N2O 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 30 and plants in a simplified way. C and N is 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 time scales on each grid cell. LPX-Bern, being a spatially explicit model, reacts differently in different regions to changes in climatic and environmental drivers (Fig. 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 5 organic matter according to their accessibility to microbial action (Averill and Waring, 2018). Potential N2O release by lichens and bryophytes is also not modelled (Porada et al., 2017). Only net N2O production during denitrification is considered, while gross production and consumption of N2O during denitrification (Schmid et al., 2001b) are not explicitly modelled. In future work, N2O may be treated as an intermediate product of denitrification to explicitly simulate N2O consumption as well as the related ratio of N2O to N2 emissions, i.e., the yield factor. 10 Uncertainties in modelled N2O 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 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. Yet, a fraction of 6 % of NPP is directly transferred to a pool with a short overturning 15 time to represent root exudates. The by far largest contribution 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 N2O 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, 2010a). Most frequently, empirical 20 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 extratropical regions (Gerber et al., 2010) or on the optimization of plant C investment into resource acquisition (Fisher et al., 2010). Meyerholt et al. (2016) describes asymbiotic BNF as a function of temperature, shading and soil moisture. 25 The implementation of such different parameterizations in LPX-Bern would likely lead to different estimates for deglacial changes in BNF and N2O 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 N2O 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 30 result is expected for BNF depending on evapotranspiration, because globally-averaged evapotranspiration changes little over the deglaciation in our standard simulation (Tab. A.1) 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 (Tab. A.1). Overall, one might expect a similar or even larger increase in BNF and N2O 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. 5 The N loss in LPX is predominantly driven by local gaseous loss from denitrification, with a much smaller role for fire, leaching, and minor contributions from volatilization of NH3, 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 NH4 + available for aerobic nitrification and of NO3available for anaerobic denitrification within a grid 10 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 CO2.
LPX-Bern was applied in an earlier study to simulate climate-N2O feedbacks under global warming (Stocker et al., 2013) 15 and results for global emissions for the period from 16 to 10 ka BP are presented in Schilt et al. (2014) LPX-Bern forced by TraCE-21kyr output represents the reconstructed emissions reasonably well, while differences between reconstructed and simulated global land N2O emission remain. The modelled deglacial increase in global N2O emissions is at the lower bound of the reconstructed range. 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 CO2 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 N2O 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 5 N2O 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 N2O emissions. Regarding the time scales of response, results of a sensitivity simulation demonstrate that LPX-Bern is able to represent abrupt, multi-decadal scale change in global terrestrial N2O emissions as reconstructed for past rapid warming events. Further, factorial simulations demonstrate a complex interplay of different driving factors for deglacial N2O emissions with considerable synergistic or antagonistic 10 interactions related to changes in temperature, precipitation, and CO2. Remarkably, about a third of the simulated global increase in terrestrial N2O emissions over the past 21,000 was counteracted by land loss due to flooding of previously exposed shelf regions.
We applied LPX-Bern to analyse the consequences of alternative C-N cycle hypotheses. In the model's standard setup, the N 15 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 favorable growth conditions and resulted in a higher flow of N through land ecosystems and increased N2O production from terrestrial ecosystems as well as increased carbon storage on land. The increase in N2O emissions, linked to this increase in N source, is further amplified by an increase in the yield of N2O per N converted by nitrification and denitrification. 20 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 N2O emissions, in clear conflict with the ice core data. The increase in N2O emissions attributable to the increase in the model's N2O 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 25 the decrease in N2O emissions in this run with 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 30 increasing N2O 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 N2O emissions (Zhu et al., 2016;Saikawa et al., 2013;Tian et al., 2015;Huang and Gerber, 2015;Zaehle and Friend, 2010b;Olin et al., 2015;Goll et al., 2017) reveal similar or alternative mechanisms and whether these models are able to represent the reconstructed terrestrial N2O emissions over the past 21,000 years. 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
Part I (Fischer et al., 2019) and II of this study present three novel elements to gain insights into the functioning of the global 5 carbon-nitrogen cycle: First, a record of the N2O isotopic history over the past 21,000 years, complementing the existing records over parts of the termination (Schilt et al., 2014); second, the first reconstructions of marine and terrestrial N2O emissions from the Last Glacial Maximum (LGM) to the preindustrial period as obtained by deconvolving the N2O and isotope ice core records using a robust, established method; third, the application of a dynamic global vegetation model to simulate deglacial terrestrial N2O emissions and to attribute deglacial N2O emissions changes to different regions and governing 10 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 N2O emission record could potentially be explained with a rapid adjustment of N cycling to the climate and CO2-driven acceleration of the C cycle, but they do not exclude the possibility that alternative explanations linked to changes in N2O yield could be important.

15
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 and Prentice, 2017;Wieder et al., 2015a) and feedbacks between climate change, land carbon, and terrestrial greenhouse gas emissions (Arneth et al., 2010;Tian et al., 2018) for improved projections. Here, we evaluated LPX-Bern and found that the model is able to simulate the evolution in global terrestrial N2O emissions over the past 21,000 years in reasonable agreement with the ice core emission data. This 20 adds confidence to projections of the climate-N2O feedbacks with this model . Model results for regional changes and specific processes and forcings await confirmation by other studies.
The ice core data presented in part I reveal a highly dynamic terrestrial nitrogen cycle. Remarkably, substantial changes in global terrestrial N2O emissions were realized within maximum 200 years, and possibly faster, in response to past northern 25 Hemisphere warming events. Reconstructed N2O emissions changed by up to 1 TgN yr -1 and within less than two centuries at the onset of the NH warming events around 14.6 and 11.7 ka BP. These time scales are relevant for 21 st 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 N2O emissions will also adjust on the global scale in the coming centuries The ice-core data and the N2O emission reconstruction in combination with models provide a new window to integrate across 30 systems and relevant time scales 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.

5
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 TgN yr -1 and land biosphere N2O emissions of 5.9 TgN yr -1 . The difference in parameter values and in key model outcomes between the version v1. 4N, used here, and version (v1.4) are summarized in Tab. A.1. 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 10 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 N2O emissions per unit N converted by nitrification (RN2ON), denitrification (RN2ODN), and leaching (RN2ONL) are re-calibrated in order to get preindustrial N2O emissions from the land biosphere of 5.9 TgN 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. We note that 15 corresponding results and figures for version v1.4 are given in the manuscript published in Biogeosciences Discussion (https://www.biogeosciences-discuss.net/bg-2019-118/). Transpiration and interception (mm yr -1 ) -2.5 (-1 %) -5.8 (-2 %)        Bengtsson, G., Bengtson, P., and Mansson, K. F.: Gross nitrogen mineralization-, immobilization-, and nitrification rates as a function of soil C/N ratio and microbial activity, Soil Biology and Biochemistry, 35, 143-154, 2003. 30