Warming and redistribution of N inputs drive long-term increase in terrestrial N2O emission factor

38 Anthropogenic nitrogen inputs cause major negative environmental impacts, including 39 emissions of the important greenhouse gas N 2 O. Despite their importance, changes in 40 terrestrial N loss pathways driven by global change and spatial redistribution of N inputs 41 are highly uncertain. We present a novel coupled soil-atmosphere isotope model (IsoTONE) 42 to quantify terrestrial N losses and N 2 O emission factors from 1850-2020, initialised using 43 a global dataset of natural soil δ 15 N, and optimized with a tropospheric timeseries of N 2 O 44 isotopic composition using a Bayesian framework. 45 N inputs from atmospheric deposition caused the majority (51%) of anthropogenic 46 N 2 O emissions from soils in 2020. Long-term growth in emissions was driven by fertilization 47 and deposition, however biological ﬁxation caused subdecadal variability in emissions. N 2 O 48 emission factors (EF) show large spatial variability due to climate and soil parameters. The 49 mean eﬀective global EF for N 2 O (weighted by N inputs) was 4.3 ± 0.3% in 2020, much higher 50 than the land surface area-weighted mean (1.1 ± 0.1%). Climate change and redistribution 51 of fertilisation have driven an increase in global EF over the past century, which accounts 52 for 18% of the anthropogenic soil ﬂux in 2020. Predicted increases in fertilisation in 53 emerging economies will accelerate N 2 O-driven climate warming in coming decades, unless 54 targeted mitigation measures focussing on fertiliser management and reduced N deposition 55 are introduced. 56


Introduction
Nitrous oxide (N 2 O) is a long-lived greenhouse gas and a key stratospheric ozonedepleting substance (Ravishankara, Daniel, & Portmann, 2009;Tian et al., 2020).The atmospheric N 2 O mole fraction has increased from ∼270 nmol mol −1 in the preindustrial era to >332 nmol mol −1 today (Tian et al., 2020;WMO, 2020).The primary global source of N 2 O is production during N cycling by microbiota in soils, which also releases NO and N 2 .N cycling and thus N gas production are strongly enhanced by direct and indirect anthropogenic N inputs (Park et al., 2012;Tian et al., 2020Tian et al., , 2019;;World Bank, 2019).
Agricultural fertilisation accounts for around two thirds of anthropogenic N inputs (Hurtt et al., 2020), with the remainder contributed by biological N 2 fixation and deposition of NO x and NH 3 .Most anthropogenic N is not incorporated into crops or soils but lost to the environment (Sebilo, Mayer, Nicolardot, Pinay, & Mariotti, 2013), representing huge monetary losses for the agricultural sector and causing a cascade of environmental problems (Fowler et al., 2013;Galloway et al., 2008;Roy, Hammond Wagner, & Niles, 2021).In the coming decades, N inputs are expected to grow in line with increasing agricultural production, and to shift towards tropical regions and emerging economies as strict N pollution controls are introduced in many developed regions (Leitner et al., 2020;Lett & Michelsen, 2014;Peng et al., 2020;Verma & Sagar, 2020;W. Wang et al., 2019).
Effective management of N fertiliser to achieve high N use efficiency is key to balancing food production with environmental protection and thus reducing N 2 O emissions (Venterea et al., 2012;Wagner-Riddle, Baggs, Clough, Fuchs, & Petersen, 2020).However, mitigation is challenging, and N 2 O emissions are currently exceeding the highest predicted scenarios (Tian et al., 2020).
N is lost from terrestrial ecosystems through several major pathways: Microbial and abiotic N gas production in soils; leaching of N species; and ammonia volatilization.Despite their importance, expected changes in N losses in the coming century are not well known (Galloway et al., 2008;Gruber & Galloway, 2008).Nitrification (aerobic) and denitrification (anaerobic) are the main processes emitting N gases (NO, N 2 O and N 2 ) from soils.On the global scale, the impacts of climate change on N-gas production processes are poorly known: Warming is generally expected to enhance microbial activity, however interactions between factors such as N availability, plant growth, and precipitation changes are poorly constrained.Moreover, it is unknown if increased nitrification or denitrification rates would manuscript submitted to Nature Geoscience in fact lead to increased N gas production (Hartmann, Barnard, Marhan, & Niklaus, 2013;Inatomi, Hajima, & Ito, 2019;Li et al., 2020).The proportion of N lost by leaching globally is not expected to change significantly over time in response to warming due to the contrasting effects of increased N mineralization and reduced moisture availability (He et al., 2018;Mas-Pla & Menció, 2019;Reichenau, Klar, & Schneider, 2016;Tian et al., 2020).However, leaching losses and predicted responses to warming vary widely between different regions depending on soil, hydrological and ecosystem parameters, and may also be strongly affected by precipitation regime changes (Frank et al., 2015;Stuart, Gooddy, Bloomfield, & Williams, 2011).Moreoever, increasing atmospheric CO 2 generally enhances plant growth and N uptake, thus impacting the availability of N in soils for different loss pathways (Eisenhauer, Cesarz, Koller, Worm, & Reich, 2012;Mitchell, Mitchell, Driscoll, Franklin, & Lawlor, 1993).Process models have been used to simulate N loss pathways in a changing climate (e.g.Giltrap and Ausseil (2016); Xu-Ri and Prentice (2008); Xu-Ri, Prentice, Spahni, and Niu (2012)).These models used generally require large amounts of input data and parameterisations as well as high computing power, which makes it difficult to constrain model parameters with observations using inversion frameworks, and complicates investigations of global or long-term emissions.Top-down modelling efforts can give robust estimates of global emissions for recent decades (R. L. Thompson, Ishijima, et al., 2014;R. L. Thompson et al., 2019;R. L. Thompson, Patra, et al., 2014), however without the incorporation of isotopes, these approaches cannot provide mechanistic information.
The isotopic composition of soil N (δ 15 N soil ) has been proposed as an integrated indicator of N loss partitioning in natural systems (Bai, Houlton, & Wang, 2012;Craine, Brookshire, et al., 2015;Craine, Elmore, et al., 2015;Houlton & Bai, 2009).Leaching of soluble N species (eg.NO − 3 ) involves very low isotopic fractionation, while losses through ammonia volatilization and N gas production strongly favour 14 N and thus cause 15 N enrichment in the remaining soil N pool (Bai et al., 2012).Observations have shown that mean global δ 15 N soil is elevated relative to N inputs, reflecting significant production of gaseous N species (Bai et al., 2012;Houlton & Bai, 2009), although the relationship may be unpredictable in some regions due to N immobilization (Craine, Elmore, et al., 2015).Previous studies have used δ 15 N soil models to constrain N losses at individual sites or globally for natural ecosystems (Bai et al., 2012;Houlton & Bai, 2009), however this approach has not been applied to estimate temporal changes in N loss pathways.Furthermore, the availability of complementary datasets to validate previous δ 15 N soil models has been limited to short-term measurements from individual sites.Recent advances in spectroscopic isotope instrumentation have delivered high precision long-term time series of background tropospheric N 2 O mixing ratio and isotopic composition (Harris et al., 2017;L. Yu, Harris, Henne, et al., 2020), allowing an integrated view of N 2 O sources and sinks.These results have been used to estimate total anthropogenic N 2 O emissions based on two-box models of the atmosphere (Sowers, Rodebaugh, Yoshida, & Toyoda, 2002;Toyoda et al., 2013;L. Yu, Harris, Henne, et al., 2020), however this approach cannot give spatially resolved information on N 2 O sources.
We aim to gain new insight into the N loss processes underlying global N 2 O emissions and their spatiotemporal patterns, in order to understand how N losses are changing under the influence of anthropogenic activities and climate change.We use an artificial neural network to estimate a global isoscape of natural soil δ 15 N.This is used to initialise a soil module to simulate spatially resolved terrestrial N losses via leaching, volatilization and gas production pathways.N 2 O emissions from the soil module are released to a twobox atmospheric module to simulate N 2 O mixing ratio and isotopic composition from the preindustrial era to the present day.The novel coupled model framework, 'IsoTONE', is optimized within a Bayesian framework using a high precision time series of N 2 O mixing ratio and isotopic composition from several background sites, as well as estimates of N 2 O emission factors from the Global N 2 O Database.

