Constraining N2O emissions since 1940 using firn air isotope measurements in both hemispheres

. N 2 O is currently the third most important anthropogenic greenhouse gas in terms of radiative forcing and its atmospheric mole fraction is rising steadily. To quantify the growth rate and its causes over the past decades, we performed a multi-site reconstruction of the atmospheric N 2 O mole fraction and isotopic composition using new and previously published

Abstract. N 2 O is currently the third most important anthropogenic greenhouse gas in terms of radiative forcing and its atmospheric mole fraction is rising steadily. To quantify the growth rate and its causes over the past decades, we performed a multi-site reconstruction of the atmospheric N 2 O mole fraction and isotopic composition using new and previously published firn air data collected from Greenland and Antarctica in combination with a firn diffusion and densification model. The multi-site reconstruction showed that while the global mean N 2 O mole fraction increased from (290 ± 1) nmol mol −1 in 1940 to (322 ± 1) nmol mol −1 in 2008, the isotopic composition of atmospheric N 2 O decreased by (−2.2 ± 0.2) ‰ for δ 15 N av , (−1.0 ± 0.3) ‰ for δ 18 O, (−1.3 ± 0.6) ‰ for δ 15 N α , and (−2.8 ± 0.6) ‰ for δ 15 N β over the same period. The detailed temporal evolution of the mole fraction and isotopic composition derived from the firn air model was then used in a two-box atmospheric model (comprising a stratospheric box and a tropospheric box) to infer changes in the isotopic source signature over time. The precise value of the source strength depends on the choice of the N 2 O lifetime, which we choose to fix at 123 years. The average isotopic composition over the investigated period is δ 15 N av = (−7.6 ± 0.8) ‰ (vs. air-N 2 ), δ 18 O = (32.2 ± 0.2) ‰ (vs. Vienna Standard Mean Ocean Water -VSMOW) for δ 18 O, δ 15 N α = (−3.0 ± 1.9) ‰ and δ 15 N β = (−11.7 ± 2.3) ‰. δ 15 N av , and δ 15 N β show some temporal variability, while for the other signatures the error bars of the reconstruction are too large to retrieve reliable temporal changes. Possible processes that may explain trends in 15 N are discussed. The 15 N site preference (= δ 15 N α − δ 15 N β ) provides evidence of a shift in emissions from denitrification to nitrification, although the uncertainty envelopes are large.

Introduction
The rise in nitrous oxide (N 2 O) since pre-industrial times has contributed significantly to radiative forcing (Forster et al., 2007). Over the past 4 decades, the N 2 O mole fraction has increased by 0.25 % per year, reaching 324 nmol mol −1 in 2011 (Ciais et al., 2013). Therefore, the understanding of the biogeochemical cycle of N 2 O is important for a reliable assessment of future climate change. In addition, the destruction of N 2 O in the stratosphere provides an important source of nitrogen oxides (NO x ), which contribute to stratospheric ozone depletion (Ravishankara et al., 2009;Crutzen, 1979;McElroy et al., 1971).
Natural sources of N 2 O are microbial processes in soils and oceans, which produce N 2 O during nitrification and denitrification (Bouwman et al., 2013;Löscher et al., 2012;Santoro et al., 2011;Galloway et al., 2004;Pérez et al., 2001;Yung and Miller, 1997;Kim and Craig, 1993). The increase in N 2 O since pre-industrial times (hereafter referred to as "anthropogenic" increase) has been attributed largely to increased microbial production, resulting from the increased use of nitrogen fertilizers in agriculture. Industry (especially nylon production) and fossil fuel combustion present a smaller contribution to the anthropogenic source (Davidson, 2009;Kroeze et al., 1999;Mosier et al., 1998). N 2 O is primarily destroyed in the stratosphere via UV photolysis (90 %) and reactions with excited oxygen atoms (10 %) (Minschwaner et al., 1993), with a minor N 2 O fraction removed by surface sinks (Syakila et al., 2010).
Estimates of the total N 2 O source strength from various bottom-up and top-down studies suggest a mean value of roughly 17 Tg year −1 N equivalents at present. However, the range in both approaches is large, especially for bottom-up estimates, which range between 8.5 and 27.7 Tg year −1 N, whereas top-down estimates range between 15.8 and 18.4 Tg year −1 N (Potter et al., 2011, and references therein). Besides the total source strength, the contributions of individual source processes are also poorly constrained. Due to the long steady-state lifetime of N 2 O in the atmosphere (123 years; SPARC Lifetimes Report, 2013), temporal and spatial gradients are small, making it difficult to resolve localized sources.
Measurements of the isotopic composition of N 2 O may help to constrain the atmospheric N 2 O budget. The N 2 O molecule is linear (NNO) and the two N atoms are chemically distinguishable; as a consequence they tend to attain different isotopic compositions. Beyond oxygen (δ 18 O, δ 17 O) and average δ 15 N av ("bulk") signatures, N 2 O also displays site-specific 15 N isotopic information. Site preference (δ 15 N sp ) is defined as the difference in δ 15 N av between the central (2, µ, or α) and terminal positions (1, τ , or β) of N atoms in N 2 O (Kaiser et al., 2002;Brenninkmeijer and Röckmann, 2000;Yoshida and Toyoda, 2000), i.e., δ 15 N sp = δ 15 N α − δ 15 N β . For consistency with many recent publications in the field, we adopt here the nomenclature from Yoshida and Toyoda (1999), α and β, for the two positions.
The different sources and sinks of N 2 O are associated with characteristic fractionation processes leading to different isotope ratios. For example, microbial sources emit N 2 O that is depleted in 15 N and 18 O relative to the tropospheric background. N 2 O that returns from the stratosphere after partial photochemical removal is enriched in both heavy isotopes (Yoshida and Toyoda, 2000;Rahn and Wahlen, 1997;Yung and Miller, 1997;Kim and Craig, 1993). Stratospheric N 2 O also has a high 15 N site preference compared to tropospheric N 2 O. The observed enrichment is caused by kinetic isotope fractionation in the stratospheric sink reactions (Kaiser et al., , 2002Park et al., 2004;Röckmann et al., 2001;Yoshida and Toyoda, 2000;Rahn et al., 1998).
The multi-isotope signature of N 2 O adds useful constraints on its budget. In particular, when the isotopic composition of tropospheric N 2 O is combined with the fractionation during its removal in the stratosphere, the isotopic composition of the global average source can be determined (Ishijima et al., 2007;Bernard et al., 2006;Röckmann et al., 2003b;Kim and Craig, 1993).
The temporal variations of the N 2 O isotopic composition are difficult to quantify on a short timescale because of its long residence time in the atmosphere. Longer timescales can be reconstructed by using air trapped in Arctic and Antarctic firn and ice which provides a natural archive of past atmospheric composition. The firn phase is the intermediate stage between snow and glacial ice, which constitutes the upper 40-120 m of the accumulation zone of ice sheets. Within the firn, air exchanges relatively freely in the upper layers and with the overlying atmosphere (convective zone). With increasing depth the air pores shrink in size due to firn compaction, and air mixes primarily via slow diffusion in the diffusive zone. At densities larger than ≈ 815 kg m −3 , air is permanently trapped in closed bubbles in the ice and totally isolated from the atmosphere. The precise age range of air that can be retrieved from polar firn between the surface and bubble close-off depends on site-specific characteristics like temperature, accumulation rate, and porosity and typically ranges from several decades to 120 years.
For N 2 O, a number of studies have reported isotope measurements from different Arctic and Antarctic firn drilling sites, showing a steady decrease in the heavy isotope content of N 2 O over the past decades Ishijima et al., 2007;Bernard et al., 2006;Röckmann et al., 2003b;Sowers et al., 2002). A more recent study by Park et al. (2012) attempted to reconstruct the long-term trends in N 2 O isotopic compositions and its seasonal cycles to further distinguish between the influence of the stratospheric sink and the oceanic source at Cape Grim, Tasmania, demonstrating that isotope measurements can help in the attribution and quantification of surface sources in general.
Taking into account the long atmospheric lifetime of N 2 O and the fact that both hemispheres are well mixed on annual timescales, it is reasonable to assume that the results from these studies are representative of the global scale. However, care needs to be taken because small differences in the diffusivity profiles of the firn column lead to large effects on the isotope signature (Buizert et al., 2013). Interestingly, for atmospheric methane (CH 4 ), another important greenhouse gas, a recent multi-site analysis of its carbon isotopic compo-  Bernard et al. (2006) and Ishijima et al. (2007). 2 Data retrieved from Bernard et al. (2006). 3 Data retrieved from Röckmann et al. (2003b).
sition showed large differences among reconstructions from different sites (Sapart et al., 2013). In particular, firn fractionation effects related to diffusion and gravitational separation and their implementation in models (Buizert et al., 2012) have large effects on the reconstructed signals. Small differences in the diffusivity profiles of the firn column lead to large effects on the isotope signatures. Therefore, more robust results may be obtained by combining isotope information from a number of different sites in a multi-site reconstruction, including a critical evaluation of diffusivity profiles.
Here we combine new N 2 O isotope measurements from the NEEM site in Greenland with previously published firn air N 2 O isotope records from four different sites from Greenland and Antarctica to reconstruct records of the N 2 O isotopic composition over the last 70 years. We use the multigas firn transport model developed by the Laboratoire de Glaciologie et Géophysique de l'Environnement and Grenoble Image Parole Signal Automatique (LGGE-GIPSA) to obtain an atmospheric scenario that is constrained by and consistent with all individual sites (Allin et al., 2015;Wang et al., 2012;Rommelaere et al., 1997). We then use an isotope mass balance model to infer the changes in the isotopic signature of the N 2 O source over time to investigate possible changes in the source mix.

