Atmospheric Nitrous Oxide Variations on Centennial Time Scales During the Past Two Millennia

The continuous growth of atmospheric nitrous oxide (N2O) is of concern for its potential role in global warming and future stratospheric ozone destruction. Climate feedbacks that enhance N2O emissions in response to global warming are not well understood, and past records of N2O from ice cores are not sufficiently well resolved to examine the underlying climate‐N2O feedbacks on societally relevant time scales. Here, we present a new high‐resolution and high‐precision N2O reconstruction obtained from the Greenland NEEM (North Greenland Eemian Ice Drilling) and the Antarctic Styx Glacier ice cores. Covering the N2O history of the past two millennia, our reconstruction shows a centennial‐scale variability of ~10 ppb. A pronounced minimum at ~600 CE coincides with the reorganizations of tropical hydroclimate and ocean productivity changes. Comparisons with proxy records suggest association of centennial‐ to millennial‐scale variations in N2O with changes in tropical and subtropical land hydrology and marine productivity.


Introduction
Improved knowledge of greenhouse gas-climate feedbacks is required to understand past and future climate changes. Nitrous oxide (N 2 O) is a particularly important greenhouse gas with a global warming potential 260 times greater than that of CO 2 for a time horizon of 100 yr (Myhre et al., 2013). With the regulation of chlorofluorocarbon emissions, N 2 O is becoming the most important ozone-destroying substance in the stratosphere (Ravishankara et al., 2009). Nonetheless, the processes governing atmospheric N 2 O variability still remain elusive.
Microbial production in soils and the ocean are dominant sources of atmospheric N 2 O (Butterbach-Bahl et al., 2013). The atmospheric N 2 O concentration is primarily regulated by the balance of these terrestrial and oceanic sources (Freing et al., 2012) with the tropical stratospheric upwelling and photolysis in the stratosphere (Khosrawi et al., 2013;Olsen et al., 2001;Prather et al., 2015). The use of nitrogen fertilizer and other industrial activities in the modern era have added extra reactive nitrogen to the global ecosystems, causing a continuous increase in atmospheric N 2 O levels since the industrial revolution (Galloway et al., 2008). In addition, natural emissions are thought to be sensitive to future climate change (Battaglia & Joos, 2018;Denman et al., 2007;Martinez-Rey et al., 2015;Voigt et al., 2017).
Nitrification and denitrification are considered to be the two major pathways for natural N 2 O production (Butterbach-Bahl et al., 2013; Ji et al., 2015). In the presence of molecular oxygen, ammonium (NH 4 + ) can be oxidized to nitrite (NO 2 − ) and nitrate (NO 3 − ) (nitrification), with N 2 O as a by-product. In anoxic or suboxic conditions, denitrification reduces NO 3 − or NO 2 − into dinitrogen gas (N 2 ) in a stepwise fashion, with N 2 O as an intermediate product. Because these microbial processes are highly dependent on temperature and oxygen levels in both marine and terrestrial environments, the feedbacks between N 2 O production and climate conditions are important.
Our understanding of climate-N 2 O feedbacks relies, in part, on paleoatmospheric records and modeling studies. In this context, air bubbles trapped in polar ice cores provide a unique archive for reconstructing and testing ancient atmospheric composition changes through climate history. Several detailed ice core N 2 O records have been developed over the past few decades (Fischer et al., 2019;Flückiger et al., 2002Flückiger et al., , 2004Schilt et al., 2010;Spahni et al., 2005). These records capture long-term variations, such as glacial-interglacial cycles, millennial-scale variations during the Holocene, and the N 2 O responses to specific abrupt climatic events, including the preboreal transition at the end of the last glacial period and earlier Dansgaard-Oeschger events (Fischer et al., 2019;Flückiger et al., 1999Flückiger et al., , 2002Schilt et al., 2010Schilt et al., , 2014Sowers et al., 2003;Spahni et al., 2005). Although previous studies reported pronounced covariability of N 2 O with northern hemispheric temperature on glacial-interglacial and millennial scales (Flückiger et al., 2002;Schilt et al., 2010), available ice core N 2 O records for the Holocene have not been sufficiently consistent to allow for an examination of smaller changes on submillennial time scales. Specifically, the past two millennia are thought to be affected by natural climate variability and the growing human influence. Centennial-scale variations are not well resolved in existing ice core N 2 O records (MacFarling Meure et al., 2006;Prokopiou et al., 2018), apart from the strong atmospheric N 2 O increase over the past 200 yr (Flückiger et al., 1999;Machida et al., 1995;Sowers et al., 2002).
To improve our understanding of the key drivers of atmospheric N 2 O variability on centennial time scales, we present new high-resolution (~15 yr) ice core N 2 O records from both the Greenland NEEM (North Greenland Eemian Ice Drilling) and the Antarctic Styx Glacier ice cores, which cover the past two millennia, and investigate the underlying control mechanisms.

