Marine N2O emissions during a Younger Dryas-like event: the role of meridional overturning, tropical thermocline ventilation, and biological productivity

Past variations in atmospheric nitrous oxide (N2O) allow important insight into abrupt climate events. Here, we investigate marine N2O emissions by forcing the Bern3D Earth System Model of Intermediate Complexity with freshwater into the North Atlantic. The model simulates a decrease in marine N2O emissions of about 0.8 TgN yr−1 followed by a recovery, in reasonable agreement regarding timing and magnitude with isotope-based reconstructions of marine emissions for the Younger Dryas Northern Hemisphere cold event. In the model the freshwater forcing causes a transient near-collapse of the Atlantic Meridional Overturning Circulation (AMOC) leading to a fast adjustment in thermocline ventilation and an increase in O2 in tropical eastern boundary systems and in the tropical Indian Ocean. In turn, net production by nitrification and denitrification and N2O emissions decrease in these regions. The decrease in organic matter export, mainly in the North Atlantic where ventilation and nutrient supply is suppressed, explains the remaining emission reduction. Modeled global marine N2O production and emission changes are delayed, initially by up to 300 years, relative to the AMOC decrease, but by less than 50 years at peak decline. The N2O perturbation is recovering only slowly and the lag between the recovery in AMOC and the recovery in N2O emissions and atmospheric concentrations exceeds 400 years. Thus, our results suggest a century-scale lag between ocean circulation and marine N2O emissions, and a tight coupling between changes in AMOC and tropical thermocline ventilation.