Results and discussion
2.1 Terrestrial N 2 O emissions 12 key parameters in the IsoTONE model were optimized using a Markov Chain Monte Carlo (MCMC) approach with 120 000 iterations, described in detail in Section A.3 and Table A.1.Total terrestrial N 2 O emissions from the optimised model (Table 1) agree well with previous results, providing confidence in the novel isotopic basis of the model.Total terrestrial N 2 O emissions (Figure A.12) for 1860 and 2010 were 5.3±0.4 and 12.6±0.7 Tg N 2 O-N a −1 respectively, showing good agreement with 1860 and 2010 estimates of 6.3±1.1 and 10±2.2Tg N 2 O-N a −1 from the N 2 O Model Intercomparison Project (Tian et al., 2019).
The range of 12.6 to 13.9 Tg N 2 O-N a −1 modelled for 2007-2016 additionally agrees with a recent global metaanalysis estimating total terrestrial emissions of 12.2-23.4Tg N 2 O-N a −1 for the same period (Tian et al., 2020), and with the mean of 12.9 Tg N 2 O-N a −1 from a metaanalysis by Scheer, Fuchs, Pelster, and Butterbach-Bahl (2020).Total global soil NO emissions for 2010 were estimated to be 13.7±3.9Tg NO-N a −1 , agreeing well with the high end of estimates from a metaanalysis suggesting global soil NO emissions of 1.8-12.3Tg NO-N a −1 (Pilegaard, 2013).
The spatial distribution of emissions from fertilisation, deposition and fixation for natural and anthropogenic soils is shown in Figure  2021)) using data from 123 ground-based sites (Figure 1), with significant differences only seen in small isolated regions.The IsoTONE framework assumes that land use changes have had a minor impact on N 2 O emission factors compared to the pre-existing variability in emission factors driven by climate and soil parameters -this assumption is supported by the good agreement between inversion and IsoTONE emission estimates.The largest differences are seen in tropical South America, Africa and Australia, which may be due to the scarcity of atmospheric monitoring stations available to constrain inversion estimates in these regions (R. Thompson, 2021).Differences may also relate to extensive cland use change and ultivation of N-fixing soybean crops in South America.Furthermore, fundamental differences in N cycling and immobilization in tropical and Arctic regions may affect the accuracy of the IsoTONE model in these regions, as most experimental and field studies were conducted in temperate soils (Craine, Brookshire, et al., 2015;Voigt et al., 2020).Total N 2 O emissions agree very well between the two models, with 13.3±0.1 Tg N 2 O-N a −1 and 13.9±0.9Tg N 2 O-N a −1 predicted for 2019 by the CAMS inversion and IsoTONE respectively.Table 1.Changing characteristics of the total and anthropogenic terrestrial N2O flux at the beginning of the anthropocene (1850) and through the past century.
The natural flux and the area-weighted EF do not vary temporally.Growth rate is the 10-year mean growth rate centred on the year of interest.The area-weighted EF is the mean global EF calculated using the areas of grid cells as weights; the N-input weighted EF is calculated using the total N inputs of grid cells as weights.The incorporation of isotopic composition means that, unlike previous global models, the IsoTONE framework distinguishes between different N 2 O production pathways (Section A.5).The model results furthermore suggest that on average, laboratory-measured fractionation factors for these pathways are only expressed at ∼55% on average in soils (frac ex = 0.55±0.05,Table A.1).This is a global average which we expect to vary widely between individual sites depending on soil structure, moisture, and the depth of N 2 O production and consumption.For example, soils with small N cycling and gas diffusion rates would be expected to have higher isotopic fractionation during N loss processes compared to soils total nitrification rate, which is not estimated in IsoTONE, so the results cannot be directly compared.
Globally, just 30% of fertiliser N appears to be available for N cycling (fert EF red, Table A.1) -the remaining fertiliser N is primarily incorporated into harvest, and may also be immobilised via increased soil storage, or lost through leaching pathways which are not explicitly modelled.This is consistent with previous metaanalyses, which show that 15-70% of fertiliser N is taken up by plants, a significant proportion could remain in soils, and 5-25% is unaccounted for in plant or soil pools and thus lost to the pathways modelled in this study (Corbeels, Hofman, & Van Cleemput, 1999;Gardner & Drinkwater, 2009;Jenkinson, Poulton, Johnston, & Powlson, 2004;Sebilo et al., 2013).Therefore, emission factors for fertiliser emissions can be significantly lower than for other N inputs, with implications for the applicability of EF measurements from agricultural sites.Mean global fertiliser N incorporation could be increasing as plant growth and thus N use is enhanced in response to increasing atmospheric CO 2 , however water and nutrient limitations will also play a role in regulating plant growth (Tian et al., 2020) -this potential effect is not currently captured in IsoTONE and should be a focus of future model studies.

2.2
The anthropogenic N 2 O budget: Inputs, losses and trends The anthropogenic N 2 O soil flux in 2020 was estimated to be 7.1±0.6Tg N 2 O-N a −1 , close to the highest projected emission scenario (RCP8.5)estimate for 2020 (IPCC, 2014;Tian et al., 2020) (Figures 2 and A.11). 1.7±0.4,3.6±0.3and 1.8±0.6Tg N 2 O-N a −1 of the anthropogenic flux were contributed by soil emissions from fertilisation, deposition and fixation N inputs respectively (Table 1), with an additional 1.7 Tg N 2 O-N a −1 from non-soil anthropogenic emission sources (emissions from EDGAR for categories 1A1, 1A3b, 2B and 6, Section 4.2.2).This agrees very well with a recent ensemble analysis, which estimated a total anthropogenic flux of 7.3 (4.2-11.4)Tg N 2 O-N a −1 (Tian et al., 2020).We find that deposition N accounts for 41±14% of all anthropogenic emissions in 2020 compared to 19±12% direct emissions from fertilisation, 21±15% from enhanced fixation, and 19% from non-soil sources: Deposition N inputs clearly contribute the majority of anthropogenic emissions.This finding is in contrast to the results of Tian et al. (2020), who report that direct N 2 O emissions from fertilisation are dominant.These contrasting results are due in part to the classification of all emissions above the 1850 baseline (eg.enhanced N 2 O from natural sites due to warming, fixation and deposition) being classified as anthropogenic in this study, whereas in previous studies some or all of these processes are not considered in the calculation of the anthropogenic burden.Moreover, our results suggest that emission factors could be significantly underestimated from field measurements (see Sections 2.1 and 2.3) -particularly at agricultural sites where a significant proportion of N is removed through harvest or immobilization as described by fert EF red -which may lead to some underestimation of deposition emissions in different model frameworks.The 10-year mean growth of total anthropogenic N 2 O emissions was between 0 and 0.04  et al., 2019;Tian et al., 2020).This caused the rate of change for the tropospheric N 2 O burden to peak at 0.026 nmol mol −1 a −2 (Figure 2).After 2015, the growth rate of emissions strongly decreased, meaning the 10-year mean growth for 2010-2020 was 0.14 Tg N 2 O-N a −1 a −1 and thus within then normal range for the last half-century (Figure 2).Fluctuations in the growth rate for tropospheric background N 2 O mixing ratio reflected the growth rate manuscript submitted to Nature Geoscience in emissions driven by different input categories.Variability in fixation N inputs dominates subdecadal interannual variability in total terrestrial N 2 O emissions, such as the 2010-2015 peak, contributing >70% of variability within decadal bins (Figure 2).Fixation inputs are 2-10 times more variable than other input types at subdecadal timescales, thus accounting for their key role in driving variability in N 2 O emissions.In contrast, changes in fertilisation inputs drive the changing growth rate of anthropogenic emissions at timescales larger than 1-2 decades.Increases in both fertilisation and deposition are responsible for the strong and constant increase in N 2 O emissions and tropospheric background mixing ratio over the last century.These findings could be impacted by uncertainties in the input datasets, which should be investigated further using targeted isotopic and modelling approaches at the site and regional scale in regions with particularly high uncertainty.
The isotopic composition of the total anthropogenic source reflects changing emission processes (Figure A.11). δ 15 N of the anthropogenic source is higher than the mean soil source (Table 1) and shows a strong decreasing trend, which may indicate an increasing dominance of agricultural emissions (Prokopiou et al., 2017).During the growth rate acceleration of 1945-1980, SP of the anthropogenic source was relatively constant, however after 1980 it strongly increased.These fluctuations were also observed by Prokopiou et al. (2017), who used a two-box model to interpret N 2 O source isotopic composition from firn air data, and found similar but more uncertain values for anthropogenic source isotopic composition.
Variability in SP and δ 15 N are in opposite directions, and thus unlikely to be caused by