Samples and Gas Chronology
Styx Glacier ice and NEEM ice were used for N 2 O measurements. The climatic information of the ice coring sites is listed in supporting information Table S1. The 210.5 m long Styx Glacier ice core and firn air samples were obtained during the Korean ice core drilling campaign in 2014-2015 (Han et al., 2015). The gas age scale of the Styx Glacier record was obtained by synchronizing its CH 4 record to the West Antarctic Ice Sheet (WAIS) divide CH 4 record on the WD2014 age scale (Buizert et al., 2015). The uncertainty of the correlation is~±20 yr . The Δage (difference in ice age and gas age) of the Styx Glacier ice was estimated to be~318 yr . The gas chronology of the NEEM main ice core was also synchronized to the WD2014 age scale to solve the observed age discrepancy of~20 yr between the NEEM-2011-S1 and WD2014 age scales (Rhodes et al., 2013) ( Figure S1). The NEEM-2011-S1 core is a shallow parallel core to the NEEM main core. The estimated age difference between the main ice core and S1 ice core is less thañ 5 yr (Sigl et al., 2015), smaller than the age uncertainty resulted from wiggle matching CH 4 records (~±20 yr). No high-resolution CH 4 data exist below 408.96 m (291 CE) (Rhodes et al., 2013). We estimated the gas age for the deeper NEEM ice by linear extrapolation of the depth-gas age relationship at depth of 338-408 m.

Ice Core N 2 O Measurements
We used a high-precision measurement method  to analyze both the NEEM and Styx Glacier ice and generate a composite N 2 O record. Ice core samples are cut into subsamples (~20 g) in a walk-in freezer at a temperature of −20°C to prevent melting. For each depth, samples are duplicated or quadruplicated to estimate the reproducibility of measurements. Ice pieces are placed in glass flasks sealed to a Conflat Flange with a copper gasket for sealing at a high vacuum level. The flasks are then submerged into a prechilled ethanol bath (< −75°C) and evacuated to remove the ambient air in the flask and any contaminants on the ice core surface, for 50 min. After sufficient evacuation, gases in the ice core bubbles are liberated by melting. The N 2 O mixing ratio of the liberated air in the headspace is measured by an Agilent 7890B Gas Chromatograph equipped with a Micro-Electron Capture Detector. Before and after the sample measurement, we analyzed a standard N 2 O gas with a concentration of 329.9 ppb from the National Oceanic and Atmosphere Administration NOAA-2006A N 2 O scale. Due to the high solubility of N 2 O in water, we repeated the freezing-melting cycle to liberate the N 2 O trapped during the extraction procedure. The additional freezing-melting step results in a~6% correction to the total N 2 O concentration.