Firn air sampling
New firn air samples added in this study to the total dataset were collected in 2008 and 2009 during the firn campaign (Buizert et al., 2012) as part of the North Eemian Ice Drilling programme (NEEM) in Greenland (77.45 • N,51.06 • W). These data are combined with existing firn air data from four other sites. Information on the locations is provided in Table 1. The firn air collection procedure is described in detail by Schwander et al. (1993). Here a brief description is presented. Essentially a borehole is drilled in the firn to a certain depth and then the firn air sampling device is inserted into the borehole. The device consists of a bladder, a purge line, and a sample line. When the sampling device reaches the desired depth the bladder is inflated to seal the firn hole and isolate the air below the bladder from the overlying atmosphere, and air is pumped out from the pore space below the bladder. Continuous online CO 2 concentration measurements are performed to verify that no contamination with contemporary air occurs during the extraction procedure. After the contaminating air has been pumped away, firn air is collected in stainless steel, glass, or aluminum containers.

N 2 O isotope analysis
The firn air samples from NEEM are analyzed for N 2 O isotopocules at the Institute for Marine and Atmospheric research Utrecht (IMAU). The N 2 O mole fraction and isotopic composition are measured using continuous flow isotope ratio mass spectrometry (IRMS). The method is described in detail by Röckmann et al. (2003b). Here only a brief summary is given. The firn air sample (333 mL) is introduced into the analytical system at a flow rate of 50 mL min −1 for 400 s. After CO 2 is removed chemically over Ascarite, N 2 O and other condensable substances are cryogenically preconcentrated. After cryo-focusing the sample the remaining traces of CO 2 and other contaminants are removed on a capillary GC column (PoraPlot Q, 0.32 mm internal diameter (i.d.), 25 m). The column is separated into a pre-column and an analytical column. This setup eliminates interferences from other atmospheric compounds that have much longer retention times. Finally the sample is transferred to the IRMS via an open split interface. For the new NEEM samples reported here, each firn air sample has been measured five times. Before and after each sample we measured five aliquots of air from a reference cylinder with known isotopic composition and mole fraction for calibration purposes. δ 15 N av values are reported with respect to air-N 2 , while δ 18 O refers to Vienna Standard Mean Ocean Water (VS-MOW). As a laboratory reference gas we used an atmospheric air sample with an N 2 O mole fraction of 318 nmol mol −1 and δ values of (6.4 ± 0.2) ‰ for δ 15 N av vs. air-N 2 , and (44.9 ± 0.4) ‰ for δ 18 O vs. VSMOW. The intramolecular δ 15 N av values of the air standard are δ 15 N α = (15.4 ± 1.2) ‰ and δ 15 N β = (2.7 ± 1.2) ‰. The calibration of the intramolecular distribution follows Toyoda and Yoshida (1999). Typically the 1σ standard deviations of replicate sample measurements are 0.1 ‰ for δ 15 N av , 0.2 ‰ for δ 18 O, and 0.3 ‰ for δ 15 N α and δ 15 N β .