Introduction
Nitrous oxide (N 2 O) is an atmospheric trace gas that contributes to the greenhouse effect and to the destruction of stratospheric ozone. It is emitted to the atmosphere from natural processes in terrestrial and marine systems. Ice core reconstructions have uncovered natural variations of atmospheric N 2 O between ∼195 and 290ppb over the past 800 000years, synchronous with variations in climate (Flückiger et al 1999, Sowers et al 2003, Spahni et al 2005, Schilt et al 2010a, 2010b, 2013. Of special interest are N 2 O variations during rapid climate events such as the Dansgaard-Oeschger (DO) events (Dansgaard et al 1984, Stocker and Johnsen 2003, Schilt et al 2010b, Kindler et al 2014 or the Younger Dryas (YD, 12.8 to 11.7 kyr BP) as their analyses may provide insight on climatic mechanisms and their timing. For example, the YD Northern Hemisphere (NH) cold swing is associated with a large-scale reorganization of ocean circulation, thought to have been provoked by freshwater release from ice sheet melting, shifts in the Inter Tropical Convergence Zone (ITCZ) and precipitation patterns, and variations in atmospheric N 2 O on the order of 30 ppb (Flückiger et al 2004, Schilt et al 2014. Only recently, attribution of emission changes to terrestrial and marine sources has become possible from ice core reconstructions of N 2 O's isotopic composition and the fact that marine and terrestrial N 2 O emissions have distinct isotopic signatures (δ 18 O and δ 15 N) (Sowers et al 2003, Schilt et al 2014. Since the reconstruction by Schilt et al (2014), new measurements of N 2 O and its isotopic composition have emerged and reconstructions of N 2 O emission changes from terrestrial and marine sources over the past 21 000 years are now available .
Here we use the improved isotope inversion by Fischer et al (2019) to constrain marine N 2 O emissions changes over the YD interval. The data shows that terrestrial and marine sources both increased between the Last Glacial Maximum and preindustrial (PI) by 1.7 and 0.7 Tg N yr −1 , respectively . Marine sources show intermittent drops into the Heinrich Stadial 1 (HS1) and into the YD on the order of 0.5 Tg N yr −1 (Schilt et al 2014, which within dating uncertainties occur in parallel to reductions in the Atlantic Meridional Overturning Circulation (AMOC) as recorded in marine sediment records (McManus et al 2004). Although the interpretation of the relative phasing of AMOC changes and marine emission reconstructions using N 2 O isotopes is hampered by independent age scales and the resolution in both marine sediment and ice core records , the timing and strength of reconstructed emissions provide important benchmarks for model evaluation and process understanding.
Only few paleo modeling studies have investigated marine N 2 O emissions, with conflicting results. Goldstein et al (2003) propose that marine and terrestrial emissions both contributed about equally to atmospheric N 2 O variations, whereas Schmittner and Galbraith (2008) suggest that marine changes were the dominant mechanism driving atmospheric N 2 O variations during past rapid climate events. The new ice core based emission reconstructions (Schilt et al 2014 resolve this dispute for the last termination, showing that both terrestrial and marine emissions are involved in the glacial/ interglacial as well as in the stadial/interstadial N 2 O change, and allow more direct model evaluation based on emission changes. In light of the new emission reconstruction, we re-investigate past changes in marine N 2 O emissions at the example of the YD within the Bern3D Earth System Model of Intermediate Complexity, while deglacial changes in terrestrial N 2 O emissions are addressed elsewhere . The goal of this study is to gain insight into mechanisms and the relative timing of changes in the AMOC and marine N 2 O emissions and atmospheric N 2 O in response to YD-type marine reorganizations. To this end, we test the model sensitivity of marine N 2 O emissions to a freshwater pulse in the North Atlantic under 14 kyr BP boundary conditions. We investigate how ocean circulation changes affect different N 2 O production pathways through changes in remineralization and O 2 concentrations. We employ a new, observation-calibrated N 2 O parameterization that includes N 2 O production as an O 2 -dependent byproduct from nitrification and as an obligate intermediate product of denitrification (Battaglia and Joos 2018a). Given the interactions between the marine carbon cycle, N 2 O, oxygen, and climate, understanding of changes in atmospheric N 2 O offers an additional constraint on past climate events.

Marine N 2 O parameterization and analysis
Marine N 2 O production, consumption, and emissions follow the parameterization developed in Battaglia and Joos (2018a, section 2.2) and the parameter values are chosen according to their best consistent parameter set (Battaglia and Joos 2018a, table 1). This novel parameterization explicitly considers, in contrast to other global parameterizations, N 2 O production and consumption by denitrification (Babbin et al 2015). Parameter values were, for the first time, tuned to both water column and surface N 2 O data in a Bayesian framework yielding a confidence range (±1 SD) for modern N 2 O emissions from 3.0 to 6.1 TgN yr −1 , smaller than the published range (1.8-9.45 Tg N yr −1 ) (Ciais et al 2013). Modeled mean O 2 concentrations in each grid cell are taken to represent a distribution of O 2 concentrations and micro-environmental conditions at the subgrid scale, allowing aerobic and, comparatively slower, anaerobic remineralization of organic matter to co-occur within the same grid cell. Aerobic remineralization sets in at a mean O 2 of ∼4.4 mmol m −3 and anaerobic remineralization is phasing out at ∼7.4 mmol m −3 . N 2 O production is coupled to the aerobic/anaerobic remineralization flux below the euphotic zone with an O 2 -dependent yield for nitrification and with a stoichiometric factor for denitrification. N 2 O consumption from denitrification is modeled with first order kinetic, constrained by organic matter availability and stoichiometry.
For analysis, N 2 O produced from nitrification is assigned to two contributions. The 'α-term' scales linearly with aerobic remineralization fluxes below 75 m depth. The 'β-term' is also coupled to aerobic remineralization fluxes below 75 m depth and its yield decreases exponentially with increasing O 2 concentrations (Battaglia and Joos 2018a).
Ideal age is simulated as an additional tracer to estimate ventilation times in the model. The tracer is set to zero at the surface and increases by 1 unit time per unit time evolved in the ocean interior. Preformed PO 4 3and four separate N 2 O tracers are included as in Battaglia and Joos (2018a) to assist interpretation of results and to attribute N 2 O inventory changes. Three tracers track the inventories associated with nitrification (α-term and β-term) and denitrification and one tracer tracks the N 2 O solubility component.

Experimental design
An idealized YD-type experiment is conducted, starting from 14 kyr BP boundary conditions. For spin up, the model was initialized with an available spin up for PI conditions and re-equilibrated for 30 000 years under 14 kyr BP conditions (see supplementary for details available online at stacks.iop.org/ERL/14/ 075007/mmedia). A triangular-shaped freshwater flux, equivalent to 2.1 m global sea level rise, is released over 2000 years to the northern North Atlantic (figure 1(a)). Duration and magnitude of the release was chosen to approximately match the duration of the YD cold period (12.8-11.7 kyr BP; Rasmussen et al 2014) and to provoke a reduction and recovery in AMOC as evident in marine proxy data ( . This simulation and a 14 kyr BP control run are integrated for 10 000 model years. The duration of the freshwater forcing, and similarly of associated perturbations, is about 900 years longer than the observed duration of the YD. Note that we intentionally did not over-tune the freshwater forcing of the AMOC to optimally match the proxy records, as there is no unique solution to achieve a reasonable AMOC reduction and AMOC responses vary across models. In addition, climate boundary conditions are kept constant. Thus, model runs are not intended to exactly reproduce reconstructed marine N 2 O emissions and associated time scales, as these time scales are dependent on the assumed AMOC forcing in the model. However our model can provide an understanding of the overall phase relationship between AMOC shutdown and marine N 2 O emissions and on the absolute magnitude of N 2 O production changes in the ocean in response to a substantial reduction of the AMOC.

Global response to circulation changes
In response to the idealized freshwater forcing into the North Atlantic region (figure 1(a)), the AMOC slows down from 19.6 to 5.8 Sv over the course of ∼1000 years (figure 1(c)). After the freshwater forcing peaks, the AMOC re-invigorates, shows a short-term overshoot before it recovers to its previous steady state value. We define recovery of the AMOC to be just before the overshoot, a level that is close to the equilibrium AMOC value before and after the freshwater forcing. With reduced northward heat transport, the response in the model is a rapid cooling in the North Atlantic region and more gradual warming in the Southern Hemisphere (not shown).
The absolute temperature change as reconstructed for Greenland is strongly underestimated by the model (figure 1(b)). This is a common bias also present in other simulations with EBMs (Schmittner and Galbraith 2008) and GCMs (Bozbiyik et al 2011), presumably owing to missing atmospheric dynamics. The timescale of the temperature changes is roughly consistent with available ice core temperature reconstructions (Kindler et al 2014).
Several marine proxy observations suggest an AMOC decline and recovery during the YD Modeled marine N 2 O emissions decrease by 0.82 Tg N yr −1 over ∼1000 years, and rise back to their steady state values following the changes in AMOC with delay ( figure 1(d)). This causes a maximum drop in atmospheric N 2 O by ∼20 ppb (figure 1(e)). Reconstructions of marine N 2 O emissions over the YD (Schilt et al 2014 show a comparable shape and amplitude to the modeled N 2 O emission changes ( figure 1(d)). Note that a maximum of reconstructed marine N 2 O emissions occurs significantly before the onset of the YD and must be related to other causes.
Global changes in net N 2 O emissions are tightly coupled to changes in N 2 O production and occur almost synchronously (figure 1(e)). Most of the perturbation in net N 2 O production occurs in the upper 500 m and is communicated to the atmosphere on decadal time scales, particularly in upwelling regions.
Next, we assess the contribution of the three N 2 O production terms to the overall change in production ( figure 1(f)). There are two N 2 O production terms associated with nitrification. On the global scale, changes from the β-term (0.46 Tg N yr −1 ) dominate N 2 O production changes, followed by the α-term (0.28 Tg N yr −1 ), while changes from denitrification are small. In other words, changes in O 2 in the low-latitude thermocline affecting nitrification (β-term) are responsible for most of the N 2 O production change as further detailed in section 3.3.
Reduced overturning circulation and deep ocean ventilation lead to higher N 2 O concentrations in the deep Atlantic and Indo-Pacific ( figure 2(d)). As a result, the global marine inventory increases from 686 to 702 Tg N (not shown) despite reduced N 2 O production and negative concentration anomalies in the thermocline and a solubility driven N 2 O loss. The imbalance between global emission and production remains small as this 16 Tg N inventory change is realized over several centuries. Figure 1(i) displays the lag in years for changes in N 2 O emissions and atmospheric N 2 O as a function of corresponding percentage changes in AMOC relative to their peak change. Compared to changes in AMOC, it takes ∼300 years longer for marine N 2 O emissions to reach a ∼25% decrease. The lag decreases subsequently such that the lag of the emission maximum is only 40 years. These lags are primarily caused by the time required to build up, on average positive, O 2 anomalies in the upper ocean ( figure 1(g)), which directly drives nitrification and denitrification. We recall that the N 2 O yield per unit of oxidized organic matter decreases with increasing O 2 (Goreau et al 1980, Battaglia andJoos 2018a). The time scales for the O 2 anomalies are linked to those of thermocline ventilation and ideal age (figure 1(h)) in regions outside the North Atlantic region, where O 2 consumption rate changes are small ( figure 2(a)). The exact timing of average changes is complicated by regionally distinct, and partly offsetting changes in O 2 and ventilation (figures 2(e), (f)). On the other hand, the decrease in export production and remineralization of organic matter and associated N 2 O production (figure 1(f)) goes more or less hand in hand with the decrease in AMOC. Hence, the time for building up an O 2 anomaly after a change in ventilation and regionally differentiated changes tend to lead to a substantial lag Figure 2. Spatial anomalies in response to freshwater hosing into the North Atlantic of different N 2 O production pathways (left) and tracer concentrations (right) at the time of peak decline in marine N 2 O production. (a) ΔN 2 O production from the α-term of nitrification; multiplication with 9.40 gC (mg N) −1 yields anomalies in aerobic remineralization below 75 m in gC m −2 yr −1 and these correspond to anomalies in export production of organic matter, except in areas with denitrification, (b) ΔN 2 O production from the β-term of nitrification, (c) Δnet N 2 O production from denitrification. Anomalies in (d) N 2 O, (e) O 2 , and (f) ideal age along a section through the Atlantic (25 • W), westwards across the Southern Ocean (58 • S) and into the Pacific (175 • W). between N 2 O production and AMOC change, while the rapid response in export productivity tends to favor a short lag during AMOC decline.

Lag between atmospheric N 2 O, marine emissions and inventory, and circulation change
During recovery, global N 2 O emissions show a larger lag with respect to the recovering AMOC than during decline. In particular, the lag to full recovery increases towards 700 years. As for the build up, it takes also time to remove the O 2 anomaly ( figure 1(g)) after a change in ventilation. In addition, it takes several centuries to re-establish surface nutrient concentrations and organic matter export and associated N 2 O production after AMOC recovery ( figure 1(f)). Hence, changes in nitrification in response to O 2 as well as in response to organic matter remineralization are both delayed, enforcing a larger lag between AMOC and N 2 O emission anomalies during recovery. The modeled AMOC overshoot has little impact on export production and the ventilation in low and midlatitude regions, and, in turn, on N 2 O.
The lag of atmospheric N 2 O with respect to changes in AMOC is longer and shifted by roughly another ∼100 years compared to lag times of marine N 2 O emissions ( figure 1(i)). This is dictated by the atmospheric life time of N 2 O of ∼120 years.
Turning to the deep ocean, there are millennialscale legacy effects associated with the AMOC collapse and recovery. The timescale of peak increase in N 2 O inventory follows that of ideal age (figure 1(h)) and O 2 (figure 1(g)) in the deep ocean. At the time when modeled temperatures over Greenland have largely recovered, these model variables just reach their peak perturbation. Remineralized PO 4 3increases synchronously with decreasing O 2 in the deep, whereas a negative anomaly in preformed PO 4 3-(difference between the solid and dashed line in figure 1(g) divided by the Redfield ratio) starts to develop towards the end of the AMOC decline (figure 1(g)). Anomalies in ideal age, N 2 O and O 2 as well as in preformed and remineralized PO 4 3are widespread in the deep ocean and it takes more than 3000 years for these anomalies to disappear in the deep Pacific. PO 4 3-(figure 1(g), dashed line) reequilibrates, on average, substantially faster than these tracers as the anomalies in remineralized and preformed PO 4 3tend to cancel each other as the simulation progresses.
3.3. Spatial response to circulation changes Next, we turn to spatial anomalies at the time of peak decline in marine N 2 O production. The α-term of N 2 O production from nitrification scales linearly with changes in aerobic remineralization below 75 m. The changes in the α-term of N 2 O production (figure 2(a) and the blue line in figure 1(f)) therefore also reflect changes in aerobic remineralization and, outside regions with subsurface anoxia, export production of particulate and dissolved organic matter. Strong declines in export and aerobic remineralization are simulated for the North Atlantic region in response to less overturning and increased nutrient limitation in the surface ocean, consistent with earlier model results (Marchal et al 1998). In the Southern Ocean, slight increases in export and remineralization are simulated in response to increased temperatures and a slight retreat of sea ice. The Pacific and Indian Oceans show small changes in remineralization The β-term of N 2 O production from nitrification is most susceptible to changes in the oxygen minimum zone (OMZ) regions ( figure 2(b)). The simulated decreases result from increased O 2 in the upper ocean (figures 1(g), 2(e)). Positive anomalies in the β-term are modeled for some near-coastal regions in all basins due to increased aerobic remineralization. Decreases are largest in the equatorial OMZ of the Atlantic basin (0.22 Tg N yr −1 ), followed by changes in equatorial OMZs in the Indian (0.13 Tg N yr −1 ) and Pacific Ocean (0.11 Tg N yr −1 ). The β-term exhibits a modest, millennial-scale overshoot towards the end of the simulation caused by slightly decreased O 2 concentrations in the upper ocean of the Pacific and Indian Oceans at that time ( figure 1(f)).
The comparably small changes in net N 2 O production by denitrification show strong spatial heterogeneity ( figure 2(c)). These changes result from a complex interplay of changes in anaerobic remineralization, O 2 and marine N 2 O concentrations and are difficult to disentangle. Integrated globally, denitrification moderately contributes to the total N 2 O production changes.
The spatial pattern of O 2 anomalies is generally inversely related with those of ideal age and N 2 O (figures 2(d)-(f)). In the upper ocean, positive anomalies in O 2 and negative anomalies in ideal age and N 2 O are widespread as already discussed above. In contrast, negative anomalies in O 2 and positive anomalies in age and N 2 O are simulated in the deep. Largest changes are simulated in the deep Atlantic, but the anomalies spread into the entire deep ocean and are large in bottom waters. The collapse of the AMOC leads to a more sluggish ventilation of the deep Atlantic and the deep Indo-Pacific yielding older water mass age and leaving more time for O 2 consumption and N 2 O production by aerobic organic matter remineralization.

Discussion
We focused on the marine contribution to atmospheric N 2 O variations during the YD cold period. We benefit from a recent advance in ice core research: the separate reconstruction of marine (and terrestrial) N 2 O emissions over the deglaciation from N 2 O isotope data (Schilt et al 2014. This permits, unlike for similar glacial events (Schilt et al 2010a), a direct comparison between simulated and reconstructed marine emissions. The ice core data show that changes in marine and terrestrial emissions both contributed significantly to the atmospheric N 2 O variations. However, terrestrial N 2 O emissions responded much faster than marine emissions to the rapid warming at the end of the YD . The NH warm-cold-warm swing around the YD can be viewed as a template for similar rapid climatic events and N 2 O fluctuations during glacial periods, where a NH warm interstadial (the Bølling-Allerød) is interrupted by a short NH cold stadial (YD), followed by a DO-type rapid NH warming.
Overall, the Bern3D model is able to represent the timing and magnitude of the variation in marine N 2 O emissions around the YD. Differences between modeled and reconstructed emissions remain but cannot be avoided given the idealized setup of any freshwater hosing experiment. Bern3D was calibrated towards modern N 2 O observations (Battaglia and Joos 2018a) and applied here without further modifications. Model N 2 O production and consumption by nitrification and denitrification depends on O 2 and organic matter remineralization and accounts for differences in reaction stoichiometry between nitrification and denitrification. The model was forced with a freshwater pulse into the North Atlantic which resulted in a near-collapse and recovery of the AMOC, accompanied by adjustments in low-latitude thermocline ventilation. This in turn caused large N 2 O production and emission anomalies in upwelling systems of the Atlantic and Pacific and in the tropical Indian ocean. The largest contribution to global emission changes stems from the low-latitude ocean and is linked to anomalies in O 2 , while about a third of the production anomalies is attributed to changes in organic matter export, mainly in the North Atlantic. Similar transient negative production anomalies as simulated for YD peak conditions are projected in global warming simulations with Bern3D by 2100 CE (Common Era) (Battaglia and Joos 2018a, their figure 8(e)). The good datamodel agreement for the YD period lends support to the Bern3D projections of marine N 2 O emissions and its calibration to modern observations. Earlier marine modeling studies addressed past climate events and N 2 O and yielded conflicting findings which can now be resolved by the isotope-based reconstruction of marine versus terrestrial emissions (Schilt et al 2014. Goldstein et al (2003), using a two-dimensional ocean-atmosphere model, suggested, in agreement with the ice core reconstruction and our results, that about half of the atmospheric N 2 O perturbation during the YD is caused by changes in terrestrial emissions, whereas Schmittner and Galbraith (2008) argued that terrestrial emissions are negligible for past DO-type N 2 O variations. Their model yielded, similar to the Bern3D, a large reduction in N 2 O production in eastern upwelling systems and in the tropical Indian Ocean, but simulated, in contrast to the Bern3D, a slight increase in the northern North Atlantic in response to an AMOC collapse. Uncertainties in current models are linked to the simplified N 2 O parameterization, their coarse spatial resolution and to the simplicity of the atmospheric energy balance models which do not simulate dynamic changes in the ITCZ as associated with the YD and DO events.
The relative timing of N 2 O and ocean circulation changes is different during AMOC slow down than during its recovery in our model. During slow down, the change in marine N 2 O emission and production is somewhat delayed compared to the slowing AMOC as it takes time to build up positive O 2 anomalies in the low-latitude thermocline. The initial lag exceeds 200 years, but is then reduced to less than a few decades towards peak decline. In contrast, global N 2 O emissions have only approached about half of their equilibrium value at the time when the AMOC has resumed its old strength and it takes more than 400 additional model years towards full emission recovery. Legacy effects may play a role for a sequence of DO events as deep ocean adjustment time scales are multicentury to millennial. The modeled atmospheric N 2 O perturbations are further delayed compared to emissions owing to the 120 year atmospheric life time. as well as the limited resolution of marine records but also of the ice core N 2 O isotope record do not allow to draw conclusions on the phasing with a precision better than a few centuries . Accordingly, our results also provide an independent model-based constraint for the relative age scales between marine sediment cores and ice core gas records. In our model the decrease in atmospheric N 2 O due to marine emissions follows the slowing down in AMOC roughly within 200-300 years. However, the influence of other processes additional to AMOC change, model uncertainties and the particular choice of the freshwater forcing may affect this timing relationship.

Conclusion
In summary, the reconstructed decrease in marine N 2 O emissions during the YD is mainly explained by N 2 O production changes in the low-latitude thermocline in response to altered ventilation and O 2 concentrations after a slowing down of the AMOC. The model results suggest a century-scale lag between a decrease in the AMOC and the associated atmospheric N 2 O perturbation and a slow recovery of marine N 2 O emissions after its resumption. Taking the model results at face value, this would imply that the strong decrease in marine N 2 O emissions reconstructed from N 2 O concentrations and their nitrogen isotopic composition in ice cores for the onset of the YD (Schilt et al 2014 should be delayed relative to the onset of the AMOC reduction by a few hundred years. Vice versa, increases in marine N 2 O emissions during the rapid warmings into the Bølling-Allerød and the preboreal would suggest an increase in AMOC preceding these changes. While the isotope inversion of the ice core N 2 O record provides a quantitative reconstruction of marine and terrestrial emission changes it is still limited in the accuracy of the amplitude of the emission changes by the relatively widespread of marine and terrestrial isotope source signatures and in terms of the exact timing of marine emission changes by the resolution of the available isotope records (Fischer et al 2019). Moreover, the independent age scales of marine and ice core records do not allow to constrain the phase relationship between AMOC and marine N 2 O emission changes to better than a few centuries based on observations. In addition, atmospheric feedbacks such as changes in equatorial upwelling through changes in the position of the ITCZ in response to the rapid warmings are not considered in our model. Accordingly, the sensitivity of our lag analysis on these caveats requires further studies in the future. The good agreement between the amplitude of modeled and reconstructed emission changes for the YD climate event provides support to our model-based interpretation and to projections (Battaglia and Joos 2018a, 2018b) of future marine N 2 O emissions and related O 2 -N 2 O-climate feedbacks.