Data Quality
Both ice cores have fairly high accumulation rates (NEEM: 0.22 and Styx Glacier: 0.13 m ice equivalent/yr; supporting information Table S1), resulting in a small smoothing effect on gas records due to gas diffusion and gradual bubble close-off processes in the firn layer (transition zone from snow to ice on top of the ice sheet). The estimated widths of the gas age distribution at half height from firn densification models for both ice cores are smaller than 40 yr (Buizert et al., 2012;Jang et al., 2019). Thus, both ice cores can resolve centennial-scale changes in atmospheric composition. We estimated the uncertainty of the N 2 O concentration by replicating measurements with series of adjacent eight samples within~20 cm intervals (corresponding to 1-2 yr in mean age change) and obtained pooled standard deviations of 3.4 ppb for Styx Glacier ice and 2.8 ppb for NEEM ice (supporting information Table S2). These are greater than the analytical uncertainty of 1.5 ppb which is determined from measurements of Styx replicate ice samples within 5 cm depth intervals , indicating that N 2 O concentration in the ice core varies in cm scales. Alteration of N 2 O concentration by in situ microbial activity in the ice (Miteva et al., 2016;Rohde et al., 2008) could contribute the N 2 O variations and may explain the greater scatter in longer sample interval. However, we do not believe that in situ production of N 2 O can explain the centennial-scale variations we describe here, because the variations in N 2 O measured in the two cores we studied are remarkably coherent (Figure 1a) on centennial time scales. The two cores have different gas age-ice age (Δage) differences (~310 yr for Styx, Jang et al., 2019, and~188 yr for NEEM, Buizert et al., 2012), which means that dust variations would be offset differently with the gas record in each core. Instead, the agreement indicates that these centennial changes represent the true atmospheric signal. Owing to the higher anthropogenic emissions in the Northern Hemisphere, current atmospheric N 2 O levels are 0.8 ppb higher in the Northern Hemisphere than in the Southern Hemisphere (Ishijima et al., 2009). However, the uncertainty of our ice core records from Greenland and Antarctica does not allow to detect such a small N 2 O gradient for the preindustrial period. To focus on the joint variability in both cores, we develop a composite N 2 O record by averaging the NEEM and Styx Glacier N 2 O records, which were initially interpolated onto the same time axis. Our new data show subtle centennial-scale changes, which were not well resolved in previous ice core N 2 O records (Fischer et al., 2019;Flückiger et al., 1999Flückiger et al., , 2002MacFarling Meure et al., 2006;Schilt et al., 2010) (supporting information Figure S2).

Results
As shown in Figure 1, the composite N 2 O records exhibit two distinctive features during the late preindustrial Holocene: (1) a minimum near 600 CE and (2) centennial-scale fluctuations of about 10 ppb. The 600 CE minimum was previously observed in the Law Dome ice core N 2 O record, although the uncertainty of the N 2 O data (± 6.5 ppb) was not sufficient to resolve the magnitude of the N 2 O decrease (MacFarling Meure et al., 2006). In our measurements, N 2 O from both NEEM and Styx Glacier confirms that a local minimum resides near 600 CE, thereby consolidating the finding of a substantial change in the land and/or ocean N 2 O fluxes (Figure 1a). The N 2 O concentration decreases from~265 ppb at 200 CE to~257 ppb at 600 CE.  (Schilt et al., 2014). Gray shading area denotes ±1σ of the mean. For (b) and (c), pale orange shading after 1800 CE denotes the period considered to be strongly affected by anthropogenic emission.

Global Biogeochemical Cycles
This~8 ppb decrease corresponds to average reduction of~0.3 Tg N yr −1 in N 2 O emissions. This value is calculated with a simple mass burden to mixing ratio (4.79 Tg N/ppb) (Prather et al., 2015) and the assumption that changes in the stratospheric N 2 O sink are negligible (see supporting information).
Centennial-scale N 2 O fluctuations are also prominent (Figure 1b). To further examine centennial-scale variations in the N 2 O budget, we estimated emissions from the concentration data using a two-box model (Figure 1c). To estimate N 2 O emissions, Monte Carlo simulations incorporating uncertainties in N 2 O concentration and lifetime were run for each time interval (see supporting information for details). We find that emission ranges from 10.2-11.7 Tg N yr −1 on average, similar to the current estimates of the natural N 2 O flux from both top-down and bottom-up approaches (Davidson & Kanter, 2014) (10-12 Tg N yr −1 ). As expected, the atmospheric N 2 O concentration shows a delayed response of several decades after the change in the N 2 O emission (Figures 1b and 1c).