Drivers of spatiotemporal variability in N 2 O emission factor
Soil moisture, N content, mean annual precipitation and soil bulk density (21, 9, 5, 4% of variability respectively) were the main parameters causing broad geoclimatic gradients in N gas production and N 2 O EF between different regions (Figures 3 and A.13), consistent with previous laboratory and field results (Cantarel et al., 2012;Smith, Thomson, Clayton, Mctaggart, & Conen, 1998;Yao et al., 2010;Zhang et al., 2019).Mean annual temperature was important for total N gas production, but not for EF (5 vs. 0.2% of variability in f gas and EF).3).This builds upon the observation of Tian et al. (2020), who showed that the mean global EF for N 2 O from agricultural soils is significantly higher than the IPCC default value.Measurements of mean annual EFs using chambers and similar methodologies are far more time intensive and expensive than measurements of soil δ 15 N, thus isotopic modelling of the N cycle can be used to understand regional variability in   & Köchy, 2012).Aridity index at a spatial resolution of 30 arc seconds was taken from the Global Aridity Index and Potential Evapotranspiration Climate Database (v2) (Trabucco & Zomer, 2019).Global soil pH at >60 000 sites worldwide was taken from the database compiled by Slessarev et al. (2016).The Worldwide Organic Soil Carbon and Nitrogen Dataset (Zinke, Stangenberger, Post, Emanual, & Olson, 2004) was used to estimate soil organic nitrogen content with data from >4 000 sites rasterised using linear interpolation with scipy.Total fertiliser N inputs were taken from the Land Use Harmonization Database (LUH) with the historical dataset used for 1800-2015 (LUH2 v2h) and the future forcing dataset for 2015-2020 (LUH2 v2f) (Hurtt et al., 2020) and deposition and fixation N inputs (annual) (Houlton, Wang, Vitousek, & Field, 2008;Kowalczyk et al., 2006;Peng et al., 2020).All datasets were converted to the model grid with 0.5 • ×0.5 • resolution for -180 • to +180 Data for δ 15 N soil for >6 000 soil samples from natural (non-agricultural) sites was compiled by Craine, Elmore, et al. (2015).This database also includes ancillary data such as MAT, MAP, carbon and nitrogen concentration, C:N ratio, soil texture, density and pH, and site coordinates and elevation, although not all parameters are available for all samples.
Geographical coverage was improved by adding a further 748 samples, including 112 samples from Australia (BASE database; Bissett et al. (2016)) and 392 from Africa -regions which were underrepresented in the Craine database.The majority of data was from near-surface soil (0-50 cm depth).Within this depth range, no significant effect of depth on δ 15 N soil was seen in the data.The sample-specific reported ancillary data was compared to the global gridded datasets (Section 4.1.1)at the same geographical locations to gapfill the ancillary data (Section A.1). Relationships between ancillary data and δ 15 N soil measurements were then used to predict δ 15 N soil using an artificial neural network (ANN) with the Python package Keras (Chollet (2015); details in Section A.1).Using the ANN, global δ 15 N soil on an 0.5 • ×0.5 • grid was estimated in order to drive the soil isotope module (Section 4.2.1).
A bootstrapping approach was used to estimate uncertainty in gridded δ 15 N soil .