Modeling trace gas transport in firn
In firn air, the interstitial gas is not yet isolated in closedoff bubbles, so diffusion processes and gravitational separation alter mole fractions and isotope ratios over time. Thus, firn air measurements cannot be used directly to derive the atmospheric history of trace gas signatures. Over time, atmospheric compositional changes are propagated downwards into the firn based on the diffusivity of the atmospheric constituent in question. Firn air diffusion models take these effects into account and thereby allow reconstruction of changes in the atmospheric composition from the firn profile.
In this study we use the LGGE-GIPSA firn air transport model to reconstruct the temporal evolution of N 2 O mole fraction and isotopic composition from the measured firn profiles (Allin et al., 2015;Witrant et al., 2012;Wang et al., 2012;Rommelaere et al., 1997).
In the "forward version" of LGGE-GIPSA, a physical transport model uses a historic evolution of atmospheric N 2 O mole fractions to calculate the vertical profiles of mole fractions in firn. For the isotopocules, further simulations are performed separately to calculate their respective vertical profiles. Important parameters needed to constrain the model are the site temperature, accumulation rate, depth of the convective layer, and close-off depth, together with profiles of firn density and effective diffusivity. The latter parameter is determined as a function of depth for each firn-drilling site by modeling the mole fractions in firn for trace gases with wellknown atmospheric histories (Buizert et al., 2012;Rommelaere et al., 1997;Trudinger et al., 1997). A multi-gas constrained inverse method  is used to calculate the effective diffusivity of each site for each specific gas. It is noteworthy that diffusivity is not constrained equally well at all sites because different sets of constraints (e.g., number of available reference gases) are used at different sites and because of different depth resolutions.
A Green function approach, as presented by Rommelaere et al. (1997) and used for halocarbon trend reconstruction by Martinerie et al. (2009), with an extension for isotopic ratios and revised to take into account the scarcity of the measurements (Allin et al., 2015;Witrant and Martinerie, 2013;Wang et al., 2012), is used to assign a mean age and age distribution to a certain depth.
Due to the long N 2 O residence time in the atmosphere, the global variability of the isotopic composition of N 2 O is very small and no significant variations between individual background locations have been detected so far (Kaiser et al., 2003). In particular, the isotope ratio difference between Northern Hemisphere and Southern Hemisphere tropospheric air is expected to be only −0.06 ‰ (based on an interhemispheric mole fraction gradient of 1.2 nmol mol −1 , Hirsch et al., 2006, and an isotope ratio difference of −15 ‰ between average source and average tropospheric isotopic delta values). These differences are within the uncertainties of the firn air measurements used here, and therefore the data from the Northern Hemisphere and Southern Hemisphere are combined into a single dataset without including an interhemispheric gradient.
With the multi-site reconstruction method, we used the measurements from six firn air drillings at five sites (NEEM-09, NEEM-EU-08, BKN-03, NGRIP-01 Bernard , DC-99, DML-98) to constrain our model and determine a set of atmospheric reconstructions that fits all sites. Data from Ishijima et al. (2007) and Sowers et al. (2002) (NGRIP-01 Ishijima , SP-01, and SP-95, respectively) were not included in our multi-site reconstruction because no data for δ 15 N α and δ 15 N β were published for those sites. These datasets were used for independent validation of δ 15 N av and δ 18 O.
To quantify the isotope fractionation due to diffusion and gravitational settling within the firn, a forward firn transport model simulation was carried out with a realistic N 2 O mole fraction scenario (based on the Law Dome record, MacFarling Meure et al., 2006), but with a constant isotopic N 2 O history. This allows determination of the role of transport isotope fractionation occurring in the firn, in the absence of isotopic changes in the atmosphere. The results are used to subtract the firn fractionation effects from the measured signals, which can then be used to assess the atmospheric history. Compared to the signal, the effect of firn fractionation is minor for δ 15 N av but important for δ 18 O, especially at the lower accumulation rates in the Southern Hemisphere (see Appendix A).
The deepest firn data from each site provide constraints furthest back in time and the oldest air samples that are included in the inversion are from the DML-98 and DC-99, which extend the reconstruction of atmospheric N 2 O back to the early 20th century (Röckmann et al., 2003b). At the same time, the correction for isotopic fractionation in firn is most uncertain for the deepest samples, where strong differences between individual firn air models have been reported (Buizert et al., 2012).