Discussion: N 2 O and Paleoclimate
N 2 O production relies on the metabolic processes of specific groups of microorganisms. Both terrestrial and marine N 2 O sources are expected to respond to changes in the large-scale climate conditions. Environmental factors such as temperature and oxygen availability are believed to have significant impacts on the regulation of N 2 O production. In the case of terrestrial N 2 O sources, which account for~70% of the global N 2 O emissions, the water-filled pore space (WFPS) is one of the key factors controlling oxygen availability and N 2 O production in the soil environments (Davidson et al., 2000). An optimum WFPS of 70-80% in the soil increases microbial activity and N 2 O emissions (Bouwman et al., 2013). Tropical forest soils are the main terrestrial N 2 O source (Stehfest & Bouwman, 2006;Xu et al., 2017), and the rainfall variations throughout the monsoon season have been suggested as the key factor influencing the interannual variability of the N 2 O emission (Werner et al., 2007). The soil temperature is another important climate factor that can enhance terrestrial microbial N 2 O production. Model studies suggest that on decadal to multidecadal time scales, temperature is the main climate parameter driving changes in terrestrial N 2 O emission (Xu et al., 2012;Zaehle et al., 2011). Oceanic N 2 O production is tightly associated with high primary productivity regions. Ocean systems are considered to contribute approximately 30% of natural N 2 O emissions (Voss et al., 2013), and high productivity regions contribute a greater portion of the marine N 2 O emissions. The eastern tropical Pacific (ETP) and Arabian Sea (AS), where considerable N 2 O accumulation occurs on the upper boundary of widely developed oxygen minimum zones (OMZs) (Hamersley et al., 2007;Ito & Deutsch, 2013), are well-known major oceanic N 2 O source regions. The thermocline depth and strength of the upwelling activity in these regions may influence the nutrient supply that sustains the primary productivity and, in turn, the N 2 O production rate. The dissolved O 2 in upper ocean can affect N 2 O production and emission variability since both denitrification and nitrification are sensitive to dissolved O 2 , closely tied to aerobic remineralization of organic matter (Battaglia & Joos, 2018;Martinez-Rey et al., 2015).
The local N 2 O minimum at 600 CE coincides with an overall weakening of the tropical and subtropical monsoon and lower productivity in major oceanic N 2 O source regions, as indicated by a variety of palaeoclimate proxy records (Figure 2 and supporting information Figure S3). The N 2 O level continuously decreased until 600 CE and then increased until 1200 CE. The timing of this N 2 O minimum is coincident with reduced precipitation in the tropical land N 2 O source areas, including subtropical China (Wang et al., 2005), Southeast Asia (Steinke et al., 2014;Wurtzel et al., 2018), India (Anderson et al., 2002(Anderson et al., , 2010Gupta et al., 2003), and central America (Bhattacharya et al., 2015;Curtis et al., 1996). Because the optimum conditions for microbial N 2 O production requires a high level of soil moisture (~70% to 80%) to achieve low O 2 availability, a weakening in the monsoon strength in the tropical land source areas may decrease N 2 O production. The hydroclimate proxies in South America are not sufficiently consistent to confirm changes in the rainfall patterns at 600 CE; however, some of them suggest an abnormally drier climate (Bird et al., 2011;Kanner et al., 2013;Novello et al., 2012), illustrating a potential regional contribution of hydroclimate changes in South America to lower N 2 O emissions in this period. There are a few available hydroclimate proxy records in central Europe (Fohlmeister et al., 2012) and central Africa (Shanahan et al., 2009;Wolff et al., 2011), resolving the last 2,000 yr; however, the climate in these areas became wetter or showed no significant changes in this period, possibly contributing to the increase in the regional N 2 O emission. The marine sediment records from AS and the equatorial eastern Pacific, which account for a significant portion of the oceanic N 2 O emissions, show an overall reduction in productivity near 600 CE, which would further lead to a reduction in the N 2 O emission ( Figure 2).
To deduce the extent and patterns of the hydrological changes on land, we compared our N 2 O composite with selected temporally well-resolved hydroclimate proxies (Figure 3). Among the cave speleothem proxy records (δ 18 O of speleothem calcium carbonate) in China, representing the strength of the Eastern Asian Summer Monsoon (EASM), the lower-latitudinal cave speleothem records from Dongge cave (Wang et al., 2005) (Zhang et al., 2008) caves (33.32°N) records do not share this pattern (supporting information Figure S4). The Wanxiang cave records (Zhang et al., 2008) indicate an even stronger EASM between 400 and 800 CE. This latitudinal discrepancy in the EASM proxies is suggested to result from changes in the width of the tropical rainfall belt with the migration of the Intertropical Convergence Zones (Denniston et al., 2016). Thus, the weakening of the EASM strength was possibly constrained to the tropical and subtropical regions.
Like the tropical and subtropical western Pacific regions, whose hydrology is largely affected by the monsoon, the precipitation in the Amazon River basin is also greatly controlled by South American Summer Monsoon (SASM). Speleothem records from Huagapo cave (11.27°S, 75.79°W) in the Eastern Andes show a weakening of the SASM during the N 2 O minimum at 600 CE (Kanner et al., 2013) (Figure 3d); however, other SASM proxies are not sufficiently consistent to confirm its weakened strength at 600 CE ( Figure 2 and supporting information Figure S3). Although our interpretation of these hydroclimate proxies is limited by the discrepancies among the proxy records, a large portion of the tropical and subtropical monsoon regions appear to indicate a change toward drier conditions near 600 CE.
What caused these globally diminishing monsoon activities at 600 CE is unclear; however, hydroclimate change and the associated reduced soil moisture content in the tropical N 2 O source regions likely contributed to the reductions in the terrestrial N 2 O emissions. During this period, the northern hemispheric temperature is believed to be lower than average (Ljungqvist, 2010;Mann & Jones, 2003;Moberg et al., 2005) (supporting information Figure S5), which could further subdue the terrestrial N 2 O emission by suppressing microbial activity. Previous work suggested that northern hemispheric temperature change caused millennial-scale N 2 O changes in the Holocene, showing parallel variation with N 2 O (Flückiger et al., 2002;Schilt et al., 2010). However, a quantitative estimation of the temperature change is difficult owing to the spatially dispersed temperature proxy records.
The weakening of the ISM prevents upwelling and reduces the primary productivity in the AS, likely leading to decreased N 2 O production owing to the reduction in organic matter input, which is required for denitrification. The abundance of the planktonic foraminifera Globigerina bulloides in a sediment core from AS (Anderson et al., 2002(Anderson et al., , 2010Gupta et al., 2003) (Figure 3e) shows a minimum at 600 CE, reflecting the reduction in the primary productivity, consistent with the reduction in the marine N 2 O flux from AS. Similarly, the nitrogen isotope ratios (δ 15 N) of organic matter in a sediment core off Peru suggest a reduction in the primary productivity and weakened OMZ development in ETP near 600 CE, followed by a gradual strengthening (Salvatteci et al., 2014) (Figure 3f). OMZs are the regions where active nitrogen loss (both denitrification and anammox) take place, contributing~4% of the estimated oceanic N 2 O production (Battaglia & Joos, 2018). In addition, the combination of aerobic remineralization of organic matter followed by nitrification results in strong regional N 2 O production where OMZs are well developed (Yang et al., 2020). Thus, a reduction in the OMZ in ETP and AS could also lead to a decrease in oceanic N 2 O emissions. In summary, both terrestrial and marine palaeoclimate records suggest that the reorganization of O records from Huagapo cave (P00-H1, P09-H2 cores), Peruvian Andes (Kanner et al., 2013). (e) G. bulloides content in Arabian Sea sediment cores (Anderson et al., 2010;Gupta et al., 2003) (filled circles: RC2730, RC2735 cores, Anderson et al., 2010; open circles: 723A core, Gupta et al., 2003). (f) δ 15 N of organic matter in off Peru sediment core (Salvatteci et al., 2014). (g) Combined tree ring-based 31 yr running variance of interannual-scale variability ENSO records (red line) with their ±1σ uncertainty (pale red shading) (Liu et al., 2017). N 2 O flux is plotted on same panel (gray line). (b-d) For clarity, highly resolved records are smoothed with a moving average technique with a time window of 50 yr to highlight multidecadal variations (thick solid lines). Pale blue shading represents the period of the lowest N 2 O level during the last two millennia. Correlation coefficients between each proxy and N 2 O concentration (b-f) or N 2 O flux and ENSO variance (g) are denoted.
tropical rainfall and wind patterns may have contributed to the marked decrease in N 2 O fluxes into the atmosphere near 600 CE.
This period was thought to be affected by changes in the periodical North Atlantic cooling events, known as Bond events (Bond et al., 2001). The North Atlantic cooling can be linked to the changes in Asian southeast monsoon and AS productivity (Gupta et al., 2003). Thus, the North Atlantic cooling can be associated with the subdued productivity of N 2 O near 600 CE. Dongge cave records may indicate a possible linkage between Asian monsoon strength and the North Atlantic cooling (Wang et al., 2005). Possible changes in climate during this time period include Pacific sea surface temperature (SST) gradient changes (Kanner et al., 2013) and northern hemispheric temperature changes (Kathayat et al., 2017), both of which could contribute to the reorganization of the tropical and subtropical hydroclimates.
On centennial time scales, regional monsoon proxies and oceanic productivity records are not highly correlated with our N 2 O flux estimates. The low correlation among the N 2 O flux and climate proxies might be because of age uncertainty of proxies and our N 2 O data, which can reach several decades in age difference between the N 2 O and the palaeoproxy data. Tropical rainfall change in such a short-term temporal scale may be spatially inconsistent and possibly lead to the lack of predominance in the N 2 O emission. The only meaningful correlation for the N 2 O flux was achieved with the annually resolved central Pacific El Niño-Southern Oscillation (ENSO) proxies ( Figure 3g). ENSO variance records based on Taiwan tree ring δ 18 O show in-phase correlation with the calculated N 2 O flux (Liu et al., 2017) (r ¼ 0.41). The centennial-scale decrease in the N 2 O flux corresponds to a suppressed ENSO variance, whereas the periods of high ENSO variance coincide with enhanced N 2 O emission. A correlation coefficient r of 0.41 is a conservative value, given the relative age uncertainty between the records. While we cannot address the physical mechanism between the central Pacific ENSO and N 2 O flux changes, it may reflect the relationship between the N 2 O production and climate variability in low-latitudinal regions.
The discussion above focuses on N 2 O production changes rather than changes in sinks related to photochemical decomposition in the stratosphere. The fluctuations in our N 2 O records appear to match with proxies of solar activity changes. N 2 O concentrations are in phase with total solar irradiance reconstruction data (Steinhilber et al., 2009), showing low levels of N 2 O during the solar minima events of the past millennium (supporting information Figure S6). A recent study reported that the N 2 O lifetime varies with solar activity (Prather et al., 2015); the lifetime of N 2 O decreases in solar maxima (while it increases in solar minima), indicating more rapid N 2 O loss and decrease of N 2 O during the solar maxima. However, our observation is a positive correlation between N 2 O and solar activity. This may imply that the changes in N 2 O loss due to solar activity are not large enough to counterbalance changes in N 2 O sources during the past two millennia.
To better understand the physical linkage between climate and N 2 O changes, further numerical modeling experiments are required as well as chemical and isotopic data. The isotopic composition of N 2 O molecule (δ 15 N, δ 18 O, and site preference) has been used as the footprint for distinguishing its sources and production processes (Fischer et al., 2019;Prokopiou et al., 2018;Schilt et al., 2014;Sutka et al., 2006), and it may provide a way to deduce the mechanistic explanation for the N 2 O variations. However, published N 2 O isotope data do not have sufficient temporal resolution and precision to facilitate the investigation of the relative contribution of marine and terrestrial sources or microbial production mechanisms (nitrification vs. denitrification) on submillenial time scales (Prokopiou et al., 2018). Future studies of N 2 O isotopologues and higher-resolution N 2 O isotope records will be very useful to improve our understanding of the connection between N 2 O emissions and regional and large-scale climate.

Conclusion
We present the submillennial-scale N 2 O variations during the past two millennia with highly resolved ice core records and evaluate their possible causes. Our data show variability of~10 ppb on a centennial time scale. A pronounced local minimum at~600 CE coincides with the changes in tropical hydroclimate and ocean productivity. Both terrestrial and marine N 2 O production rate changes in response to climate are likely explanations for the atmospheric N 2 O concentration changes. Disentangling the underlying changes in terrestrial and marine N 2 O production rate changes requires further analysis of N 2 O isotopes and targeted modeling experiments with state-of-the-art Earth system models.