Atmospheric data
An tropospheric background time series of N 2 O mixing ratio and isotopic composition since the preindustrial era formed the core set of observational data used to constrain simulation results in this study (Figure 5).The background mixing ratio, for the purposes of this study, is defined as the concentration of a given species when the impact of local or recent sources and sinks is absent -also known as the baseline concentration (IUPAC, 2019).The primary dataset comprised tropospheric background measurements of N 2 O mixing ratio, δ 15 N and N-isotope SP from the Jungfraujoch High Alpine research station and the Cape Grim Air Archive (CGAA), as well as firn air sampled during the 2018 East GReenland Ice coring Project (EGRIP) campaign (L.Yu, Harris, Henne, et al., 2020;L. Yu et al., 2022) covering the period from 1982 to the present day, all measured at Empa, Switzerland.Additionally, mixing ratio, δ 15 N and SP data from CGAA (1979CGAA ( -2005) ) and firn air  from Park et al. (2012) were used.These were corrected to the Empa dataset scale using the mean offset in the overlap period, with offsets of -1.1❻ and -1.0❻ found for firn and CGAA data for δ 15 N, and +1.5❻ and +1.3❻ for firn and CGAA data for SP.To extend the isotopic timeseries further back, ice core and interstitial snowpack air measurements from the Greenland Ice Sheet Project II (GISP II) (Sowers et al., 2002) were also used, covering the period from 1785-1990.This dataset showed no offset for δ 15 N compared to the Empa dataset, and did not include SP measurements.Other available datasets (e.g.

Parameters optimised in MCMC
Observations Model results

Legend:
Non-soil anthropogenic and oceanic N 2 O emissions

Soil nitrogen module
The soil N module of IsoTONE estimates fractionational losses of N to different pathways based on δ 15 N soil (Section 4.1.2).The module is based on the equations used by Houlton and Bai (2009) and Bai et al. (2012), which state that the δ 15 N soil reflects the balance between N inputs (i), with a relatively constant rate and isotopic composition, and variable N losses.Loss pathways are divided into leaching (L), with low isotopic fractionation; NH 3 volatilization, with high isotopic fractionation but a low contribution to total losses (N H 3 ); and abiotic and microbial gas production (G), with high and variable isotopic fractionation.At steady state, inputs must equal outputs for both total amounts and isotopic composition: manuscript submitted to Nature Geoscience where ε is the fractionation factor ( k15 N/k14 N expressed in permil, where k is the reaction rate) for the respective loss pathway and δ 15 N i was set as -1.5❻ (Bai et al., 2012).Assuming steady state, the fraction of N loss accounted for by leaching (f L ), NH 3 volatilization (f N H3 ) and gas production (f G ) can be calculated from δ 15 N soil .The steady state assumption is only valid for natural ecosystems, thus only δ 15 N soil measurements from natural sites were used to initialise the model.Unlike Bai et al. (2012), 14 N and 15 N were traced separately in Eq. 2, whereby the rate of 15 N loss via a pathway is the rate of 14 N multiplied by the fractionation factor, for example, for leaching: The rate of 15 N addition was calculated analogously.Tracing isotopes separately facilitated calculation of the stepwise impact of N losses on the remaining N reservoir isotopic composition ('Rayleigh fractionation' (Mariotti et al., 1981)), and necessitated an iterative solution to Eq.
2. Rayleigh fractionation describes how the isotopic composition of reactant and product pools change as a reaction proceeds depending on the fractionation factor and the fraction of reactant consumed, assuming both reactant and product pools are well mixed (Mariotti et al., 1981).
We assumed that N loss processes are primarily linked with the inorganic N pools (i.e.NO − 3 and NH + 4 ), and estimated fractionation factors related to the two species.
Fractionation during N mineralization was assumed to be minor (Möbius, 2013).Soil organic N was not considered, given that the rates of turnover processes and their associated isotopic fractionation effects are very scarcely studied at the global scale.Fractionation during leaching (ε L ) was estimated to be 1❻ (Bai et al., 2012).Fractionation during NH 3 volatilization (ε N H3 ) was estimated to be -17.9❻based on the equation in Stern, Baisden, and Amundson (1999).Fractionation during N gas production (ε G ) was estimated based on the relative contributions of the two dominant microbial processes: • Fractionation for N 2 O production during nitrification of NH + 4 (ε nit ) was estimated as −56.6 ± 7.3❻ with SP of 29.9 ± 2.9❻ (Denk et al., 2017).
In have recognised that soil gas is not a homogeneous, well-mixed environment -instead, gas production and consumption occurs in pores that are only partially connected relative to the rate of microbial processes (Bai & Houlton, 2009;Wen et al., 2016).This means that Rayleigh fractionation processes occur within all pores and microenvironments, and thus the effective net fractionation may be significantly lower than calculated from measured fractionation factors, often referred to as 'underexpression' (Bai & Houlton, 2009).We therefore introduce a parameter called the fractionation expression factor (frac ex) to the model such that net or effective fractionation = ε × frac ex for all modelled processes.
Thereby, frac ex = 1 would reflect a well-mixed soil gas reservoir where effective fractionation is equal to expectations from laboratory measurements, and frac ex = 0 would reflect a completely closed soil gas environment where all reactions go to completion within separate pores and thus no effective fractionation is observed.The initial value of frac ex was set to 0.7 and it was optimized in the MCMC framework (Section 4.4) assuming a uniform error distribution between 0.3 and 1.0.In reality, different values for frac ex may be expected for different processes and environments, in particular depending on soil texture and pore structure, however this cannot be determined within the available model-data framework.
The soil model was run on an 0.5 • × 0.5 • grid from -60 • to 80 To solve Eq. 2 for each grid cell, the soil module was run over four iterations (n = 4).For the first iteration, f G was set at 0.1 for all gridcells.f N H3 was estimated using the CABLE model (Kowalczyk et al., 2006;Peng et al., 2020).f N H3 is nearly always <0.05, therefore the model results are not highly sensitive to the parameterisation of f N H3 .f L was then determined as 1 − f G − f N H3 .Partitioning of total gas losses (f G ) into N 2 O, N 2 and NO (f N2O , f NO , f N2 ) was estimated with a sigmoid fit of WFPS to available experimental measurements (Section A.2, Figure A.5), with mean global WFPS from CABLE used to constrain partitioning for each grid cell in the model.Similarly, denitrification (f denit ) and nitrification (f nit ) contributions to N 2 O production were estimated with a sigmoid fit to available experimental data (Section A.2, Figure A.6). Mean WFPS will not describe all variability in gas partitioning, given the highly variable and non-linear nature of N gas emissions, however it provides the best available estimate based on the current status of experimental and modelling research.An overall ε G was estimated based on f denit and f nit and the fractionation factors for the individual processes multiplied by frac ex, with a mean value of -30❻.This is lower than in Bai et al. (2012) where a range of -16 to -20❻ was assumed for ε G without explicit consideration of microbial pathways -measurements made since the publication of Bai et al. (2012) have shown that fractionation is much larger than previously estimated (Denk et al., 2017).To achieve steady state, 10 cycles of N addition and removal were conducted within each iteration.First, N inputs for 14 N and 15 N were added to the initial soil N pool.N was then removed via leaching, volatilization, and gas production for 14 N and 15 N.The steady state soil pool δ 15 N value at the end of each iteration (δ 15 N ss ) was compared to the δ 15 N value for the gridcell (δ 15 N soil ) to find an improved f G for the next iteration according to: Following four iterations, δ 15 N soil and modelled δ 15 N ss agreed within 0.01❻ and final values of f G and f L were accepted.
In addition to the calculation of f G and f L , the soil module estimated N 2 O production rate and isotopic composition for each grid cell.The production rates for 14 N-N 2 O and 15 N-N 2 O were calculated from the fractionation factors for each pathway multiplied by frac ex; both 15 N α and 15 N β were traced separately to consider N 2 O SP.For example, the 14 N reaction rate for nitrification N 2 O production was estimated as f G × f N2O × f nit and the corresponding reaction rate for 15 N as: Analogously, rates for each isotopic variant were also calculated for denitrification production and reduction.The production of N 2 was assumed to represent the consumption of N 2 O via complete denitrification, as no other significant N 2 sources are known, thus the extent of N 2 O reduction could be estimated and its isotopic effect calculated.N 2 O production and consumption were estimated in each step of the steady state calculation to give a final N 2 O flux and isotopic composition (δ 15 N and SP) for each gridcell to pass on to the atmospheric module described in the following subsection.

Atmospheric N 2 O module
The atmospheric module takes the fractional losses, estimated by the soil module described in the previous subsection, and uses a two-box model representing a well-mixed troposphere and stratosphere to reconstruct N 2 O background mixing ratio and isotopic composition.Absolute N emissions into the tropospheric box were calculated using the fractional estimates of losses for each grid cell from the soil module, combined with annual N inputs to each grid cell for fixation, deposition and fertilisation for the period 1800-2020.Fertiliser N inputs were taken from the LUH database and covered the whole time period (see Section 4.1.1).Deposition N inputs for 1860-2050 were from Lamarque et al. (2013).Fixation (for 1901-2100) was estimated using the CABLE model with scenario A1 as described by Peng et al. (2020), where fixation is calculated using resource optimization with temperature dependence (Houlton et al., 2008).The 1901 fixation values were used to estimate fixation in earlier years and the 1860 values for deposition in earlier years, thus assuming negligible anthropogenic influence before 1850.
Deposition and fixation inputs were assumed to generally fulfil the steady state criteria required by the soil module, however harvest N exports mean that this assumption is not valid for fertiliser N inputs.Harvested N is often returned to the soil or atmosphere via manure or deposition at a different location (Schlesinger, 2009).Previous studies have shown that a large proportion of fertiliser N is incorporated into crops and therefore removed during harvest, as well as potentially stored in soils (Corbeels et al. (1999) 2020) also report that recent increases in EF for direct soil emissions are likely due to climate change feedbacks and interactions, rather than a direct increase in EF due to increasing fertiliser application.The majority of plant N uptake is not relevant in the IsoTONE model, as it is returned to the soil system at the same rate it is removed within the steady state assumption.This is not the case for harvested agricultural plant matter -the isotopic impact of harvest N exports cannot be accounted for within the scope of this study, as the required input data (in particular the proportion of N removed in harvest per grid cell annually) is not currently available.Harvest N will be explicitly incorporated into a future version of the IsoTONE framework.
Previous studies have shown that microbial N gas production is temperature sensitive (Hu, Chen, & He, 2015;Inatomi et al., 2019;Pilegaard, 2013) and likely increasing in a warming climate by between 0.5 and 1.0 Tg N 2 O-N a −1 • C −1 , with a best estimate of manuscript submitted to Nature Geoscience Tian et al., 2019;Xu-Ri et al., 2012;Zaehle, 2013).We assume that this increase in emissions relates to a general sensitivity of microbial activity to temperature, and is therefore likely to affect production of N 2 and NO as well as N 2 O.We therefore incorporated a temperature sensitivity increase of 10% of microbial N emissions from the year 1800 per degree of warming (temp sens = 1.1±0.04)for N 2 O, N 2 and NO, which corresponds to 0.6 Tg N 2 O-N a −1 • C −1 according to the best estimate of 6.3±1.1 Tg-N a −1 for N 2 O emissions in the preindustrial era (Tian et al., 2019).N losses through leaching were consequently reduced to maintain N mass balance.
Yearly soil emissions of N 2 O, NO and N 2 were calculated for each input type for each grid cell incorporating both fert EF red and temperature sensitivity of emissions.
In addition, EDGAR (Emission Database for Global Atmospheric Research, Crippa et al. for wastewater treatment (Harris et al., 2017).Isotopic composition of natural soil N 2 O emissions was given by the soil module for each grid cell, and weighted by the emissions per grid cell to calculate overall isotopic composition of soil N 2 O emissions.The δ 15 N isotopic composition of fertiliser N inputs is estimated to be 3❻ (Bateman & Kelly, 2007;Savard et al., 2010) compared to -1.5❻ for natural N inputs (Bai et al., 2012), thus increasing the δ 15 N of emitted N 2 O.Total global emissions were found by adding the oceanic emissions (F ocean ) to the global terrestrial emissions using the flux and isotopic compositions listed in Table A.1.
Total global N 2 O was emitted into a two-box atmospheric model comprising a wellmixed troposphere and a well-mixed stratosphere, based on the model described in L. Yu, Harris, Henne, et al. (2020).Emissions were assumed to be immediately well-mixed through the atmosphere, as the time required for tropospheric mixing is estimated to be around one year (Bowman & Cohen, 1997), which is the time resolution of the model.Previous versions of this model explicitly calculate the required preindustrial terrestrial N 2 O flux to achieve steady state in the preindustrial era accounting for best estimates of tropospheric N 2 O mixing ratio, N 2 O lifetime, and stratosphere-troposphere exchange.However, in the atmosphere module of IsoTONE, preindustrial terrestrial emissions (F terr ) are prescribed by the soil module.An iterative calculation (maximum of 20 iterations) is therefore used to optimize the preindustrial oceanic N 2 O flux (F ocean ) and the troposphere-stratosphere exchange term (T to S) by shifting these two terms stepwise according to their prior uncertainty (Table A.1) using two equations until steady state is achieved.First, tropospheric N 2 O inputs are balanced against losses to the stratosphere: where MR is the N 2 O mixing ratio.The preindustrial stratospheric N 2 O mixing ratio at steady state is then calculated based on stratosphere-troposphere exchange and stratospheric N 2 O destruction (Toyoda et al., 2013;L. Yu, Harris, Lewicka-Szczebak, & Mohn, 2020): where MW N 2O−N is the molecular weight of N in N 2 O (28 g mol −1 ) and moles trop and moles strat are the moles of the troposphere and stratosphere (1.5 and 0.27×10 20 moles) respectively.This iteration meant that optimized values for both F ocean and T to S were found for each solution of the model, although they were not explicitly targeted in the inversion described in the following subsection.F ocean was not varied with time, as recent results suggest that the oceanic flux is relatively stable compared to the terrestrial flux (Tian et al., 2020).Smaller, dynamic aquatic and semiaquatic ecosystems, such as estuaries and manuscript submitted to Nature Geoscience coastal wetlands, show highly dynamic and variably N cycling and N 2 O emissions (Brase, Bange, Lendt, Sanders, & Dähnke, 2017;Moseman-Valtierra et al., 2011;Wells et al., 2018), but are too small to be resolved at the global scale of this study.
Once steady state was achieved for preindustrial N 2 O fluxes and mixing ratio, model calculations proceeded as described by L. Yu, Harris, Henne, et al. (2020) and will therefore only be briefly presented here.First, net stratospheric fractionation and the isotopic composition of the preindustrial stratosphere were found for δ 15 N and SP assuming steady state.The model was then run forwards with annual time steps using annual terrestrial emissions and isotopic composition provided by the soil module of IsoTONE.At each time step, N 2 O inputs and destruction were considered to estimate the rate of change in N 2 O mixing ratio and isotopic composition, and thus calculate a time series of mixing ratio, δ 15 N and SP for a well-mixed troposphere (L.Yu, Harris, Henne, et al., 2020).

Optimization of model N cycle parameters
Twelve key model parameters were optimised using a Markov Chain Monte Carlo (MCMC) approach (Table A.1), with several datasets used to constrain the results.The geoclimatic gradients (spatial variability) were constrained using N 2 O emission factors from the Global N 2 O Database, binned for 16 climatic zones (Section 4.1.4).Temporal variability was constrained using a combined background tropospheric timeseries of N 2 O mixing ratio and isotopic composition (δ 15 N and SP) from several different sites, averaged for 25-year blocks from 1740 to 1940 and for 2 year blocks from 1940 to 2020 to give a total of 49 data points (Section 4.1.3and Figure 5).Data was reduced to 16 climate zone EFs and 49 temporal data points to avoid strong overweighting of recent atmospheric results, which are much more frequent and less uncertain than older measurements (e.g.biweekly monitoring at Jungfraujoch station (L.Yu, Harris, Henne, et al., 2020)) but highly covariable.Observation uncertainty was defined as the standard deviation of measurements within each spatial or temporal block.Model uncertainty was set at 0.5 nmol mol −1 for N 2 O mixing ratio, 0.1❻ for δ 15 N and SP, and 0.5 for climate zone EF.The incorporation of isotopic composition in the atmospheric model gave an implicit sensitivity to the spatial distribution of emissions, as the isotopic composition of emitted N 2 O depends on the dominant emission processes in a particular gridcell, thus allowing the model to distinguish between changes in different regions and input types.
The MCMC was run with three different stepsizes: 0.75, 0.5 and 0.25.Within each iteration i of the MCMC, Gaussian parameters P G (see Table A.1) were varied according to: where r unif is a uniformly distributed random number between -1 and 1. Uniform parameters P U were varied according to: where R max and R min are the maximum and minimum of the parameter uncertainty range.
Independent values of r unif were determined for every parameter.Observation uncertainty follows a Gaussian distribution, thus observations were also varied within each iteration using Eq. 8, however r unif was determined separately only for different groups of observations, e.g.N 2 O mixing ratio, N 2 O SP.
Once parameters were varied within an iteration i, the Metropolis rule was applied to determine if parameters and observations could be accepted (Rayner, 2020).If both were accepted, the model was run for the parameter set i, and the model-observation probability was calculated for i.The Metropolis rule was then applied to determine if model-observation i could be accepted compared to i−1. 5 000 iterations were run at each step size sequentially until a total of 120 000 iterations had been run (40 000 at each step size), which was sufficient  A.1).The parameters F ocean and T to S (see Section 4.2.2) were not explicitly varied in the MCMC but were calculated to achieve steady state in the atmospheric module in each accepted iteration, thus posterior estimates for these parameters were also obtained.