Scaling of different datasets
At present, no international reference materials for the isotopic composition of N 2 O exist. Kaiser et al. (2003) and Toyoda and Yoshida (1999) linked the isotopic composition of N 2 O in tropospheric air to the international isotopes scales for nitrogen isotopes (air-N 2 ) and oxygen isotopes (either VSMOW or air-O 2 ). Our measurements are linked to a standard gas cylinder of tropospheric air with known N 2 O mole fraction and isotopic composition based on the scale of Kaiser et al. (2003) for δ 15 N av and δ 18 O values and Yoshida and Toyoda (1999) for position-dependent 15 N values. However, the reference air cylinder used for the calibration was exhausted and had to be replaced three times over the years in which the different measurements that we combine in this study were performed. Although the cylinders were carefully compared, the long-term consistency of the isotope scale could not be guaranteed because long-term isotope standards are not available (Table 2). In fact, analysis of the data from the convective zone for the different sites shows small but significant differences from the temporal trends that are well established from previously published data from the German Antarctic Georg von Neumayer station for 1990 to 2002 (Röckmann and Levin, 2005). The linear trends reported in that paper are (−0.040 ± 0.003) ‰ year −1 for δ 15 N av , (0.014 ± 0.016) ‰ year −1 for δ 15 N α , (−0.064 ± 0.016) ‰ year −1 for δ 15 N β and (−0.021 ± 0.003) ‰ year −1 for δ 18 O. Since they were derived from direct air samples (unaffected by firn fractionation), these trends can be used as a reference to re-scale the different firn air results from different dates. To do so, data from the diffusive zone (ρ < 815 kg m −3 ) for each individual site were scaled to one reference site, DC-99, taking into account the temporal differences in sampling and the model-assigned mean age of the firn air samples (see below). DC-99 was chosen as the reference site because it has the most measurements in the diffusive zone. Also, the precision of these measurements was high because high-volume cylinders were available from which many measurements could be performed and averaged. To test the sensitivity to the choice of reference site, we repeated the re-scaling using NEEM-09 as a reference, which generated almost identical results within uncertainty bars (Appendix C).
The average difference between the samples from the diffusive zone at a given site and the interpolated DC-99 results was compared to the expected temporal trend between the sampling date of each station and DC-99, using the temporal trends established by Röckmann and Levin (2005), as shown in the equations below. The effect of this scaling is that the temporal trend in the past decade is effectively forced to follow the atmospheric measurements at Neumayer station (Röckmann and Levin, 2005).
After re-scaling the firn isotopic data, we detected some individual data points that clearly deviated from the general trends. These were considered outliers, because they exceeded the 2σ error, and were removed from the dataset. All of these values are site-specific 15 N values, and specifically the following were excluded: NEEM-EU-08 hole depths  The mole fraction data that can be obtained from the NEEM air isotope measurements were substituted with more precise measurements of the N 2 O mole fraction by the Commonwealth Scientific and Industrial Research Organisation (CSIRO), the Institute of Environmental Physics, University of Heidelberg (IUP), the Centre of Ice and Climate, University of Copenhagen (CIC), and the National Oceanic and Atmospheric Administration (NOAA). In this way we combine all available N 2 O mole fraction data that narrow the uncertainty envelope but do not affect the trend.
The mole fraction data were scaled to the most recent international scale, NOAA-2006A from the CSIRO scale or the NOAA-2000 scale. Conversion of the NOAA-2000 data to the NOAA-2006A scale is done using a conversion factor available from the National Oceanic and Atmospheric Administration (NOAA) (http://www.esrl.noaa.gov/gmd/ccl/ scales/N2O_scale.html). Converting from the CSIRO to the NOAA-2006A scale, though, requires the reference cylinder details, which were not available. Instead we used a trend scenario, based on the CSIRO atmospheric scale combined with Law Dome data and assuming a constant interhemispheric gradient. This trend scenario was then compared with the data provided on the NOAA-2006A scale, and a polynomial fit was generated which was then used to convert the data to the NOAA-2006A scale. All results presented in the next section are based on the scaling procedure and removal of the outliers as described above (Appendix B).

Global N 2 O (isotope) budget calculations
The tropospheric budget is controlled by N 2 O emissions from natural and anthropogenic sources at the surface and by the exchange between troposphere and stratosphere. A simple two-box model is used to quantitatively understand the emissions and the budget changes of N 2 O. The model consists of a tropospheric N 2 O reservoir (index T) into which N 2 O is emitted from natural (E nat ) and anthropogenic (E anth ) sources. N 2 O is then transported to the stratosphere (index S) where part of it is destroyed by photochemical reactions (L), Table 3. Emission fluxes and isotopic composition of the natural and anthropogenic source results from the mass balance model. Stratospheric isotope fractionation (ε L ) used in the mass balance model, and the respective results from Park et al. (2012).
The change in the tropospheric N 2 O reservoir is given by the following mass balance equations (Allin et al., 2015): where n is the amount of air (85 % for troposphere and 15 % for stratosphere) and χ S and χ T are the mole fractions of N 2 O in the stratosphere and troposphere, respectively. Annual fluxes between the two reservoirs, F exch , are calculated based on previous estimates (Appenzeller et al., 1996;Holton et al., 1990) and given in Table 3. The loss due to strato-spheric sink is determined by where τ is the atmospheric lifetime, which we choose to fix at 123 years. The isotopic budgets are calculated by simply multiplying the reservoir sizes by the corresponding δ values of the different flux terms: Separating the l.h.s. (left-hand side) into two terms and substituting Eqs. (1) and (2) into Eqs. (4) and (5) yields the final isotope equations: where δ T is either δ 15 N av , δ 18 O, δ 15 N α , or δ 15 N β from the multi-site reconstruction as shown below. δ nat and δ anth are the isotopic composition of the natural and anthropogenic N 2 O sources, respectively (our target quantity). ε L is the isotopic fractionation factor associated with stratospheric destruction. δ S is also not known, but can be calculated using the analog from Röckmann et al. (2003b) by employing the observed apparent Rayleigh fractionation in the stratosphere (ε app ) ( Table 3). Based on this, the relative isotope ratio difference between the stratosphere and the troposphere can be calculated by Here, we used the average ε app of all lowermost stratospheric measurements from Kaiser et al. (2006) (Table 3). Note that slightly different fractionations ε app have been used in previous studies by Rahn and Wahlen (2000), Röckmann et al. (2001), and Park et al. (2012Park et al. ( , 2004. The sensitivity of the results to these differences will be examined below. Furthermore we assume that the N 2 O lifetime and ε app remained constant from the pre-industrial time to 2008; thus, the annual sink strength can be scaled down from its current value at χ T = 322 nmol mol −1 to the pre-industrial level of χ T,pi = 270 nmol mol −1 and the relative enrichment of stratospheric N 2 O relative to tropospheric N 2 O described by Eq. (8) remains constant over time. The effect of changing the N 2 O lifetime is also examined below.
Our model approach assumes that during the pre-industrial period only natural emissions occurred without any anthropogenic input. After the industrialization (≈ 1750) any increase in the emissions is considered to be due to anthropogenic input, while natural emissions remain constant. Hence, the temporal change in isotopic composition is formally due to the increase in one single "anthropogenic" source only, which may in reality also contain a natural component.

Uncertainty estimation using random scenarios
The precision of the calculated N 2 O emissions (E nat , E anth ) depends primarily on the precision of the atmospheric reconstruction of the N 2 O mole fraction (χ T ). However, taking random histories within the uncertainty envelope provided by the firn air reconstruction is not adequate to quantify the uncertainty of the atmospheric N 2 O reconstruction: the yearto-year variability of N 2 O is constrained by the N 2 O lifetime in the troposphere. Possible realistic N 2 O scenarios are scenarios that are within the confidence intervals provided by the atmospheric reconstructions, and that have realistic yearto-year variability.
Mathematically, this can be represented by an uncertainty variance covariance matrix B, where the diagonal elements (variances) are the yearly uncertainties in the atmospheric N 2 O mole fractions, and the off-diagonals are the covariances of the uncertainties of different years. The covariance between the uncertainty on the reconstruction in one year i and the uncertainty in another year j is defined as The correlation (r i,j ) is maximum between two consecutive years, and decreases as the time difference increases. We generated an ensemble of 50 random realistic N 2 O scenarios within the uncertainty envelope of the firn atmospheric N 2 O reconstruction constrained by the covariance matrix B. For each of these atmospheric N 2 O scenarios, we calculated the corresponding N 2 O emission time series. The range of emissions from these scenarios then provides a realistic estimate for the uncertainty in N 2 O emissions.
We carried out the same analysis for the different N 2 O isotopocules: for each isotopocule (δ value), we generated a covariance matrix B δ , constrained by the uncertainty ranges provided by the atmospheric reconstructions and the correlation coefficients defined in Eqs. (9) and (10) to generate a set of 50 random scenarios within the uncertainty envelopes. For each of these random scenarios, we calculated the corresponding source signature scenario, and the range in the results provides an uncertainty estimate of the isotopic source signatures. 3 Results

Mean age
The mean age of N 2 O in air sampled from different depths in the firn for all datasets that are used in this study is shown in Fig. 1. The strong change in the mean age gradient that is clearly visible in each profile reflects the transition between the diffusive and bubble close-off zones, which occurs at a specific depth and mean age for each site. Figure 1 also shows that for each site the few samples that are collected within the bubble close-off zone provide the constraints for most of the reconstructed record (for instance, at BKN-03, 50 m depth is the beginning of the bubble close-off zone). In addition to the mean age, the width of the age spectrum also increases with depth. Therefore, the temporal resolution of signals that can be reconstructed from the firn air measurements reduces with depth and approaches the one of ice core samples towards the bottom of the bubble close-off zone.
The Greenland sites (NH) have similar meteorological and glaciological conditions (Table 1); thus, the differences between the mean age profiles in Fig. 1 are small. The Antarctic sites (SH) show clear differences because the meteorological and glaciological variables differ strongly from site to site. As a result the firn-ice transition is at a different depth for each location (e.g., the firn-ice transition zone for DML-98 is located at about 73.5 m compared to about 99.5 m at DC-99).

Experimental results and multi-site reconstruction
Mole fraction and isotopic composition of N 2 O in firn air are presented vs. depth of the firn air sampling in the middle panels of Fig. 2 for the different sites. The mole fraction de-creases with depth in qualitative agreement with the known increase in N 2 O in the atmosphere over time. In contrast, all isotope deltas slowly increase with depth in the upper firn and show stronger heavy isotope enrichment in the close-off zone, both indicating heavy isotope depletion in atmospheric N 2 O with time.
The atmospheric history that has been reconstructed from these firn datasets using the multi-site inversion (using the data from NEEM-09, NEEM-EU-08, NGRIP-01 Bernard , BKN-03, DC-99, and DML-98) as described in Sect. 2.4 is shown in the left column of Fig. 2. The solid line shows the scenario that leads to the best fit with the firn data as shown in the middle panel, and the dashed lines show the upper and lower ranges of possible scenarios that would still produce an acceptable fit to the data within the uncertainty bars. Colorcoded symbols show data plotted at their respective mean age (as derived from the firn air model). When the best-fit scenario is used as input for the forward firn air model for each individual site, the model produces the vertical profiles that are shown as colored lines together with the data in the middle panels. For the sites that were included in the multi-site reconstruction, the firn profiles based on the best-fit scenarios generally match the experimental data points well, which is expected after a successful inversion procedure and with consistent datasets. The right panels in Fig. 2 show the differences between these model results and the data. For the data that were used in the multi-site inversion, the modeldata differences are generally very small, although individual firn drilling sites in some cases show small systematic deviations, in particular in the close-off zone. This means that when inversions would have been performed on individual sites, the optimal reconstructions would be slightly different. Hence, the advantage of the multi-site reconstruction is that the reconstructed scenario is constrained by all sites and all sampling depths. Despite the small differences between individual sites, the left panels show that all data fall within the uncertainty bars of the reconstructed scenario of the inversion.
From 1940 to 2008 the total changes in the δ values of atmospheric N 2 O are (−2.2 ± 0.2) ‰ for δ 15 N av , (−1.0 ± 0.3) ‰ for δ 18 O, (−1.3 ± 0.6) ‰ for δ 15 N α , and (−2.8 ± 0.6) ‰ for δ 15 N β , respectively (Fig. 2, left panels). The average linearized trends are (−0.032 ± 0.004) ‰ year −1 for δ 15 N av , (−0.014 ± 0.008) ‰ year −1 for δ 18 O, (−0.019 ± 0.015) ‰ year −1 for δ 15 N α , and (−0.041 ± 0.020) ‰ year −1 for δ 15 N β . These overall trends are slightly lower compared to previous studies that used only the data at individual sites (Ishijima et al., 2007;Bernard et al., 2006;Röckmann et al., 2003b;Sowers et al., 2002) and other studies that used data from the same period, which were not used in the present study . However, the differences are well within the combined uncertainties. We note that comparisons of average linear trends can be flawed when the firn air records have different lengths and the temporal profiles do not change linearly (see below). Trends for δ 15 N α are smaller in magnitude than for δ 15 N β , while results from Bernard et al. (2006) showed stronger changes for δ 15 N α than for δ 15 N β . However, in that study the trends were largely determined from measurements on young ice core samples with comparatively higher measurement errors and larger scatter.
Data from two sites were not included in the multi-site inversion and are used as independent validation of the reconstructed scenarios. The data points from Ishijima et al. (2007) (NGRIP-01 Ishijima , yellow) are within the range of scenarios reconstructed by the inverse model and thus independently validate our results. The δ 15 N av and δ 18 O data from Sowers et al. (2002) (SP-01 in light blue and SP-95 in blue), however, agree only for the more recent atmospheric history (Fig. 2, left panels). For mean ages before 1990 most of the points are outside the uncertainty envelopes of the multi-site reconstruction. Inter-laboratory calibration differences might be a possible explanation for the discrepancy, but the differences are not a systematic shift, and they are larger than offsets among laboratories that were established in the past (Sapart et al., 2011;Kaiser et al., 2003). In fact, the data reported by Sowers et al. (2002) were actually measured in two different laboratories with good agreement. So measurement flaws can be excluded. A possible origin of the difference could be based on the reconstruction model. Because the uncertainties in the South Pole data are large, compared to the other sites, the multi-site homogenization is more uncertain and less efficient (see Appendices A and C, and Figs. A1 and C1-C3). Sampling uncertainty should also be taken into consideration since when pumping firn air and filling the sampling flasks you could encounter uncertainties (contamination, possible leak, fractionation, incomplete flask flushing, etc.). At present though the discrepancy cannot be resolved.
To evaluate our scaling approach we repeated the multisite reconstruction using the original non-re-scaled data and re-scaled them to NEEM-09 instead of DC-99 (see Appendix C). This yielded similar results (within uncertainties) to the original reconstruction; thus, results do not depend on the choice of the site used for re-scaling. Without re-scaling, the overall change in N 2 O mole fraction and isotopic composition remained the same, but an additional decadal variability was introduced for δ 15 N av , δ 15 N α , and δ 15 N β . In addition to that, the uncertainty envelopes doubled because of the scale inconsistencies. All scaling approaches produce results that are consistent with our preferred scaling to DC-99 within the uncertainty envelopes. We conclude that scaling removed the discrepancies that would cause larger uncertainties if the original data were used instead, but the re-scaling does not introduce artificial signals (see Appendix C).
The regularization of the inversion results using a rugosity factor introduces a free parameter, which is chosen to eliminate overfitting of experimental uncertainties and which controls the smoothness of the reconstruction. The value of this parameter is set based on a robust generalized cross validation criterion, ensuring that the resolution obtained from the inverse model is similar to the experimental data while taking into account the scarcity of the measurements . A sensitivity experiment where the weight of the regularization term was increased by a factor of 10 led to nearly linear tropospheric histories within the uncertainty envelopes presented in Appendix C (Fig. C2). This combined with the fact that straight lines can be drawn within the uncertainty envelopes of the reconstructed scenarios and the sensitivity tests (see Appendix C) indicates that the isotopic trends are not significantly different from straight lines within the current uncertainties. The increase in the N 2 O mole fraction of (32 ± 1) nmol mol −1 over the reconstruction period can be explained in the mass balance model by a (4.4 ± 1.7) Tg year −1 N increase in the emissions from 1940 to 2008. The emissions increased with an increasing trend until 1975, then the annual increase continued, but at a slower rate up to 1990, and from then on the annual emissions have stayed approximately constant or even decreased slightly. The minor increase in the N 2 O mole fraction towards the end of the time series is likely not significant and does influence our reconstructions. The corresponding changes in the mole fraction are difficult to discern due to the long atmospheric lifetime of N 2 O. On average, the annual growth rate from the 1995 to 2008 period is 0.7 nmol mol −1 year −1 , corresponding to average annual emissions of 3.5 Tg year −1 N.

The temporal evolution of the N 2 O isotope signatures
The results from the isotope budget calculations are presented in Fig. 4. The left panels show the atmospheric trends. The solid black lines represent the best-fit scenarios, while the dashed black lines represent the upper and lower uncertainty envelopes of the firn air reconstructions. The magenta lines represent 50 scenarios generated randomly within the reconstructed uncertainty range, as described in Sect. 2.6. The middle panels show the temporal changes in the isotope signatures of the total N 2 O source, with their accompanied uncertainties, as calculated from the atmospheric mass balance model (Sect. 2.5). The total source is split into an assumed constant "natural" and an increasing "anthropogenic" component and the right panels show the isotopic evolution of the "anthropogenic" component. Results show that the average δ 15 N av of the total N 2 O source, over the reconstruction period, is (−7.6 ± 0.6) ‰ where the uncertainty is calculated using the 1σ uncertainty from the scenarios with respect to the mean value (magenta lines). There is no statistically significant long-term trend, but a temporal variability is observed on the decadal scale that might mask this trend. δ 15 N av first decreased from (−6.5 ± 0.6) ‰ in 1940 to (−8.5 ± 0.6) ‰ in 1965, then slowly increased again to (−6.6 ± 0.6) ‰ in 1985, followed by another decrease to (−8.5 ± 0.6) ‰ in 2008. These oscillations originate from the slightly curved trends in the isotopic reconstructions for δ 15 N av in Fig. 4 (left panels).
When the source is split into constant natural and varying anthropogenic components, the variability is projected onto the anthropogenic part and the temporal variations increase accordingly. However, the uncertainties also increase substantially, because the differences between the individual scenarios are attributed to only a small fraction of the total source.
The δ 15 N av signature of the anthropogenic source has an average value of (−18.2 ± 2.6) ‰. It initially increases (the small initial decrease is not significant) from (−21.5 ± 2.6) ‰ in 1940 to (−8.6 ± 2.6) ‰ in 1990, when it starts to slowly decrease reaching (−15.4 ± 2.6) ‰ in 2008. During the early part of the reconstruction period before 1970, when the "anthropogenic" contribution was only a small fraction of the total source, the uncertainty ranges of the source signatures are larger. Therefore, the uncertainties for the early part were excluded when calculating the 1σ uncertainties over the entire period from the generated scenarios. This applies to all anthropogenic isotope signatures. The budget calculations suggest an overall trend towards more enriched anthropogenic emissions, but the uncertainties are large. Mathematically, this trend arises from the fact that the isotope reconstructions yield relatively linear temporal isotope trends, whereas the source strength increases in a strongly non-linear fashion (Fig. 4). In the beginning of the record a small increase in the source strength needs to produce a certain absolute isotope shift, whereas a smaller increase in the source strength is needed during later years to cause a similar isotope shift. This can only be solved mathematically by a lower δ 15 N av value for the small "anthropogenic" emissions in the early part of the firn record. A constant δ 15 N av source signature would result in a small temporal change in δ 15 N av of atmospheric N 2 O in the beginning of the record and increasing isotope trends with increasing emissions, similar to the exponential curves that were fitted to the firn air data in Röckmann et al. (2003b).
The δ 18 O of the total source varies within (27.2 ± 2.6) ‰ over the entire period. δ 18 O does not show significant decadal-scale oscillations because the reconstructed scenario for δ 18 O is even more linear than the δ 15 N av scenario. For this reason, as explained above, in the best-fit scenario the δ 18 O of the anthropogenic source for the initial 30 years has a more depleted value starting with (7.7 ± 2.6) ‰ in year 1940, reaching (31.1 ± 2.6) ‰ in year 1975, and remaining around this value until 2008 (Fig. 4). However, the relatively larger uncertainty envelopes for the atmospheric history of δ 18 O actually allow scenarios with smaller δ 18 O changes in the beginning of the record and larger changes in the later period, which means that the reconstruction does not exclude a constant value for the anthropogenic δ 18 O source signature. The available dataset thus does not allow quantification of a longterm trend in δ 18 O.
For the position-dependent 15 N signatures of the total source, no significant long-term trends were detected. For δ 15 N α no decadal-scale variability is observed, whereas for δ 15 N β a temporal variability is observed similar to δ 15 N av . The uncertainty ranges for δ 15 N α and δ 15 N β are about a factor 2 greater than for δ 15 N av , which is due to the larger analytical error that leads to higher uncertainties in the scenario reconstructions. δ 15 N α varies in the range (−3.0 ± 1.9) ‰, δ 15 N β in the range (−11.7 ± 2.3) ‰.

Discussion
The N 2 O mole fraction atmospheric history from our multisite reconstruction is in agreement with recent work by Meinshausen et al. (2016), who combined all available published N 2 O data (atmospheric, firn, ice) in order to reconstruct a historical atmospheric record of the past 2000 years. It differs slightly from the one determined by Battle et al. (1996) and to a smaller extent by Machida et al. (1995). Battle et al. (1996) collected firn air data and Machida et al. (1995) used ice data. Both studies used samples from a single Antarctic site. One could argue that the difference is due to an interhemispheric difference, but it is too large to be explained by this alone. In the past, N 2 O mole fraction measurements have been reported on different calibration scales, which is likely to explain part of the differences between individual studies. Furthermore, differences in the firn air model and possible differences between sites may contribute. In our case we used measurements from five sites to constrain our model, while Battle et al. (1996) and Machida et al. (1995) used only one site. In addition, the atmospheric histories of up to nine known gases (depending on site, Witrant et al., 2012) were used to constrain diffusivity in our model, while Battle et al. (1996) only used two gases.
From the combination of the firn air reconstruction with a simple two-box model we conclude that N 2 O emissions increased from (11.9 ± 1.7) Tg year −1 N in 1940 to (16.4 ± 1.7) Tg year −1 N in 2008. This agrees, within uncertainties, with previous firn reconstruction studies from Ishijima et al. (2007) and Park et al. (2012) and bottom-up approaches using emission databases (Syakila and Kroeze, 2011;Kroeze et al., 1999). A more recent study by Thompson et al. (2014) performed inversions of atmospheric measurements for 2006 to 2008 with multiple models and reported emissions of 16.1-18.7 Tg year −1 N for 2008, which is also in agreement with our findings.
To investigate the effect the N 2 O lifetime on the N 2 O isotopic signatures (Prather et al., 2015, we performed a sensitivity study where we linearly changed the N 2 O lifetime from 123 years pre-industrially (≈ 1750) to 119 years in modern times (2008). The results are shown in Appendix D, where the effect on the emission strength and isotopic composition is discussed in detail. Results from this sensitivity study showed that the effect of a decreasing lifetime gives higher N 2 O emissions for year 2008 while keeping the same preindustrial value, confirming the sensitivity to the lifetime in line with Prather et al. (2015). This change in lifetime in the model leads to changes in the isotope signatures of the order of (2.0 ± 1.0) ‰. The lifetime effect is most pronounced for the earliest part of the record (< 1970) where the reconstruction uncertainties are larger than this systematic uncertainty.
We furthermore investigated the sensitivity to the value of F exch (stratosphere-troposphere flux) between a low and high value of 0.16 and 0.28 Tmol s −1 , respectively, following Appenzeller et al. (1996) and Holton et al. (1990), with the default value being 0.22 Tmol s −1 . As shown in Appendix D, the isotope values are not very sensitive to the changes in F exch ; the results are well within the uncertainty envelopes.
The increase in N 2 O emissions over the past decades resulted in an overall decrease in all isotopic signatures of atmospheric N 2 O with time. The isotopic signature of the total source of N 2 O (Fig. 4, middle panels) is strongly depleted in all heavy isotopes compared to tropospheric N 2 O (Table 3), which is due to the strong enrichment associated with the removal in the stratosphere. In Table 3 the isotopic composition for the pre-industrial period (≈ 1750) (δ nat,pi ) is compared with the derived anthropogenic source signature derived from our multi-site reconstruction (δ anth , averaged from 1940 to 2008). The results show that the anthropogenic source is more depleted in heavy isotopes than the natural one for all signatures, confirming results from studies prior to firn air measurements (Rahn and Wahlen, 2000), and from studies that used forward firn air modeling on measurements from individual sites Ishijima et al., 2007;Röckmann et al., 2003b). It is important to remember that we assume the natural sources to be constant, but the method itself does not provide evidence of this.
Anthropogenic N 2 O emissions are dominated by agricultural soil (70 %) with smaller contributions from automobiles, coal combustion, biomass burning, and industry. Oceanic emissions were previously assumed to be only natural. However, the latest IPCC Assessment Report (Ciais et al., 2013) for the first time separated oceanic emissions into a natural component and an anthropogenic component, e.g., due to atmospheric N deposition to rivers (Syakila and Kroeze, 2011;Duce et al., 2008;Kroeze et al., 2005). The oceanic fraction of the anthropogenic source was estimated as 1 Tg year −1 N. N 2 O emitted from agricultural soils and biomass burning is more depleted in δ 15 N av and δ 18 O than the tropospheric background (Park et al., 2011;Goldberg et al., 2010;Ostrom et al., 2010;Tilsner et al., 2003;Pérez et al., 2001Pérez et al., , 2000, while N 2 O emitted from other minor sources, such as automobiles, coal combustion, and industry, has values closer to tropospheric N 2 O values (Syakila and Kroeze, 2011;Toyoda et al., 2008;Ogawa and Yoshida, 2005a, b). An increase in strongly depleted agricultural emissions in the first part of our reconstruction, followed by a decreasing relative contribution from agriculture and increasing contributions from more enriched sources like industry, automobiles, and coal combustion, could qualitatively explain the reconstructed changes in isotope signatures of both the total source and the anthropogenic component. The global N 2 O budget study from Syakila and Kroeze (2011) indicates that agricultural emissions were 78 % of the total during the 1940-1980 period, with little input from industry, vehicle exhaust, and coal combustion. After 1980 the relative share of agricultural emissions dropped to 64 %, while the other sources increased, supporting our suggestion.
According to FAO statistics (http://www.fao.org/faostat/ en/#data/GY/visualize), emissions from synthetic nitrogenous fertilizers increased between 1961 and 1985, then stayed relatively constant or even decreased until 2000, and increased again after 2000. The reasons of the decrease between 1985 and 2000 are a shift towards organic soil cultivation in combination with more efficient agricultural methods and fertilizer use. This variation in fertilizer use qualitatively matches with the temporal evolutions of our reconstructed source signatures, but the trends in the reconstructions are likely too large to be explained by this source change only.
Although the decadal variability for δ 15 N av and δ 15 N β appears statistically significant with respect to the choice of scenarios constructed within the error bars of the firn air reconstruction, additional systematic uncertainties in this reconstruction could potentially produce such trends artificially from small undulations on the scenarios, since the emissions are related to the derivative of the trend. As it is possible to draw straight lines within uncertainty envelopes of the scenarios, the decadal variability may not be robust. An increase in the regularization term by 10 confirms that the generated scenarios are straight lines well within the uncertainty envelopes; thus, the decadal variability could be an artifact of the model (see Appendix C).
Additional evidence of potential changes in the N 2 O source composition between the pre-industrial and present atmosphere may be derived from the position-dependent 15 N signatures, quantified by the 15 N site preference. Table 3 shows that the difference in the δ 15 N av signature between the pre-industrial and anthropogenic sources derived from our reconstruction is primarily due to a change at position δ 15 N β , whereas δ 15 N α remains relatively constant. This is reflected by a larger difference in δ 15 N sp between natural and anthropogenic emissions, which could indicate a temporal change in production processes. Sutka et al. (2006) suggested that there may be two distinct classes of N 2 O sources with different δ 15 N sp . N 2 O produced during nitrification and fungal denitrification had a high δ 15 N sp of (33 ± 5) ‰ and N 2 O from denitrification and nitrifier denitrification had a low δ 15 N sp of (0 ± 5) ‰. Park et al. (2012) used these two endmembers to calculate a change in the relative fractions of these source classes over time based on their firn air data. Although this approach is strongly simplified and several other sources and factors may contribute (Toyoda et al., 2015), we use the results from our box model calculations (Table 3) in a similar way to estimate the fraction of the two source categories according to the following simple mass balance calculation: This returns a fractional contribution of the δ 15 N sp high component of (19 ± 4) % to the total pre-industrial emissions and (35 ± 11) % to the total present source. The errors were derived by propagating the errors of the δ 15 N sp endmembers and δ 15 N sp meas within the ranges stated above. We note that the errors associated with the precise isotopic composition of the endmembers are correlated if the values of δ 15 N sp for the two endmembers remain relatively constant in time. Therefore, the change in the relative fraction of the two categories is likely better constrained than the absolute values.
Splitting the total present emission strength into natural (pre-industrial, 11.0 Tg year −1 N) and anthropogenic (5.4 Tg year −1 N) components, we derive a fraction of the δ 15 N sp high component (which includes nitrification) of (54 ± 26) % for the "anthropogenic" emissions. This is another piece of evidence for agricultural sources being the main contributor to the N 2 O increase, because nitrificationdominated agricultural emissions can be associated with the δ 15 N sp high component. The temporal changes in the derived fraction of nitrification are in good qualitative agreement with the results from Park et al. (2012), who reported a change of (13 ± 5) % from 1750 to (23 ± 13) % today. However, the absolute numbers derived from our study are higher than the results from Park et al. (2012). The difference is due to the fact that different apparent isotope fractionations during stratospheric removal (ε app ) are used in the mass balance model (Table 3; Eqs. 7 and 8). In our study we used the averaged lowermost stratospheric apparent isotope fractionations from , which we consider more representative than the numbers used by Park et al. (2012). Using different values for ε app causes a shift in the isotopic source signatures from the mass balance model. The choice of this value thus adds a systematic source of uncertainty to the absolute value of the δ 15 N sp high fractions reported above (F high ).
Nevertheless, this systematic uncertainty should not alter the overall change in F high from pre-industrial to modern times and the results from our multi-site reconstruction of the isotopic composition of N 2 O thus confirm the suggestion by Park et al. (2012) that the relative importance of the high-SP component (presumably nitrification) has increased with increasing mole fraction since pre-industrial times.

Conclusions
The temporal evolution of the total N 2 O emission fluxes and the source isotopic composition have been estimated in a top-down approach using a multi-site reconstruction of N 2 O mole fraction and isotopic composition from 6 firn air samplings at 5 different Arctic and Antarctic locations in a twobox model. The results from a mass balance model constrain the source strength and suggest a total increase in N 2 O emissions of (4.5 ± 1.7) Tg year −1 N between the 1940 and 2008 due to anthropogenic processes. This agrees with previous top-down estimates, but deviates from bottom-up model estimates, which suggest higher N 2 O emission increases. A significant source of the uncertainty in top-down estimates is a possible change in the N 2 O lifetime over the reconstruction period, which we have quantified following the recent results from Prather et al. (2015).
The reconstruction of mole fraction and isotopic composition was used to investigate temporal changes in the isotopic signature of N 2 O emissions over the study period. The average total source for δ 15 N av and δ 15 N β shows no statistically significant long-term trend but possibly significant decadal scale variability. For δ 18 O and δ 15 N α of the total N 2 O source, no significant temporal changes can be detected with the present dataset because the uncertainties are large, especially in the beginning of the reconstruction period.
When the total source is split into a constant natural component and a varying anthropogenic component, the reconstruction of the δ values of the anthropogenic source indicates a significant increase in δ 15 N av from the early part to the modern part of the record. This originates from the near-linear isotope histories of the best-guess scenario, which would imply that small emissions in the early part had a similar absolute effect on the δ values to stronger emissions in the latter part. A similar effect for δ 18 O is likely but not significant given the larger uncertainties for this signature.
Nevertheless, the isotope signal in δ 15 N av may also be a signal for changing source contributions over time. Bottomup models suggest that N 2 O emitted from agricultural soils was the dominant contributor to the anthropogenic N 2 O increase in the first decades. Smaller contributions due to emissions from more enriched sources, like industry, automobiles, and coal combustion, increased. This may have contributed to an isotope enrichment of the emissions, which is not detectable within the error bars for the other isotope signatures. However, one has to be cautious with a firm interpretation of these trends since the reconstruction method itself may also induce decadal variability if the smoothness of the scenario is incorrectly constrained.
Results from the mass balance model yield an increase in 15 N site preference between the pre-industrial and modern total N 2 O sources. When this trend is evaluated with a simplified two-endmember mixing model, the results suggest an increase in nitrification sources relative to denitrificationrelated sources over the industrial period.
Data availability. A supplementary data table is available with the new NEEM-EU-08, NEEM-09 N 2 O mole fraction and isotopic signature firn measurements. Results from the N 2 O mole fraction and isotopic signature atmospheric reconstruction are also included. The respective raw data used in this study are available from the original referenced data providers upon request.

Appendix B: Data processing
In this study isotope deltas (δ) are used to denote the relative 15 N / 14 N and 18 O / 16 O ratio difference of N 2 O in firn air with respect to a standard reference,  Kaiser et al. (2008), assuming a constant 17 O excess of 0.9 ‰.
There is a disagreement between reported trends of the position-dependent δ 15 N av values reported in the literature from firn air on the one hand and archived air samples on the other hand Ishijima et al., 2007;Bernard et al., 2006;Röckmann and Levin, 2005;Röckmann et al., 2003b;Sowers et al., 2002). In principle the temporal trend measured directly on archived air samples should be fully consistent with top firn air samples of the various datasets, which were collected over a decade or more, since the air in the diffusive zone is not very old. However, this is not the case. Using the high-precision determination of the temporal trend of the N 2 O isotope signatures on archived air samples from Röckmann and Levin (2005) as reported in Sect. 2.4, we rescale the different firn profiles to match this trend in the diffusive zone by interpolating the measurements from the diffusive zone of all sites to DC-99 (δ INT ). By using the firn model-assigned mean age of each sample, the maximum age difference from diffusive zone to surface corresponds to age = DCt−t 0 = 10 years. Below you can find the equations used: where m is the slope connecting the two points we want to interpolate. The applied scaling (δ Final ) is given in   Figure C2. Results of the firn data evaluation (similar to Fig. 2) using the data re-scaled to the NEEM-09 site. Colors as in Fig. C1. Figure C3. Sensitivity test to the regularization term increased by a factor of 10. Reconstructed atmospheric scenarios (left), corresponding fit of the firn data (center), and model data discrepancies (right). The best reconstructed scenarios are shown as the black continuous lines, with model-derived uncertainties (2σ ) in dashed lines. Colors as in Fig. C1. Figure C4. Comparison of the atmospheric reconstructions between different re-scaling methods. Solid and dashed green lines are the scenarios from data re-scaled to DC-99 used in this study. Solid red lines are the best-case scenario for the non-re-scaled data and solid blue lines are the best-case scenarios from the data re-scaled to NEEM-09. The latter data series is shifted because of a calibration offset. When this is corrected for, the data superimpose the green lines, as expected. For the default calculations with the mass balance model a constant lifetime for N 2 O was used. A recent study from Prather et al. (2015), though, highlighted that top-down model calculations are sensitive to changes in the N 2 O lifetime. To quantify the effect on our results we performed a sensitivity test where we linearly changed the N 2 O lifetime from pre-industrial to modern times from 123 years in 1700 to 119 years in 2008. We also included runs with the absolute mean value changes in the assumed mean lifetime. The results are shown in Figs. D1 and D2 above.
In Fig. D1 the N 2 O atmospheric budget is re-calculated and compared with the results when the constant lifetime of 123 years is used. In year 1940 the N 2 O emission is (12.3 ± 2.7) and (17.0 ± 1.7) Tg year −1 N in year 2008, with a total increase of (4.7 ± 1.7) Tg year −1 N. When keeping the lifetime constant, the results for the same years are (11.9 ± 1.7) and (16.4 ± 1.7) Tg year −1 N with a total increase of (4.5 ± 1.7) Tg year −1 N. In addition, when also looking into the absolute mean value changes in the assumed mean lifetime, we only observe a vertical shift of the scenar- ios that do not affect the temporal change. This shows that there is a sensitivity to the choice of lifetime for our mass balance model on the N 2 O atmospheric budget, as was indicated by Prather et al. (2015). The N 2 O source isotopic signature shows no significant change, with the choice of lifetime giving similar average source values for all source signatures to when using a constant lifetime of 123 years.
Sensitivity tests were also performed on the F exch parameter which gives us the annual fluxes between the two reservoirs (stratosphere-troposphere). Following Appenzeller et al. (1996) and Holton et al. (1990) the value was tested at a low and high value of 0.16 and 0.28 Tmol s −1 , respectively, with the one used in the paper being 0.22 Tmol s −1 . Results are shown in Figs. D3 and D4 above. In Fig. D3 (middle panel) the atmospheric budget is recalculated and compared to the optimal scenario values. In the bottom panel the air returned to the troposphere from the stratosphere is presented (F exch ). It is clear that when a low F exch value is chosen, then less N 2 O is returned to the troposphere. Contrarily, when a higher F exch value is used, more N 2 O is returned.  Figure D4. Left panels: historic evolution of δ 15 N av , δ 18 O, δ 15 N α , and δ 15 N β in N 2 O as derived from the firn air reconstruction. The solid black line represents the best-fit scenario, while the dashed ones represent the respective uncertainties as determined by the reconstruction method. Magenta lines represent the emissions that are required to produce the magenta N 2 O histories in the left panels. Middle and right panels: isotope signatures of the total emitted N 2 O and anthropogenic source, respectively, assuming high (0.28 Tmol s −1 ) F exch in light green and a low (0.16 Tmol s −1 ) value in yellow. The solid black line represents the result for the best-fit reconstruction, while magenta lines represent the results for the individual scenarios from the top as used in the main paper.
F exch choice has little effect on the isotopic signature results as shown in Fig. D4 and is mainly limited to the earliest part of the record (> 1970) where the reconstruction uncertainties are larger. While it is expected that when the F exch value is low, the isotopic results will be more enriched compared to higher F exch , in our case this is not clear from the test. The overall averaged values have a less than 2 % difference compared to the chosen (optimal) scenario and results of total averaged source and anthropogenic isotopic signatures are well within agreement with combined uncertainty errors in both total source and anthropogenic signatures, respectively.
Thus, we conclude that while the flux is indeed sensitive to the F exch choice value, the isotopic composition is not.