Estimating uncertainty in posterior model results
The uncertainty in posterior parameters show in Table 1 was estimated as the standard deviation of all accepted results.To estimate uncertainty in the posterior model results, 100 iterations of the model were run using 100 randomly selected sets of accepted parameters.
The standard deviation of results from all 100 iterations was used to estimate the uncertainty in the final model results.Standard error propagation was used to estimate uncertainty in all subsequently calculated values, eg.ratios and sums across time or space.
5 Data and code availability statement δ 15 N soil data from >7000 samples from natural (non-agricultural) sites was compiled as described in Section 4.1.2,however ancillary data for these measurements (e.g.MAT, MAP, pH) was incomplete.Ancillary data reported for the point measurements was therefore compared to the global gridded datasets (Section 4.1.1).Mean annual temperature (R 2 =0.96, slope=0.98),precipitation (R 2 =0.87, slope=0.67)and pH (R 2 =0.60, slope=0.62) agreed well between the point and gridded datasets.Soil carbon showed moderate agreement (R 2 =0.3, slope=0.35),while soil N and C:N ratio showed no correlation between the point and gridded datasets.The relatively poor for soil C and N can be attributed to the high spatial variability of these parameters; additionally, the gridded datasets refer to organic C and N while the point data does not always report whether organic or total C and N were measured.A complete ancillary dataset was compiled by gapfilling missing point values of MAP, MAT and pH with gridded dataset values.For C, N and C:N only gridded values were used, as point and gridded datasets were inconsistent and point values are likely less comparable between sites.Other relevant parameters, such as bulk density, which were not reported or very sparsely reported in the point measurements were taken from the gridded datasets.
An overview of the relationships between ancillary data and δ 15 N soil was gained using single linear correlation analyses and a multivariate principal components analysis (Figure A.1).The strongest positive linear correlations were found between δ 15 N soil and MAT, N and bulk density, and the strongest negative correlations between δ 15 N soil and C, C:N, elevation and aridity index.Only parameters with <1000 missing values were used for the principal components analysis, therefore elevation and soil texture were not included.PC1 was dominated by the broad geoclimatic gradient between wetter and drier areas (based on MAP), which did not relate strongly to δ 15 N soil , in contrast to previous expectations (Amundson et al., 2003).PC2 (24.6% of variability) showed a relationship between MAT, N, C, C:N and δ 15 N soil evident from the single linear correlations, thus reflecting the impact manuscript submitted to Nature Geoscience The δ 15 N soil data is not globally representative and many large regions have very sparse data coverage, thus interpolation of the data to a global grid is not possible.An artificial neural network (ANN) was therefore used to predict δ 15 N soil based on geoclimate parameters and thus estimate a global grid of δ 15 N soil values using the 'Keras' package (Chollet, 2015) in Python.A sequential ANN was compiled with: 1. 9 input nodes; 2. A dense layer with 16 nodes using the rectified linear unit (ReLU) activation function; 3. A dense layer with 16 nodes using a hyperbolic tangent (tanh) activation function, and L2 regularization with a regularization parameter of 0.01; 4. Identical to layer 3. but using a linear activation function; 5. A batch normalization layer; 6.A dense output layer with one node using a ReLU activation function.
Weights were initialised with Xavier initialization using a uniform distribution (default initialization in Keras).The Adam optimizer was used with a learning rate of 0. A bootstrapping approach was used to find the ANN model best able to predict the observation data and to estimate the uncertainty in predicted δ 15 N soil .First, the relevant uncertainty in the δ 15 N soil measurements was determined; this relates primarily to the degree to which a single measurement is representative of soils in a grid cell, and is much higher than the measurement uncertainty (usually <0.2❻ for IRMS measurements).The The final ANN was able to reproduce the observations with a slope of 1.04, R 2 of 0.41 and wRMSE of 2.6❻, thus exhibiting significantly more explanatory power than a multiple linear regression, which was only able to reproduce data with R 2 of 0.27 and RMSE of 3.1❻.
The role of specific parameters in the ANN was examined using two methods: i) using the ANN to predict δ 15 N soil with all inputs set to their mean values, except for the parameter -25-manuscript submitted to Nature Geoscience of interest, to understand how much of the variability in δ 15 N soil is contributed by the parameter, and ii) shuffling all data for the parameter of interest before predicting δ 15 N soil , to measure the reduction in goodness-of-fit.MAT was the most important parameter in the ANN, accounting for 30% of variability in δ 15 N and an RMSE reduction of 0.8❻, followed by absolute longitude (14.5%, 0.3❻), latitude (9.5%, 0.2❻), pH (9.8%, 0.07❻) and bulk density (8.7%, 0.16❻).The importance of MAT as well as soil composition agrees with previous studies by Amundson et al. (2003) and Craine, Elmore, et al. (2015).N 2 O EFs are predicted to be highest at low pH, consistent with the observations in the from the 'Global N 2 O Database' (Hu et al., 2015;Inatomi et al., 2019).Both water availability and bulk density impact O 2 availability, which is known to be key in determining nitrification and denitrification rates and thus N gas emissions (Hu et al., 2015;Pilegaard,

2013
).The observations are inconclusive regarding the impact of soil N availability on EF.The correlation analysis shows a relationship between fertilisation and EF only for the Pearson correlation coefficient, which is potentially more affected by outliers, but no relationship between EF and soil N availability.Soil N and fertilisation do not plot strongly against PC1 or PC2; fertilisation shows a strong positive relationship and soil N availability shows a weak negative relationship to EF along PC3.This is consistent with previous studies, which show that although more N availability clearly leads to more N gas production, the relationship between N availability and N gas EFs is unclear (Butterbach-Bahl et al., 2013;Cui et al., 2021;Inatomi et al., 2019;R. L. Thompson et al., 2019).Tian et al. (2020) report increasing EF for direct soil emissions over recent decades that is attributed to climate change feedbacks and interactions as well as spatiotemporal variability of climate, rather than a direct increase in EF due to increasing fertiliser application.In order to compare modelled EFs to the EF data, which came from a range of agricultural ecosystem types with varying management practices, measurement methods and manuscript submitted to Nature Geoscience Experimental data to constrain the relative contributions of nitrification and denitrification to N 2 O emissions is sparse; only 45 measurements were available, from grassland (Harris et al., 2021;Ibraim et al., 2019;Wolf et al., 2015) and rice paddy (Verhoeven et al., 2019) systems.Although the data was not representative of global ecosystem types, a wide range of WFPS was covered.A sigmoid curve fit the data well (Figure A.6): Unexpectedly, denitrification makes a slightly lower contribution to N 2 O emissions at higher WFPS, because a significant amount of N 2 O is reduced to N 2 before emission at high WFPS.However, the scarcity of data means that this parameterisation is highly uncertain, and further laboratory and field measurements of the partitioning between pathways are needed to better constrain the roles of nitrification and denitrification.Furthermore, the denitrification to nitrification ratio is likely also dependent on agricultural management and fertiliser application, estimation of which is currently beyond the scope of this model.
Recent developments in spectroscopic isotope measurements offer a promising method for online, in situ quantification of nitrification and denitrification N 2 O emissions (Harris et al., 2021;L. Yu, Harris, Lewicka-Szczebak, & Mohn, 2020), which offers the potential for improvements to the parameterisation as new data from a wider range of sites becomes available.
manuscript submitted to Nature Geoscience A.3, A.5 and A.6), N2O denit (dotted line) and N2Onit (dashed line) as a proportion of total N gas emissions using the posterior fit described in Eqs.A.3, A.5 and A.6.

A.3 Model optimization
120 000 iterations of the MCMC framework were run and 1 667 solutions were found.The MCMC approach led to a strong improvement in model-observation agreement (Section A.3), with RMSE decreasing from 14.1 to 1.7 nmol mol −1 for N 2 O mixing ratio, 1.5 to 0.5❻ for δ 15 N, 1.5 to 1.1❻ for SP, and 5.4 to 1.8% for climate-zone emission factor (Figure A.9).The improvement in SP is lower than for δ 15 N because measurements are more uncertain relative to the trend, and the timeseries is much shorter.Prior observation EFs were highly uncertain, and therefore EFs were the only observation parameter showing a significant difference between the prior and the posterior.Posterior observation EFs were higher than prior observations EFs, which is likely due to the extremely dynamic nature of  worse, with modelled EFs were higher than observations.Emissions are highly variable and dynamic in the wetland and tropical regions (Arias-Navarro et al., 2017;Kort et al., 2011;Rowlings, Grace, Scheer, & Liu, 2015), thus it is very challenging for observational studies to accurately capture annual emissions, which may account for this mismatch in these regions.Total N losses to nitrification were estimated as the sum of NO emissions and nitrification N 2 O emissions per grid cell, and denitrification total N losses were similarly estimated as the sum of N 2 emissions and denitrification N 2 O emissions (Bai et al., 2012) (Figure A.10).Summing across all grid cells therefore allowed calculation of total nitrification and denitrification N losses over time, as well as estimation of the proportion of total N losses from each pathway, and the contribution of each pathway to N 2 O emissions.The proportion of total N inputs lost to both pathways began decreasing around 1945, when fertilisation began to increase, as a significant amount of fertiliser N is removed and relocated via harvest.
Since 1990, nitrification and denitrification loss fractions have slightly increased, as N inputs move from temperate to tropical regions which favour N gas production (Figure A.10).
Total N losses to nitrification were somewhat lower than to denitrification, totalling 25.2±1.8 and 21.2±1.7 Tg-N a −1 in 2020 respectively.These are lower than suggested by a recent studies (Scheer et al., 2020;Schlesinger, 2009)  The anthropogenic N 2 O flux was calculated by assuming that all emission increases since 1850 are anthropogenic (Figure A.11).Most of the input datasets did not include data before 1850, so it was not possible to calculate a mean pre-anthropogenic baseline for either mixing ratio or isotopic composition over several decades prior to 1850, or to account for possible anthropogenic influences before 1850.Although natural emissions would show some variability due to natural climate fluctuations, natural emission variability is expected to be much lower than the key anthropogenic influences of climate change and agriculture (Tian et al., 2020).This method of estimating the total anthropogenic N 2 O flux includes both direct and indirect emissions, as well as increased emissions from natural ecosystems due to anthropogenic climate change feedbacks, and is thus not directly comparable with many other studies that report only the flux directly accounted for by anthropogenic N inputs.
The isotopic composition of the anthropogenic source was found using: where δ represents either δ 15 N or SP of N 2 O, 'total' designates the total source isotopic composition (directly modelled as the sum of fluxes across all gridcells), subscripts 'a' and 'n' refer to anthropogenic and natural sources respectively, and f a is the fraction of N 2 O from anthropogenic sources.Anthropogenic source isotopic composition was found for years where anthropogenic emissions contribute at least 10% of the total flux.

Figure 1 .
Figure 1.Modelled total terrestrial N2O fluxes from the IsoTONE model (a) and from CAMS (b) (R. Thompson (2021); regridded to 0.5 • using linear interpolation) for the year 2019.Both panels use the same logarithmic colour scale.The difference between the IsoTONE and inversion estimates is shown in c; areas where the difference is significant compared to the uncertainty are highlighted with black outlines.
with high diffusivity.In 2020, 60% of N 2 O emitted from soils (7.5±0.4Tg N 2 O-N a −1 ) was produced from denitrification and 39% (4.7±0.8Tg N 2 O-N a −1 ) from nitrification.The contribution of nitrification to soil N 2 O emissions has decreased very slightly from 40% in 1850 to 39% in 2020 (Figure A.10).The spatial distribution of nitrification-N 2 O from IsoTONE agrees very well with the global map produced by Pan, Lam, Wang, Mosier, and Chen (2021), with high nitrification in areas such as the Amazon basin, sub-Saharan Africa, Europe and coastal Australia.However the model of Pan et al. (2021), centers around

Figure 2 .
Figure 2. Temporal evolution of N inputs and N2O emissions, and growth rates of N2O emissions and N2O mixing ratio.The top panel shows annual inputs for fertilisation, deposition, fixation, and total N used in the model (see Section 4.1 for data sources).The second panel shows the anthropogenic flux broken down into N input categories of fertilisation, deposition and fixation, estimated by assuming that all increases in N2O emissions for all input categories (fixation, fertilisation, deposition) after 1850 are due to anthropogenic influences.The shaded areas indicate the 1σ uncertainty.The third panel shows the growth rate of N2O emissions (Tg N2O-N a −2 ) from each input category calculated over 3 (pale) and 10 (dark) year windows.The bottom panel shows the growth rate in modelled N2O tropospheric mixing ratio (left axis; purple solid line) as well as the modelled change in mixing ratio growth rate (right axis; blue dotted line).
Figure 2).Between 1945 and 1980 the 10-year mean growth rate in emissions increased to 0.12±0.02Tg N 2 O-N a −1 a −1 due to rapid growth in fertiliser inputs, after which it stabilised at around 0.1 Tg N 2 O-N a −1 a −1 , corresponding to a steady rate of change of around 0.015 nmol mol −1 a −1 a −1 in tropospheric background N 2 O growth rate.The growth rate of anthropogenic emissions was particularly high between 2010 and 2015 (3-year growth rate up to 0.3 Tg N 2 O-N a −1 a −1 ) compared to the average of the last 50 years (0.09±0.06 Tg N 2 O-N a −1 a −1 ), in agreement with two recent studies (R. L. Thompson

Figure 3 .
Figure 3. Modelled N2O emission factor: a) Global gridded N2O emission factor (upper) and the 1σ uncertainty (lower); b) Relationship between mean annual precipitation (MAP), mean annual temperature (MAT; point colour) and N2O emission factor for each grid cell.
EFs to facilitate upscaling and extrapolation of data from traditional methods.Modelled EFs were compared to the compiled values for croplands used in the metaanlysis ofCui et al. (2021) by finding the mean and standard deviation of all EFs reported byCui et al. (2021) in each gridcell.Field EF measurements are laborious and challenging, thus relatively few data points are available -of the 55 478 gridcells with a valid modelled EF, only 179 have one or more EF measurements.The agreement between modelled values and observations was relatively good, with a Spearman correlation coefficient of 0.4 (p < 0.01) indicating moderate agreement.However, the slope of 0.15 suggested that measurements may consistently underestimate EFs due to insufficient measurement frequency and duration and the importance of 'hot spots' and 'hot moments' for annual N 2 O emission totals (Butterbach-Bahl, Baggs, Dannenmann, Kiese, & Zechmeister-Boltenstern, 2013; McDaniel, Kaye, & Kaye, 2014).The impact of climate change -in particular changing precipitation patterns -on N 2 O emissions in rapidly developing regions like sub-Saharan Africa and India is not yet captured by observations, and should be a focus of future studies.Understanding the large-scale impact of moisture availability and other climate parameters on N loss processes through targeted measurements campaigns and model development will be key to predicting interactions between climate change-driven precipitation changes and N 2 O emissions.The proportion of N 2 O lost to different pathways shows clear temporal changes, particularly for fertilisation N inputs, which show an increase in mean EF of N 2 O from <4% prior to 1940to >5% in 2020(Figure A.14).This has two major causes: Changing spatial distribution of N inputs and climate warming feedback (Figure4a).The shift in fertilisation from North America and other temperate regions (EF N2O of 4-5%) towards emerging economies in warmer regions with higher EFs (EF N2O >7.5%), in particular China, has caused additional emissions of 0.5±0.3Tg N 2 O-N a −1 since 1940 (Figure4c).This effect is only seen for fertilisation N inputs: Deposition and fixation inputs show vary minor spatiotemporal manuscript submitted to Nature Geoscience changes (Figure4c).Deposition inputs are stabilising in China due to vigorous controls on N pollution (G.Yu et al., 2019), and changes in N 2 O from deposition in coming decades are likely to be of minor importance compared to fertilisation-and climate-driven emissions.The impact of climate warming (Figure4b) has led to an increase in total microbial gas production of around 3% between 1850 and 2020(Figure A.14), resulting in additional emissions of 0.8±0.4Tg N 2 O-N a −1 , and a consequent decrease in leaching losses.The parameterised 'warming' impact will be a combination of feedbacks driven by both changing precipitation and increased temperature, which cannot be separated in the current IsoTONE framework.The combined increase of 1.3 Tg N 2 O-N a −1 from spatial input changes and climate warming represents 18% of soil and 15% of total anthropogenic N 2 O emissions in 2020.As climate warming continues, and fertiliser use increases in many tropical and subtropical regions, the mean EF of N 2 O and thus the growth rate of tropospheric N 2 O will accelerate, unless significant efforts are focussed on N 2 O mitigation and increased fertiliser N use efficiency in emerging economies(Tian et al., 2020).Moreover, the impact of potential unknown climatic feedbacks, in particular non-linear responses associated with extreme events such as drought and flooding, should be a key focus for both measurement and model studies.

Figure 4 .
Figure 4. Impact of changing spatial distribution of fertilisation as well as climate warming on N2O emissions from 1940 to 2020: a) The average change in annual N2O emissions due to the shifting spatial distribution of fertiliser N inputs, calculated using the 1940 spatial distribution of fertilisation scaled to the 2020 total fertilisation quantity, and comparing this to the actual 2020 emissions.b) The average change in annual N2O emissions due to climate warming and the temperature dependency of N2O production, found by calculating 2020 emissions with and without the impact of warming on EF since 1940.c) The proportion of anthropogenic N inputs from fertilisation, deposition and fixation in 1940, 1980 and 2020 accounted for within bins defined according to the EF for N2O.20 bins were used; 19 were distributed evenly between the parameter minimum and the mean + 3 standard deviations; the highest bin was for all data >mean + 3 standard deviations.

Figure 5 .
Figure 5. Tropospheric background N2O mixing ratio and isotopic composition (δ 15 N and SP).Scale offsets between different datasets have been corrected using overlapping periods to match the Empa dataset as described in the text.Mean data at 25-year intervals (1740-1940) and 2-year intervals are shown as orange points with black outline; the standard deviation of the data averaged for each point is shown as the orange shaded region.

ParameterisationsANNFigure 6 .
Figure 6.A schematic view of the IsoTONE model and inversion structure, showing different input and output data types.Optimized variables (green) are defined in Table A.1.Sources of input and observation data are given in Section 4.1.Abbreviations: MAT = Mean Annual Temperature; WFPS = Water-Filled Pore Space; CABLE = Community Atmosphere Biosphere Land Exchange land surface model; MCMC = Markov Chain Monte Carlo; fG, fL, fNO, fN2O, and fN2 refer to the proportion of N lost through leaching and N gas production, discussed in Section 4.2.1.
contrast to bacterial and abiotic denitrification, fungal denitrification results in N 2 O with a high and variable SP.The partitioning of N 2 O production into fungal denitrification as compared to bacterial nitrification and denitrification (shown in Figure A.6) is currently not known and therefore cannot be parameterised within this model.Fungal denitrification is a minor source of N 2 O, expected to contribute less than 10-15% globally (Rohe et al., 2020), and will be included in a future version of this model as data to drive parameterisations becomes available.These fractionation factors refer to the immediate production or consumption of N substrates and N 2 O, and thus do not account for the complexity of processes leading to the net isotopic composition of emitted N 2 O in real environments.Several studies

(
2020);Olivier and Berdowski (2001)) gridded emissions for the major non-soil anthropogenic N 2 O emission categories (1A1 = power industry, 1A3b = road transport, 2B = chemical processes, and 6 = wastewater treatment) were added to estimate total terrestrial N 2 O emissions.Isotopic composition was estimated as δ 15 N = 3.9±2.9and SP = 17.6±0.5 for power industry, δ 15 N = -7.2±1.2 and SP = 10.0±4.3 for road transport, δ 15 N = -8.3±10.6 and SP = 3.3±5.5 for chemical processes, and δ 15 N = -11.6±12.7 and SP = 10.0±5.7 Nature Geoscience to achieve stable results with no significant difference between posterior parameters by step size, and no change in posterior parameters following more iterations.All tested and accepted parameter sets were saved; tested parameter sets were used to check coverage of the parameter uncertainty space (Figure A.7), and accepted parameter sets were used to find posterior results and uncertainties for the parameters (Table Figure A.1.Relationships between δ 15 N soil and ancillary data.a) Color map showing linear correlation coefficients (P = Pearson, S = Spearman) indicated with the color bar between δ 15 N soil and the specified parameter.b) Principal components analysis of the combined δ 15 N soil and ancillary data; soil parameters are shown in blue and climate parameters in orange.
0005 and root mean squared error (RMSE) as a loss function.During training, a validation split of 25% was used to minimise overfitting, and early stopping was implemented based on the minimum loss function for validation data with a patience level of 30 epochs.A model checkpoint was used to save the model that achieved the lowest RMSE for the validation data during the training period.Training was carried out over a maximum of 80 epochs with a batch size of 200, and data shuffling was activated.As samples were not distributed evenly across different regions, sample weighting to reduce the relative importance of samples from highly represented regions such as Europe was needed.Samples were grouped into 30 clusters according to latitude and longitude using the KMeans function from the 'sklearn' package.A sample weight vector was calculated as 1 √ n where n is the number of members in a cluster.The sample weights were used to weight the loss function during the ANN training.
representation uncertainty was estimated to be 1.2❻, based on the mean of the standard deviation of all measurements placed within each grid cell for the 461 grid cells with more than 2 measurements.250 bootstrap iterations were run, and within each iteration, data for δ 15 N soil was perturbed according to the representation uncertainty multiplied by normally distributed random numbers.δ 15 N soil and the chosen input parameters (latitude, absolute longitude, MAP, MAT, soil N, soil C, soil C:N, pH, bulk density, and aridity index) were randomly split into training ( 2 3 of data) and testing subsets.Data was normalised to a range of 0-1 for each parameter, and all rows containing NAs were removed.6000 rows of data were used for training and 3000 for testing the ANN in each iteration.The trained ANN was used to predict the test dataset values of δ 15 N soil .The predicted and observed test δ 15 N soil were compared using the weighted root mean square error (wRMSE) and the R 2 of a weighted linear fit within each bootstrap iteration, using 1 √ n where n is the number of members in the k-means cluster size as the weight for each point.The ANN achieving the best fit (lowest wRMSE and highest R 2 ) on the test data was used for the final prediction of δ 15 N soil , and the standard deviation of the predicted δ 15 N soil from all bootstrap iterations was used to estimate the uncertainty.

Figure
Figure A.2. a) Global gridded (0.5 • × 0.5 • grid) δ 15 N soil estimated using an artificial neural network.Note the colour scale is non-linear to better show the narrow range of most values.b) Estimated uncertainty in gridded δ 15 N soil .
Figure A.3.Relationships between N2O fluxes and emission factors (EF) and ancillary data.a) Color map showing linear correlation coefficients (P = Pearson, S = Spearman) for fluxes and EFs indicated with the color bar.b) Principal components analysis of the combined soil flux, EF and ancillary data; soil parameters are shown in blue, climate parameters in orange, and N-cycle parameters in red.
Figure A.4. Mean emission factors for different climate zones based on measured EF data binned according to MAT and soil bulk density.a) Mean EFs for the different climate zones; MAT (mid-point of bin) is plotted on the x-axis and line colours refer to the bulk density bin.The standard deviation of values in each bin is shown with the shaded region and the number of values in each bin is annotated.b) Global map showing the spatial coverage of the 16 bins: MAT increases with bin number, and bulk density increases sequentially through bins 1-4, 5-8, 9-12 and 13-16.
Figure A.5. Parameterisation of N2O/NO/N2 as a fraction of total N gas emissions based on water-filled pore space.The upper panels show experimental data (includes all data used by Bai et al. (2012) as well as additional data from Harris et al. (2021); Ibraim et al. (2019); Verhoeven et al. (2019); Wolf et al. (2015); Xia et al. (2019)).Prior estimates of the parameterisation using sigmoid fits to the data are shown as dotted lines, with the prior 1σ uncertainty range shown with pale shading.Posterior fits following the MCMC optimization are shown as solid lines with the posterior 1σ uncertainty range in darker shading.Eqs.A.1, A.2 and A.3 describe the prior and Eqs.A.3, A.5 and A.6 the final fits.The bottom panel shows the proportion of total N gas contributed by NO, N2O and NO according to the parameterisation used by Bai et al. (2012) (dashed line), the prior fits (dotted line) and the posterior fits (solid line with uncertainty range shaded).

Figure A. 6 .
Figure A.6.Parameterisation of denitrification and nitrification contribution to N2O emissions based on water-filled pore space.The upper panel shows experimental data for denitrification N2O as a proportion of total N2O emissions from Harris et al. (2021); Ibraim et al. (2019); Verhoeven et al. (2019); Wolf et al. (2015), with a sigmoid fit to the data.The lower panel shows N2O (Eqs. 40 000 iterations were run for each step size 0.25, 0.5 and 0.75, with 1 760 (4.4%), 76 (0.2%) and 34 (0.1%) solutions accepted respectively(Figure A.7 and Table A.1).There were no significant differences in the means of accepted values for the different stepsizes, and no change in the final values if only a subset of the accepted values randomly sampled, showing that the posterior parameters were robust.The posterior parameters showed a strong reduction in uncertainty compared to the prior for most parameters (Figure A.7 and Table A.1), although τ PI /τ PD , τ PD and temp sens -which are all primarily constrained by N 2 O mixing ratio increase -cannot be easily separated during model-data assimilation.F ocean , τ PD and T to S showed the strongest correlations between accepted solutions and dominated PC1 in the PCA (Figure A.8), again reflecting their role in regulating N 2 O mixing ratio increase in the atmosphere submodule.Addition of an independent tracer of T to S, or very high precision measurements from multiple sites in separate hemispheres, may help constrain these parameters in future model studies.PC2 reflected the role of temperature increase, while PC3 was dominated the chosen fractionation factors, which drive loss partitioning in the soil submodule through the parameter frac ex.Both the correlation and PCA analysis show that most of the optimized parameters -particularly those related to isotopic composition -can be sufficiently discriminated by the model-data assimilation.Parameters only constrained by mixing ratio cannot be distinguished easily, reflecting the importance of isotopic data for constraining the global N cycle.
Figure A.7. Tested and accepted model solutions from the MCMC optimization of a coupled soil-atmosphere model of N2O emissions and isotopic composition.Parameter descriptions and references for prior estimates are given in Table A.1.
Figure A.8. Relationships between accepted solutions of parameters optimized with MCMC.a) Color map showing Pearson correlation coefficients between the specified parameters as indicated with the color bar.Only correlations significant at p < 0.01 are shown.b) Principal components analysis of all accepted solutions; soil and emissions module parameters are shown in orange, atmospheric parameters in blue, and isotopic parameters in red.
, which estimated total terrestrial denitrification N losses of 90-135 Tg-N a −1 .Our estimate of the mean terrestrial denitrification product ratio (R N 2O = N2O denit N2O denit +N2 ) of 27±4% in 2020 (26±3% in 1850) agrees very well with the estimate of 23% based on a metaanalysis of experimental studies by Schlesinger (2009).Wetland and estuarine soils are not represented in the δ 15 N database used as the basis for modelling in this study, suggesting these regions may be hotspots for complete denitrification and N 2 production, accounting for the discrepancy in modelled total denitrification.Total denitrification N losses are highly uncertain due to the technical challenge of making accurate annual measurements of N 2 fluxes, thus model validation extremely difficult (e.g.

Figure
FigureA.6).An improved quantitative understanding of the drivers of denitrification, as well as measurements of soil δ 15 N in wetland and estuarine environments, will be key to achieving robust model estimates of total global denitrification.
Figure A.11. Anthropogenic and natural N2O emissions (top panel) and isotopic composition (δ 15 N and SP) of natural and anthropogenic terrestrial N2O emissions.The anthropogenic source was estimated by assuming that all increases in N2O emissions for all input categories (fixation, fertilisation, deposition) as well as all changes in isotopic composition after 1850 are due to anthropogenic influences.The shaded areas indicate the 1σ uncertainty.

Figure A. 12 .
Figure A.12. Modelled terrestrial (natural + anthropogenic) N2O soil emissions from different input types.The left panels show total emissions for inputs from biological N fixation, N deposition, and N fertilisation for the year 2020; all three subpanels have the same colour scale.The right panels show anthropogenic emissions for each input type in 2020, estimated by subtracting total emissions for the year 1800 from total 2020 emissions.All subpanels have the same colour scale.
Total gas production accounts for the largest proportion of N losses in warm, dry regions, whereas N 2 O EFs are highest in non-desert tropical regions, in particular sub- • latitude and -180 • to 180 • longitude.The δ 15 N grid (Section 4.1.2;Figure A.2a) was first initialised in each model run by adding 5% of the uncertainty in δ 15 N (Figure A.2b) multiplied by a normally distributed random number (µ = 0, σ = 1).For each grid cell, the soil model was initialised with an available soil nitrogen pool of size 1 (unitless) with δ 15 N of 0❻, and a soil N 2 O pool of of size 1 (unitless) with δ 15 N and SP of 0❻.The 14 N input rate (k i ) was set at 1 (unitless) and the 15 N input rate (k i,15 ) calculated using δ 15 N i of 0.5❻ (based on Bai et al. (2012); Beyn, Matthias, Aulinger, and Dähnke (2015); Unkovich (2013); Zhang, Liu, Fangmeier, Goulding, and Zhang (2008)).Changing the unitless pool sizes affects only the number of iterations until steady state is reached in the model, and not the final result of the model.Reducing δ 15 N i by 1❻ increased the calculated mean global f G by ∼ 1%, while increasing δ 15 N i by 1❻ increased the calculated mean global f G by ∼ 5%, with very little impact of spatial distribution of f G .
(Craine, Elmore, et al., 2015)uscript, IsoTONE model code (written in Python) will be available open access on Github.Gridded datasets calculated in this study (eg.δ 15 N soil , f gas , N 2 O emission factors, gridded δ 15 N and SP of emissions) will be released together with the model code.δ15N point data not present in the(Craine, Elmore, et al., 2015)dataset will be published as a compilation within 1-2 years, once associated research projects of the individual data owners are completed.Tropospheric background N 2 O isotopic data collected at Empa and used for model optimization will be released open access in 2022 together with the associated manuscript (L.Yu et al., 2022).Using isotopes to trace the effects of climate extremes on N 2 O emissions and the nitrogen cycle in managed grasslands' and the Swiss National Science Foundation (SNF) project 163075 'Assessment of the global N 2 O budget based on seasonal and longterm isotope measurements at Jungfraujoch and the Cape Grim Air Archive'.An OECD Cooperative Research Program for Sustainable Agricultural and Food Systems (OECD-CRP) grant for the project 'Identifying drivers of N 2 O emissions in a changing climate' supported Eliza Harris' research stay at the University of Melbourne to carry out this work.Longfei Yu was supported by the EMPAPOSTDOCS-II program, which receives funding from the European Union's Horizon 2020 research and innovation program under the Marie Sklodowska-Curie grant agreement number 754364.Geoclimatic gradients in δ 15 N soil and in N 2 O emission factors A.1.1 δ 15 N soil The N 2 O observations at Jungfraujoch were supported by the International Foundation High Altitude Research Stations Jungfraujoch and Gornergrat (HFSJG), the Swiss Federal Office for the Environment, and ICOS-CH (Integrated Carbon Observation System Research Infrastructure), and we would like to thank Christoph Zellweger (Empa) for assistance with data collection from these sites.We would also like to thank Eric Slessarev (Lawrence Livermore National Laboratory) for provision of soil pH data, and Kate Summer and Matheus Carvahlo de Carvahlo (Southern Cross University) for assistance with sample preparation and data collation for Australian soils in the BASE database.The inverse modelling results were provided by R. L. Thompson who was funded through the Copernicus Atmosphere Monitoring Service (https://atmosphere.copernicus.eu/),implemented by ECMWF on behalf of the European Commission, and were generated using computing resources from LSCE.