Skip to main content


Front. Mar. Sci., 09 May 2017
Sec. Marine Biogeochemistry
Volume 4 - 2017 |

A Three-Dimensional Model of the Marine Nitrogen Cycle during the Last Glacial Maximum Constrained by Sedimentary Isotopes

  • 1GEOMAR Helmholtz Centre for Ocean Research Kiel, Kiel, Germany
  • 2College of Earth, Ocean, and Atmospheric Sciences, Oregon State University, Corvallis, OR, USA

Nitrogen is a key limiting nutrient that influences marine productivity and carbon sequestration in the ocean via the biological pump. In this study, we present the first estimates of nitrogen cycling in a coupled 3D ocean-biogeochemistry-isotope model forced with realistic boundary conditions from the Last Glacial Maximum (LGM) ~21,000 years before present constrained by nitrogen isotopes. The model predicts a large decrease in nitrogen loss rates due to higher oxygen concentrations in the thermocline and sea level drop, and, as a response, reduced nitrogen fixation. Model experiments are performed to evaluate effects of hypothesized increases of atmospheric iron fluxes and oceanic phosphorus inventory relative to present-day conditions. Enhanced atmospheric iron deposition, which is required to reproduce observations, fuels export production in the Southern Ocean causing increased deep ocean nutrient storage. This reduces transport of preformed nutrients to the tropics via mode waters, thereby decreasing productivity, oxygen deficient zones, and water column N-loss there. A larger global phosphorus inventory up to 15% cannot be excluded from the currently available nitrogen isotope data. It stimulates additional nitrogen fixation that increases the global oceanic nitrogen inventory, productivity, and water column N-loss. Among our sensitivity simulations, the best agreements with nitrogen isotope data from LGM sediments indicate that water column and sedimentary N-loss were reduced by 17–62% and 35–69%, respectively, relative to preindustrial values. Our model demonstrates that multiple processes alter the nitrogen isotopic signal in most locations, which creates large uncertainties when quantitatively constraining individual nitrogen cycling processes. One key uncertainty is nitrogen fixation, which decreases by 25–65% in the model during the LGM mainly in response to reduced N-loss, due to the lack of observations in the open ocean most notably in the tropical and subtropical southern hemisphere. Nevertheless, the model estimated large increase to the global nitrate inventory of 6.5–22% suggests it may play an important role enhancing the biological carbon pump that contributes to lower atmospheric CO2 during the LGM.


The global nitrate (NO3-) inventory regulates marine productivity and sequestration of carbon in the deep ocean via sinking organic matter (i.e., the biological pump). In the modern ocean, NO3- is the main limiting nutrient throughout the tropical and subtropical oceans (Moore et al., 2013) because N-loss processes cause its depletion relative to other macronutrients such as phosphate (PO43-). Therefore, significant changes to the global NO3- inventory, which is known to be sensitive to climate (Gruber and Galloway, 2008), may feedback on climate via the biological pump.

The predominant source of bioavailable fixed nitrogen (N) to the ocean is N2 fixation by specialized cyanobacteria (diazotrophs) capable of converting dinitrogen gas to ammonium to meet their nitrogen requirements for growth. N2 fixation has extra energy requirements relative to assimilating forms of fixed nitrogen associated with breaking down the triple N bond in dinitrogen, extra respiration requirements to keep the N2-fixing cellular compartment anoxic, as well as additional structural iron requirements of the N2-fxing nitrogenase enzyme (e.g., see Karl et al., 2002; Grosskopf and Laroche, 2012).

Nitrogen loss (N-loss) processes occur in oxygen deficient zones (ODZs, O2 <~10 mmol m−3) where NO3- replaces O2 as the electron acceptor during organic matter respiration. Under these suboxic conditions, forms of fixed nitrogen (NH4, NO2, NO3-) are converted to N2O and N2 gas primarily by heterotrophic denitrification and autotrophic anammox. These N-loss processes occur both in the pore-water sediments and marine water column. Previous estimates suggest that the global ratio of sedimentary to water column N-loss is between 1 and 4 (Brandes and Devol, 2002; Altabet, 2007; Codispoti, 2007), with more recent modeling efforts suggesting a narrower range of 1.3–2.3 (Eugster and Gruber, 2012; Devries et al., 2013; Somes et al., 2013).

The balance between N2 fixation and N-loss mainly determines the preindustrial global NO3- inventory. In the preindustrial Late Holocene (0–5,000 years ago), these processes were likely close to balanced (Gruber, 2008), which has also been interpreted from the stability of sedimentary nitrogen isotope records during this period (Altabet, 2007). Some estimates have suggested that present-day N-loss processes may be much larger than N2 fixation, leading to an anthropogenic ocean that could be rapidly losing N (Codispoti, 2007). However, methodological issues with historical N2 fixation measurements suggest a significant underestimation of N2 fixation by at least a factor of 2 (Mohr et al., 2010; Großkopf et al., 2012), which could explain much of the large imbalance in some previous budgets. Recent modeling efforts (Eugster and Gruber, 2012; Devries et al., 2013) have estimated global N-loss rates (120–240 Tg N yr−1) that are on the low-end of previous studies that suggest as high as 400 Tg N yr−1 (Codispoti, 2007), again suggesting previous budgets should be closer to balanced.

Isotope ratios of nitrogen (15RN =15N/14N), typically reported as delta values (δ15N = 15RN/15RN, std−1), can be used to infer changes in physical and biogeochemical processes related to nitrogen cycling (Wada, 1980). δ15NO3- is strongly fractionated by water column N-loss, preferentially removing 14N thus increasing δ15NO3- of the remaining nitrate (εWCNl = 15–30‰; Cline and Kaplan, 1975). Sedimentary N-loss has a lower isotope effect (εSedNl = 1.5–13‰; Brandes and Devol, 1997; Lehmann et al., 2007; Granger et al., 2011; Dale et al., 2014) due to more complete NO3- utilization in pore-water sediments compared to typical water column conditions. N2 fixation introduces 15N-depleted nitrogen (δ15NNfix = −1‰; Minagawa and Wada, 1986), thereby decreasing δ15NO3-, which is ~5‰ on average globally (Sigman et al., 2000; Somes et al., 2010). Phytoplankton assimilate 14N more efficiently than 15N (εNO3-upt = 5–10‰; Wada and Hattori, 1978), which increases δ15NO3- of the residual NO3- and has a large effect in N-depleted surface waters in the absence of N2 fixation. Zooplankton excrete 15N-depleted nitrogen causing their biomass to be enriched by 3–4‰ relative to their prey (Minagawa and Wada, 1984), causing the step-wise enrichment of δ15N through higher tropic levels.

δ15N measured in the water column and on sedimentary organic matter has been used to reconstruct changes in nitrogen cycling from previous climate periods (Robinson et al., 2012). δ15N records are typically interpreted qualitatively as changes to nutrient utilization in modern High Nitrate Low Chlorophyll (HNLC) regions, water column N-loss, or N2 fixation (e.g., Altabet and Francois, 1994; Altabet et al., 1995; Ren et al., 2009). However, nitrogen isotope ratios are influenced by a variety of biogeochemical processes, sometimes at the same location, including transport by the circulation and mixing (Schmittner and Somes, 2016), which complicates their ability to quantitatively constrain nitrogen cycling processes based on individual sediment cores alone.

During the Last Glacial Maximum (LGM) ~21,000 years before present, lower temperatures have presumably increased solubility of oxygen and improved oxygen supply to the upper ocean ODZs (Jaccard and Galbraith, 2012), thereby reducing water column N-loss. Sedimentary N-loss was also likely reduced during the LGM due to exposed continental shelves from lower sea level (McElroy, 1983; Christensen et al., 1987). In the modern ocean approximately half of total sedimentary N-loss is estimated to occur on these shallow continental shelves (Bohlen et al., 2012). Previous box modeling studies using sedimentary δ15N as a constraint to estimate that total N-loss increased by ~30–120% across the last deglaciation (Deutsch et al., 2001; Eugster et al., 2013; Galbraith et al., 2013). However, more realistic, three-dimensional models have not been used so far to quantify these effects and their impacts on the global nitrogen inventory.

Effects of glacial conditions on N2 fixation are more uncertain. Falkowski (1997) suggests that increased atmospheric Fe deposition to the glacial ocean stimulated additional N2 fixation that could have increased the global NO3- inventory. Broecker (1982), and more recently (Wallmann et al., 2016) propose that a higher PO43- inventory, due to reduced loss on continental shelves, could have enhanced N2 fixation during the LGM since PO43- is diazotrophs' other main limiting nutrient. On the other hand, Ren et al. (2009, 2012) interpret δ15N records from the western tropical North Atlantic and Pacific as reduced N2 fixation, in response to lower N-loss in the glacial ocean. Some specific N2-fixing Trichodesmium species may also be limited by CO2 concentrations of surface waters (Barcelos e Ramos et al., 2007; Hutchins et al., 2007), which could have reduced N2 fixation during the LGM, although other studies investigating Trichodesmium community assemblages (Gradoville et al., 2014) and unicellular diazotrophs (Law et al., 2012) found no consistent effect of CO2 on N2 fixation.

Previous box modeling studies have produced large differences in their estimates of the global NO3- inventory during the LGM, ranging from < +10% (Deutsch et al., 2004) up to +100% (Eugster et al., 2013). The lower estimates would imply negligible impacts on the carbon cycle and atmospheric CO2, whereas the higher estimates could lead to enhanced atmospheric CO2 sequestration in the ocean. However, more realistic estimates with three-dimensional models remain lacking. While changes in the oceanic nutrient inventory have been featured in hypotheses put forward to explain the glacial drawdown of atmospheric CO2 (McElroy, 1983; Christensen et al., 1987; Falkowski, 1997; Broecker and Henderson, 1998), possible effects of changes in the nitrogen inventory on glacial-interglacial CO2 cycles are sometimes ignored in recent reviews (e.g., Hain et al., 2014) and model simulations (e.g., Brovkin et al., 2007).

Here we present the first three-dimensional model with realistic LGM boundary conditions to estimate changes in N-loss and N2 fixation and their effects on the nitrogen inventory. We build upon the isotope modeling work of Schmittner and Somes (2016), herein denoted as SS16, which used a global database of sedimentary nitrogen isotopes from the LGM (Tesdal et al., 2012) to constrain idealized 3D LGM simulations. In our global 3D model-data analysis, observations are directly compared to model results at the same locations they represent, rather than defining large, basin-scale averages that previous box modeling studies had to impose. Here we produce more realistic LGM simulations by considering spatially varying atmospheric Fe deposition and sea level effects on sedimentary N-loss, which were not accounted for in SS16. We also provide the first 3D quantitative estimates of effects of hypothesized increases in the ocean's PO43- inventory on the nitrogen cycle (Wallmann et al., 2016).

Model Description and Preindustrial Results

We use a slightly modified (see Appendix) model version of SS16, which is based on the UVic Earth System Climate Model (Weaver et al., 2001) with a version of Kiel biogeochemistry (Somes and Oschlies, 2015). In the following we provide a general overview of the model components and their preindustrial results.

Physical Model

The physical ocean-atmosphere-sea ice model includes a three-dimensional (1.8 × 3.6°, 19 vertical levels) general circulation model of the ocean (Modular Ocean Model 2) with parameterizations such as diffusive mixing along and across isopycnals, eddy-induced tracer advection (Gent and McWilliams, 1990), computation of tidally-induced diapycnal mixing over rough topography including sub-grid scale (Schmittner and Egbert, 2014), as well as anisotropic viscosity (Large et al., 2001) and enhanced zonal isopycnal mixing schemes in the tropics to mimic the effect of zonal equatorial undercurrents (Getzlaff and Dietze, 2013). Background vertical mixing is reduced in the tropical subsurface (20°S–20°N, 185–565 m), consistent with microstructure observations and tracer release experiments from the tropical Atlantic (Fischer et al., 2013). A two-dimensional, single level energy-moisture balance atmosphere and a dynamic-thermodynamic sea ice model are used, forced with prescribed monthly climatological winds (Kalnay et al., 1996) and ice sheets (Peltier, 2004). The maximum value of the Atlantic Meridional Overturning Circulation (AMOC) is 17.1 Sv (Figure 1).


Figure 1. Meridional streamfunction in preindustrial control (PIctl; A–C) and Last Glacial Maximum control (LGMctl; D–F) in the Global, Atlantic, and Indian-Pacific Oceans, respectively.

Marine Biogeochemical Model

The marine ecosystem-biogeochemical model coupled within the ocean circulation includes 2 nutrients in the inorganic (NO3- and PO43-) and organic (DON and DOP) phases, 2 phytoplankton (ordinary and N2-fixing diazotrophs), zooplankton, sinking detritus, as well as dissolved O2, dissolved inorganic carbon, alkalinity, and Δ14C (Somes and Oschlies, 2015). Iron limitation is calculated using monthly surface dissolved iron fields prescribed from the BLING model (Galbraith et al., 2010). Our model was initialized with World Ocean Atlas (WOA) (Garcia et al., 2010a,b) and the Global Ocean Data Analysis Project (GLODAP; Key et al., 2004) observations and was simulated to quasi steady-state for over 8,000 years under preindustrial boundary conditions. The basin scale patterns of Δ14C, O2, and NO3- are reproduced (Figure S2).

N2 Fixation

Diazotrophs grow slower than ordinary phytoplankton due to extra energetic demands of N2-fixation (e.g., Karl et al., 2002; Grosskopf and Laroche, 2012). However, since they have no N-limitation, they can out-compete ordinary phytoplankton in surface waters that are depleted in NO3-, but still contain sufficient P and Fe (e.g., water with low NO3- from denitrification and high iron from atmospheric Fe deposition). They will consume NO3- if it exists at high concentrations (>5 mmol m−3), consistent with culture experiments (Mulholland et al., 2001). Diazotrophs consume DOP in the model when PO43- limited, which expands their ecological niche (Somes and Oschlies, 2015). Zooplankton graze diazotrophs with a lower preference compared to ordinary phytoplankton in the model, qualitatively supported by observations (O'Neil, 1999).


Water column N-loss occurs when dissolved oxygen becomes depleted in poorly ventilated ODZs. NO3- replaces dissolved oxygen as the electron acceptor during organic matter respiration and is reduced to dinitrogen gas. We use a threshold of 5 mmol m−3 O2 that sets where respiration of organic matter occurs equally between denitrification and aerobic respiration. Further above (below) this threshold, a greater fraction of aerobic respiration (denitrification) occurs with complete aerobic respiration above 10 mmol m−3 O2. Although, the ODZs are misplaced too close to the equator in the Pacific and Atlantic, and in the Indian ocean they are more intense in the Bay of Bengal rather than in the Arabian Sea in contrast to observations, the model reproduces the general regions and sizes of ODZs and global water column N-loss rates are within observational uncertainties (Bianchi et al., 2012).

Sedimentary N-loss is simulated according to an empirical transfer function based on organic carbon sinking flux to the sediments and bottom-water dissolved oxygen and nitrate. We apply the empirical function from Bohlen et al. (2012), which follows the general methodology of Middelburg et al. (1996), but uses additional observations to constrain their model. These empirical functions are computationally efficient alternatives to coupling a full sediment model. This function also accounts for nitrogen flux from sediment porewater to oceanic bottom waters, which reduces the net sedimentary N-loss by ~25% over the upper 1,000 m. Therefore, it should be considered as a net sedimentary N-loss rather than a direct sedimentary denitrification rate.

Since our coarse resolution model does not fully resolve narrow continental shelves and coastal dynamics, it underestimates sedimentary N-loss in these regions. To account for this model deficiency, we accelerate sedimentary N-loss only within the subgrid-scale bathymetry scheme (Figure S1) in the upper 250 m. The model's coarse resolution poorly represents physical processes (e.g., upwelling), which ultimately drive sedimentary N-loss, to the largest degree on unresolved narrow shelves. Therefore, we focus the sedimentary N-loss acceleration there using the equation accelSedNl = [(1FSGS) × ΓSedNl + 1] × SedNl, where accelSedNl is the accelerated sedimentary N-loss rate used in the model, FSGS is the fraction of the model grid covered by the subgrid-scale scheme (Figure S1), SedNl is the true sedimentary N-loss function from Bohlen et al. (2012), and ΓSedNl is the arbitrary acceleration factor. This formulation accelerates sedimentary N-loss the most when the subgrid-scale fraction is lowest (i.e., narrow shelf) where coastal processes are likely most unresolved in the model. It has no acceleration effect when the subgrid-scale fraction covers the entire model grid bathymetry where the model should better resolve net processes occurring over the sea floor. The sedimentary N-loss acceleration factor is chosen such that our global model rate is within the low-end uncertainty range from previous observationally constrained model estimates (Eugster and Gruber, 2012; Devries et al., 2013). This acceleration formulation increases global sedimentary N-loss by 29 Tg N yr−1 to yield a global rate of 91 Tg N yr−1 in the preindustrial control run (PIctl, see below).

Nitrogen Isotopes

Nitrogen isotope ratios are affected by inventory-altering (N2 fixation and N-loss) and internal-cycling (NO3- uptake, excretion, DON remineralization) processes in the model (Somes et al., 2010). N2 fixation introduces isotopically light atmospheric nitrogen into the ocean (δ15NNfix = −1‰), whereas N-loss fractionates strongly in the water column (εWCNl = 20‰) and moderately in the sediments (εSedNl = 6‰). NO3- uptake by phytoplankton fractionates NO3- at 6‰ in the model. Zooplankton excretion fractionates at 4‰ enriching its biomass in 15N relative to phytoplankton. DON remineralizes with a fractionation factor of 1.5‰ to reproduce upper ocean δ15N-DON observations mainly occurring within the range of 3–6‰ (Knapp et al., 2011).

The patterns and rates of N2 fixation and N-loss are generally consistent with its observational biogeochemical indicators N* and sedimentary δ15N (Figure 2). Water column N-loss produces low N* and high δ15N in the ODZs of the eastern tropical Pacific and northern Indian Ocean. Sedimentary N-loss delivers similar, but more moderate trends most notably on large continental shelves (e.g., Bering Sea, Patagonian Shelf). In the more oligotrophic tropics and subtropics, N2 fixation conversely generates high N* and low δ15N. Low NO3- utilization due to iron limitation in HNLC regions of the Southern Ocean, North Pacific, and eastern equatorial Pacific drive low δ15N, whereas high utilization in the oligotrophic causes higher δ15N in the absence of N2 fixation.


Figure 2. Preindustrial control (PIctl) annual vertically-integrated (A) N2 fixation and (B) water column plus sedimentary N-loss. Comparison of upper ocean (0–300 m) N* and seafloor δ15N of PIctl (C,D) to (E) World Ocean Atlas 2009 (WOA09) (Garcia et al., 2010b) and (F) NICOPP database (Tesdal et al., 2013) observations. (B–D) Contours of minimum dissolved oxygen (10 mmol m−3) in the water column show where water column N-loss occurs.

Global mean δ15NO3- is determined by global rates and isotope effects of N2 fixation, and water column and sedimentary N-loss (Brandes and Devol, 2002; Altabet, 2007). Sedimentary N-loss rates and its isotope effect represent the largest uncertainties for determining global mean δ15NO3- (Devries et al., 2013; Somes et al., 2013). Therefore, we took an empirical approach and chose the sedimentary N-loss fractionation factor (εSedNl = 6‰) within its estimated uncertainity range (1.5–13‰; Brandes and Devol, 1997; Lehmann et al., 2007; Granger et al., 2011; Dale et al., 2014) that sets deep ocean δ15NO3- near its observed ~5‰ value given the model predicted rates of N2 fixation, water column N-loss, and sedimentary N-loss described above.

Last Glacial Maximum Model Configurations

Similar to Schmittner and Somes (2016) and guided by the Paleoclimate Model Intercomparison Project (PMIP3; Braconnot et al., 2012), we prescribe the LGM boundary conditions of the model: lower atmospheric concentrations of the greenhouse gases carbon dioxide, nitrous oxide, and methane, changes of Earth's orbit, and the increased area and height of ice sheets (Peltier, 2004). The ocean grid bathymetry and total ocean volume remains unchanged relative to the preindustrial simulation. However, effects of reduced sea level on sedimentary N-loss are accounted for by calculating a new sub-grid scale bathymetry scheme assuming a constant 120 m sea level reduction (Figure S1). This is different from SS16, who neglected sea level effect. The marine biogeochemistry equations are identical in the preindustrial and LGM simulations with parameters given in Table 1.


Table 1. Model configurations.

Wind stress anomalies (LGM minus preindustrial) are prescribed using monthly averages of 8 models that participated in the PMIP3 (Muglia and Schmittner, 2015) since the UVic model does not include an atmospheric general circulation component (Figure S3). This is different from SS16, who used pre-industrial wind stress. As shown by Muglia and Schmittner (2015), LGM wind stress anomalies increase the Atlantic Meridional Overturning Circulation (AMOC). We further deviate from SS16 by reducing moisture diffusivity in the 2D energy-moisture balance atmospheric model across the Southern Ocean by a factor of 2. Reduced meridional atmospheric moisture flux has been suggested as an important mechanism to influence glacial stratification and carbon storage in the Southern Ocean (Sigman et al., 2007), which weakens and shoals the AMOC in our model. The resulting AMOC is weaker (~9.9 Sv maximum) and shallower (by ~500 m as indicated by shoaling of the zero streamfunction line in Figures 1B,E from ~2,700 to ~2,200 m), similar to that of SS16 due to the compensating effects of wind stress and moisture diffusivity changes. The simulations were initialized with results from SS16 and run for over 8,000 years as they approached their quasi steady-state solution.

Atmospheric Fe Deposition

We account for increased atmospheric Fe deposition during the LGM (Figure 3). LGM and preindustrial dust flux estimates were obtained from Albani et al. (2014). We assumed a 3.5% Fe content in dust, and variable solubility was used from an atmospheric model that simulates the relation between modern iron flux and solubility of iron in dust (Luo et al., 2008). We applied their relationship to LGM and preindustrial dust fluxes, thereby obtaining surface soluble atmospheric Fe fluxes (Figure S4).


Figure 3. Annual euphotic zone dissolved iron in (A) preindustrial control (PIctl) and the ratio relative to preindustrial of the LGM simulations (B) LGMloFe, (C) LGMctl, and (D) LGMhiFe based on LGM atmospheric deposition estimates (Albani et al., 2014).

To predict the surface soluble atmospheric Fe fluxes effect on dissolved Fe concentrations in our model, we analyzed the global relationship between soluble atmospheric Fe flux (atmFeflx) and dissolved Fe (dFe) in the upper ocean from a different model that includes Fe as a prognostic tracer (Nickelsen et al., 2015). The relationship that yielded the highest correlation (R = 0.33) was a power law of the form dFe = A × (atmFeflx)B where A = 6.3 × 10−4 and B = 0.24 (Figure S4). LGM dissolved Fe was calculated as dFeLGM = (atmFeflxLGM/atmFeflxPI)B× dFePI. This pragmatic approach is used due to high uncertainties associated with atmospheric dust deposition rates and Fe solubility, responsible for soluble Fe deposition estimates that differ by over an order of magnitude in different preindustrial global marine iron models (Tagliabue et al., 2016). In our LGM Fe sensitivity experiments, we changed the exponent B by the value ±0.2 from the LGM control value 0.25, which yields ± ~20% changes to global euphotic zone dissolved Fe concentrations (Table 1, Figure 3).

PO43- Inventory

We test a hypothesis that the global PO43- inventory was larger during the LGM mainly due to reduced organic phosphorus burial on the exposed continental shelves from lower sea level (Broecker, 1982; Wallmann, 2010; Wallmann et al., 2016). Since the global oceanic PO43- inventory is conserved in each model simulation, we implement this hypothesis by imposing the estimated + ~15% increase in the global PO43- inventory from Wallmann et al. (2016) in LGMctl (Table 1). Model experiments #4 LGMloPO43- and #5 LGMhiPO43- probe its sensitivity to ±10% global PO43- change relative to LGMctl.

N-Loss Rates

We perform additional simulations to investigate how different water column and sedimentary N-loss rates affect LGM nitrogen cycling and δ15N distributions. In the model, water column N-loss can be increased pragmatically by increasing the elemental ratio of oxygen consumed per organic nitrogen remineralized (RO:N, see Table 1). This leads to more oxygen loss, lower oxygen concentrations and larger ODZs in the model, thereby increasing water column N-loss. We interpret these simulations as water column N-loss sensitivity experiments rather than predictions of higher RO:N stoichiometry, however noting that Devries and Deutsch (2014) suggest this ratio is higher in more stratified, nutrient depleted surface regimes that are more prevalent features in the LGM. Model experiment #6 (LGMmidWCNl) and #7 (LGMhiWCNl) assume RO:N ratios of 11.0 and 11.4, respectively, compared to the standard value of 10.6 used in the preindustrial and LGM control simulations.

Sedimentary N-loss has been altered separately by changing its rate only on the subgrid-scale bathymetry scheme. We have applied an acceleration parameterization focused on narrow continental shelves in the upper 250 m to account for poorly resolved coastal processes (e.g., upwelling) that fuel sedimentary N-loss there (Figure S1). In the control simulations, applying an acceleration factor of 2 in this scheme increases sedimentary N-loss by 29 Tg N yr-1 (44%) in PIctl and 11 Tg N yr−1 (29%) in LGMctl. LGM model experiments #8 LGMloSedNl and #9 LGMhiSedNl demonstrate the model sensitivity to this parameterization with lower and higher factors (Table 1), respectively, which altered global rates by ± ~40% relative to LGMctl (Table 2) in this sensitivity test. Note that since LGMloSedNl does not contain the sedimentary N-loss acceleration formulation, its factor is applied directly to the sedimentary N-loss rate on the subgrid-scale scheme. These N-loss sensitivity simulations were initialized with LGMctl and run for an additional 4,000 years.


Table 2. Global biogeochemical model results.

LGM Model Results

In this section, we first provide a brief overview of the physical climate since the focus of this study is on marine biogeochemistry, noting that the physical climate (e.g., temperature and circulation) is identical for all LGM simulations. The subsequent subsection describes the LGMctl biogeochemistry relative to PIctl, followed by the LGM sensitivity experiments.

Physical Climate

The model reproduces major physical characteristics of the LGM climate, which are identical in all LGM simulations. Global average surface air temperature is 4.3°C colder in the LGMctl. The AMOC is shallower and weaker (9.9 Sv) compared to PIctl (17.1 Sv; Figure 1). Global mean Δ14C values are −165 and −225‰ in PIctl and LGMctl, respectively. Bottom waters in the LGMctl simulation are colder everywhere and saltier except for the eastern North Atlantic, where they are fresher compared with PIctl (Figure S5) due to the shoaling and weakening of the AMOC since Antarctic Bottom Water is fresher than North Atlantic Deep Water in the model. The simulated pattern of salinity changes indicates enhanced brine rejection from intensified sea ice formation in the Weddell Sea. The basin scale simulated changes in bottom water properties are qualitatively consistent with observations (Figure S5; Adkins et al., 2002; Insua et al., 2014). The simulated decrease of Δ14C below 2 km depth (−75‰, ~800 years radiocarbon age) is slightly older than a recent estimate of ~600 years based on sediment reconstructions (Sarnthein et al., 2013).

LGM Control Marine Biogeochemistry

In LGMctl, oxygen concentrations are higher above ~1,500 m and lower below (Figures 4A,B) relative to PIctl, which is consistent with O2 reconstructions (Jaccard and Galbraith, 2012). The colder LGM climate enhances sea surface solubility of oxygen that is responsible for the upper ocean increase (Meissner et al., 2005). The combination of the more sluggish AMOC (Figure 1C) and increased sea ice cover, which reduce oxygen supply to the deep ocean, and enhanced atmospheric Fe deposition (Figure 3), which fuels organic matter export production that increases subsurface oxygen consumption, are the main factors responsible for the decreased deep ocean oxygen concentrations. On the global average, there is a 21.5% reduction of the global O2 inventory (Table 2). This is similar to the 20.8% reduction resulting from SS16's best δ15N simulation, who attribute ~7% of this decrease to physics and ~14% to changes in Fe fertilization.


Figure 4. Zonal dissolved oxygen differences in the Atlantic and Indian-Pacific Oceans, respectively, of (A,B) Last glacial Maximum control (LGMctl) minus preindustrial control (PIctl) and the Last Glacial Maximum simulations minus LGM control: (C,D) LGMloFe, (E,F) LGMhiFe, (G,H) LGMloPO43-, (I,J) LGMhiPO43-, (K,L) LGMmidWCNl, (M,N) LGMhiWCNl, (O,P) LGMloSedNl, and (Q,R) LGMhiSedNl.

Higher O2 concentrations in the upper ocean reduce the volume of ODZs by 47% and water column N-loss by 78% (Figures 5B, 6). The difference in the reduction of ODZ volume compared to water column N-loss rate occurs due to the emergence of a ODZ in the deep (1,000–2,000 m) subarctic Pacific that accounts for 10% of total water column N-loss in LGMctl, which cannot be excluded because trace metal proxies at bottom waters in these locations indicate less oxygen during the LGM (Jaccard and Galbraith, 2012). Here, OM remineralization rates that drive water column N-loss are lower than in the shallower tropical ocean, which results in less water column N-loss per volume of ODZ waters. In combination with decreased sedimentary N-loss by 55% due to exposed continental shelves from lower sea level in LGMctl (e.g., Bering Sea Shelf, North Atlantic; Figure 5B), total N-loss decreases by 59% (Table 2, Figure 6). However, note that some regions at mid-latitudes in the Southern Ocean, North Pacific, and North Atlantic show increases in sedimentary N-loss due enhanced atmospheric Fe deposition fueling additional export production (Figure 5B). Although reduced total N-loss occurs primarily in the upper ocean, the resulting increased NO3- accumulates in the deep ocean due to enhanced POM export in the Southern Ocean that remineralizes and remains for longer in the older, more sluggish Antarctic Bottom Waters in LGMctl.


Figure 5. Vertically-integrated differences of N2 fixation and total N-loss, respectively, of (A,B) Last Glacial Maximum control (LGMctl) minus preindustrial control (PIctl) and the Last Glacial Maximum simulations minus LGM control: (C,D) LGMloFe, (E,F) LGMhiFe, (G,H) LGMloPO43-, (I,J) LGMhiPO43-, (K,L) LGMmidWCNl, (M,N) LGMhiWCNl, (O,P) LGMloSedNl, and (Q,R) LGMhiSedNl. Contours of minimum dissolved oxygen (10 mmol m−3) in the water column show where water column N-loss occurs.


Figure 6. (Top) Global N2 fixation (left blue bar), water column N-loss (right red bar), and sedimentary N-loss (right orange bar) for the simulations. (Bottom) Global NO3- (left green bar) and global δ15NO3- (right magenta bar) for the simulations.

N2 fixation decreases in the LGMctl simulation by 59% mainly in response to the reduction in N-loss when approaching the steady-state solution (Figures 5, 6). Only the oligotrophic eastern tropical Pacific and southern boundary of the southern subtropical gyres show subtle increases of N2 fixation due to enhanced atmospheric Fe deposition (Figure 5A). In the LGMctl model, global N2 fixation is not sensitive to increased Fe deposition because glacial estimates of increased Fe deposition applied here (Albani et al., 2014) primarily occur at high latitudes where diazotrophs are limited by temperature. However, the combined effects of the changes in sources and sinks result in a significantly increased global NO3- inventory (+22%). Note that DOP utilization by diazotrophs does not play a significant role stimulating N2 fixation in LGMctl because nitrogen limitation is reduced due to the larger NO3- inventory and thus ordinary phytoplankton outcompete diazotrophs for the majority of DOP consumption.

Atmospheric Fe Deposition Sensitivity

Atmospheric Fe deposition was increased during the LGM compared to pre-industrial times, which has a large fertilizing effect on organic matter production in the HNLC regions of the Southern Ocean, subarctic North Pacific, and equatorial eastern Pacific in the model. Additional NO3- is thus consumed at the surface and stored in the deep ocean (Figures 7A,B), where it can remain for millennia.


Figure 7. Zonal dissolved nitrate (NO3-) differences in the Atlantic and Indian-Pacific Oceans, respectively, of (A,B) Last glacial Maximum control (LGMctl) minus preindustrial control (PIctl) and the Last Glacial Maximum simulations minus LGM control: (C,D) LGMloFe, (E,F) LGMhiFe, (G,H) LGMloPO43-, (I,J) LGMhiPO43-, (K,L) LGMmidWCNl, (M,N) LGMhiWCNl, (O,P) LGMloSedNl, and (Q,R) LGMhiSedNl.

Enhanced export production at high latitudes reduces the amount of preformed surface NO3- in Subantarctic Mode Waters that flow to the tropics (Figures 7C–F). This reduces tropical productivity, volume of tropical ODZs and global water column N-loss in LGMhiFe relative to LGMloFe (Table 2). Intermediate waters (~1,000–2,000 m) in the subarctic North Pacific are near the threshold of suboxia in simulation LGMctl. Those suboxic regions intensify when enhanced atmospheric Fe deposition in LGMhiFe stimulates additional organic matter remineralization (Figure 4F), which enhances water column N-loss (Figure 5F).

Only minor increases to N2 fixation are stimulated by enhanced atmospheric Fe deposition in the model (Figure 5E). This is most pronounced in two zones in the southern hemisphere that are largely Fe limited in the preindustrial ocean: the eastern tropical Pacific and a zonal band across the edge of the southern subtropical gyres. Since diazotrophs prefer warm oligotrophic N-limited surface conditions, they grow outside the cold, high surface NO3- waters of the eastern equatorial Pacific upwelling region. Although, not typically thought to occur outside the tropics, the strong atmospheric Fe deposition originating around Patagonia (Figure 3) is sufficient to stimulate N2 fixation in our model. However, on the global average, these areas of additional N2 fixation are much smaller than areas in which N2 fixation decreased in model LGMctl relative to PIctl due to reduced N-loss (Figure 5A).

Water column δ15NO3- is strongly influenced by nitrogen cycling changes in response to enhanced atmospheric Fe deposition (Figure 8). Enhanced (reduced) surface NO3- utilization increases (decreases) surface δ15NO3- in the Southern Ocean and Subarctic Pacific in LGMhiFe (LGMloFe). The minor increase (decrease) to N2 fixation rates in the southern tropics and subtropics in LGMhiFe (LGMloFe) substantially decreases (increases) δ15NO3-. This large isotopic signal occurs because these systems switch from being driven by high NO3- utilization, which causes extremely 15N-enriched surface NO3-, to a N2 fixation isotopic regime that introduces 15N-depleted nitrogen to the surface ocean. This model prediction suggests that all relevant processes should be considered when interpreting δ15N and the change to the isotopic signal may not be linearly representative of the rate change to one specific nitrogen cycling process alone.


Figure 8. Zonal dissolved δ15N of nitrate (δ15NO3-) differences in the Atlantic and Indian-Pacific Oceans, respectively, of (A,B) Last glacial Maximum control (LGMctl) minus preindustrial control (PIctl) and the Last Glacial Maximum simulations minus LGM control: (C,D) LGMloFe, (E,F) LGMhiFe, (G,H) LGMloPO43-, (I,J) LGMhiPO43-, (K,L) LGMmidWCNl, (M,N) LGMhiWCNl, (O,P) LGMloSedNl, and (Q,R) LGMhiSedNl.

Global PO43- Inventory Sensitivity

The model simulations testing a hypothesized increase to the global PO43- inventory predict larger NO3- inventories due to enhanced N2 fixation (Table 2, Figures 6, 7G–J). PO43- is the major limiting macronutrient for N2-fixing diazotrophs. Although N-loss is reduced during the LGM, N-limitation still persists, leaving the tropical and subtropical ocean mostly N limited. Therefore, increased PO43- stimulates additional growth of diazotrophs via N2 fixation that is most pronounced in the tropical ocean. The additional organic matter production via enhanced N2 fixation and remineralization at depth stimulates additional water column N-loss, which balances the enhanced N2 fixation at a new steady-state with a larger global NO3- (Figure 6) and reduced O2 inventory (Figures 4I,J). A 10% increase in the PO43- inventory thus causes an 8% increase in the NO3- inventory (LGMctl-LGMloPO43-). However, a further 10% increase in the PO43- inventory (LGMhiPO43--LGMctl) leads to diminishing returns in terms of NO3- (4%) due to the further expansion of ODZs and N-loss. The additional N2 fixation decreases surface δ15NO3- in the tropics and subtropics, whereas the additional water column N-loss increases subsurface δ15NO3- in the Indian and Pacific (Figures 8I,J).

N-Loss Sensitivity

The water column N-loss sensitivity simulations alter rates by increasing oxygen consumption during organic matter remineralization, thereby decreasing subsurface oxygen concentrations (Figures 4K–N). Note that these simulations do not change the large-scale pattern of higher oxygen in the upper ocean and lower oxygen in the deep ocean that is supported by observational constraints (Jaccard and Galbraith, 2012), but only reduces the magnitude of this qualitative trend in the deep ocean. Enhanced water column N-loss in the expanded ODZs in LGMmidWCNl and LGMhiWCNl reduces NO3- in the eastern tropical Pacific, North Indian, and subarctic Pacific (Figures 7K–N). This significantly increases subsurface δ15NO3- since water column N-loss has a strong isotope effect (Figures 8K–N). Even though this signal is generated in the ODZs of the tropical ocean and subarctic Pacific, it eventually reaches the Southern Ocean and has a global impact (Figures 8K–N).

Sedimentary N-loss has lower local rates that are spread throughout all ocean basins so its signal is more distributed through the global ocean compared to vertically-integrated water column N-loss. Enhanced (reduced) rates in LGMhiSedNl (LGMloSedNl) decreases (increases) global NO3- (Figures 7O–R). Since sedimentary N-loss has a relatively small isotope effect, its main isotopic signal is an indirect one by stimulating additional (less) N2 fixation to balance additional (less) sedimentary N-loss in LGMhiSedNl (LGMloSedNl). This decreases (increases) surface δ15NO3- most notably in the tropical and subtropical ocean basins (Figures 8O–R).

Nitrogen Isotope Model-Data Analysis

Modeled nitrogen isotopes are compared to 61 datapoints from the recent synthesis of reconstructions from the Nitrogen Cycling in the Oceans Past and Present (NICOPP) group (Galbraith et al., 2013) for which LGM—Late Holocene differences are available. These data, based on bulk organic nitrogen measurements, were corrected for a diagenetic effect according to Robinson et al. (2012), who examined over 100 locations and found a significant correlation with water depth (~1‰/km), while noting that oxygen exposure time is likely the main driver. We use the five additional foraminifera and diatom bound measurements included in Schmittner and Somes (2016), plus 6 more from recent studies in the subarctic Pacific and central equatorial Pacific (Ren et al., 2015; Costa et al., 2016), all of which were not corrected for diagenesis. Although combining different forms of δ15N records creates uncertainty, excluding the diatom and foraminifera bound measurements does not affect which simulations performed best in terms of the statistical model-data metrics so we include them into the main analysis to improve spatial coverage.

Overall, δ15N records from LGM sediments document lower values in modern ODZs such as the eastern tropical Pacific and Arabian Sea indicating less water column N-loss, higher values in the western tropics of the North Pacific and North Atlantic indicating less N2 fixation, and higher values in high latitude HNLC regions (i.e., Southern Ocean, western subarctic Pacific) indicating enhanced surface NO3- utilization (Figures 911). These sedimentary δ15N records are compared to the δ15N values of sinking particulate organic nitrogen (PON) in bottom waters above the seafloor, which we refer to as sedimentary δ15N in the model below. Our model does not account for terrestrial sources of δ15N, which can influence records near the continental margin and potentially contribute to model-data misfits. Therefore, we focus our analysis and discussion on regions that are strongly influenced by marine nitrogen cycling processes, but note that all data points are included in the model-data statistical metrics.


Figure 9. Sedimentary δ15N of (A) LGM minus preindustrial observations and (B) Last Glacial Maximum control (LGMctl) minus preindustrial control (PIctl). Sedimentary δ15N differences of the Last Glacial Maximum simulations: (C,D) LGMloFe, (E,F) LGMhiFe, (G,H) LGMloPO43-, (I,J) LGMhiPO43-, (K,L) LGMmidWCNl, (M,N) LGMhiWCNl, (O,P) LGMloSedNl, and (Q,R) LGMhiSedNl minus PIctl and minus LGMctl are shown, respectively.


Figure 10. Comparison of LGM minus PI changes in sedimentary δ15N from observations (horizontal axis) and model results (vertical axis) of (A) LGMctl, (B) LGMloFe, (C) LGMloPO43-, (D) LGMhiFe, (E) LGMhiPO43-, (F) LGMmidWCNl, (G) LGMloSedNl, (H) LGMhiWCNl, and (I) LGMhiSedNl. Locations in the Pacific, Atlantic, and Indian Oceans are shown as circles, squares, and triangles, respectively. Locations in oxygen deficient zones (ODZ; minimum O2 < 10 mmol m−3) are shown as blue crosses and Southern Ocean (south of 40°S) are red Xs. The black line indicates a “perfect” representation of observations.


Figure 11. Regional Last Glacial Maximum minus Preindustrial δ15N differences of LGM model simulations (#1) LGMctl, (#2) LGMloFe, (#3) LGMhiFe, (#4) LGMloPO43-, (#5) LGMhiPO43-, (#6) LGMmidWCNl, (#7) LGMhiWCNl, (#8) LGMloSedNl, (#9) LGMhiSedNl in the (A) global ocean, (B) Southern Ocean (<40°S), (C) oxygen deficient zones (minimum O2 <10 mmol m−3), (D) Central Equatorial Pacific (170°W–130°W, 5°N–5°S), (E) Western Tropical North Atlantic and Pacific (0–25°N, western 50° of Atlantic and Pacific basins), and (F) Subarctic Pacific (40–65°N). “Data-masked” model (dark blue left bar) only includes locations where data exist, whereas the “full” model (light blue right bar) averages over the entire region.

LGM Control Simulation

Model simulation LGMctl, which includes enhanced atmospheric Fe deposition and assumes a +15% increase to the global PO43- inventory, reproduces the major trends in the LGM reconstructions (R = 0.45, Table 3; Figures 911). It reproduces the largest increase in δ15N in the Southern Ocean with respect to preindustrial values (Figure 11B) and the smaller increase in the subarctic Pacific HNLC region (Figure 11F). However, it predicts a slight increase (+0.7‰) to the Central Equatorial Pacific, where the data show a minor decrease (−0.2‰). The largest model-data mismatch exists in the Eastern Tropical North Pacific ODZ, where the predicted decrease of δ15N is too strong relative to observations (Figure 11C).


Table 3. δ15N model-data statistical metrics.

The change to global δ15NO3- (Table 2) and sedimentary δ15N (Figure 11A) in LGMctl relative to PIctl is small compared to the change of the ratio of sedimentary to water column N-loss, which in part determines global δ15NO3- (Somes et al., 2013). This relatively small δ15NO3- change is due to the dilution effect in N-loss zones (Deutsch et al., 2004), which increases the isotope effect of water column N-loss (per mole N removed) in LGMctl relative to PIctl as NO3- utilization decreases in the smaller ODZs. This effect decreases the isotopic signature of δ15N removed by water column N-loss from −5.1‰ in PIctl to −10.3‰ in LGMctl, thereby increasing its isotope effect by 5.2‰. This also occurs to a lesser degree in bottom waters where sedimentary N-loss occurs increasing its isotope effect by 1.1‰. Therefore, additional N2 fixation and thus a higher ratio of sedimentary to water column N-loss is required to balance the increased isotope effect from N-loss on global δ15NO3- in the LGMctl steady-state solution compared to PIctl.

The processes that lead to slightly higher deep ocean δ15NO3- in LGMctl are mainly controlled at high latitudes. In the North Atlantic, reduced N2 fixation elevates surface δ15NO3- that eventually subducts in the Atlantic meridional overturning circulation and increases deep ocean δ15NO3- (Figure 8A). In the Southern Ocean, enhanced atmospheric Fe deposition increases NO3- utilization and the isotopic signature of sinking particulate organic nitrogen (PON) that remineralizes at depth (Figures 8A,B, 9B). The emergence of a small ODZ in the deep (1,000–2,000 m) subarctic Pacific (Figure 5B), which accounts for 10% of total water column N-loss in LGMctl, contributes to elevated deep ocean δ15NO3- as well. Since most of the increased NO3- inventory in LGMctl accumulates in these old, sluggish deep waters, our modeling suggests that high latitude process also play a role determining global mean δ15NO3-, not only N2 fixation and water column N-loss in the tropics.

Atmospheric Fe Deposition Sensitivity

Lower atmospheric Fe deposition in LGMloFe compared to LGMctl reduces NO3- utilization and sedimentary δ15N in HNLC regions (Figures 9C,D), where the fit with observations deteriorates (Figures 10B,11B). The better agreement of the magnitude of the increase in δ15N in LGMctl in the Southern Ocean and Subarctic Pacific supports LGMctl's assumption for increased Fe deposition during the LGM in those regions. In fact, LGMctl's increase in NO3- utilization may be underestimated in the Southern Ocean, where simulation LGMhiFe agrees better with the observations (Figure 11B).

In the Central Equatorial Pacific, which is on the western border of the eastern Pacific HNLC region, both LGM Fe sensitivity simulations reproduce the observations better than LGMctl (Figure 11D), but for different reasons. In LGMloFe, the small increase to atmospheric Fe deposition was insufficient to increase NO3- utilization and hence sedimentary δ15N in this region. In our high-end simulation for atmospheric Fe deposition (LGMhiFe), enhanced N2 fixation introducing additional 15N-depleted nitrogen reduces the overestimated model bias of LGMctl. Until more observations in and around the central and eastern tropical Pacific HNLC become available, it will remain difficult to evaluate whether relatively low NO3- utilization or increased N2 fixation was the most important process driving the observed decreased LGM δ15N trend.

Our sensitivity model experiments in the Southern Hemisphere show a balance between lower δ15N due to enhanced N2 fixation in the subtropics and higher δ15N from enhanced surface NO3- utilization toward high latitudes. Our model predicts that increased Fe could stimulate additional N2 fixation at the poleward edge of the southern subtropics as far south as ~40°S, which counteracts with the isotope effect of enhanced NO3- utilization that increases δ15N at higher latitudes (Figures 8C–F, 9C–F). As NO3- decreases with enhanced atmospheric Fe deposition, the additional 15N-depleted nitrogen via N2 fixation contributes to a larger fraction of total sinking PON, which prevents a higher δ15N increase in LGMhiFe compared to the δ15N decrease in LGMloFe (Figure 11B). The absence of open ocean δ15N observations in the tropical and subtropical southern hemisphere makes it uncertain to validate this increase to N2 fixation there.

Global PO43- Inventory Sensitivity

In the model, a larger global PO43- inventory stimulates increases to N2 fixation, water column N-loss, and the global NO3- inventory (Figure 6). This results in increased sedimentary δ15N near ODZs and decreased δ15N outside of them in the more oligotrophic ocean where additional N2 fixation occurs (Figure 9J). The magnitude of change to δ15N in the model (Figure 10) is relatively small because N2 fixation and water column N-loss have counteracting nitrogen isotopic effects. Since these processes occur within close spatial proximity in the tropics, the input of 15N-depleted nitrogen via fixation roughly compensates the extra 15N-enriched NO3- via enhanced water column N-loss rates, preventing a larger change to the large scale δ15N signal (Figures 9G–J).

Despite large changes to the global rates of these processes, changes to simulated δ15N at the location of observations only moderately affect the model-data statistical metrics (Table 3, Figure 10). The higher water column N-loss in LGMhiPO43- better reproduces δ15N in tropical ODZs (Figure 11C). However, the additional water column N-loss in the subarctic Pacific overestimates δ15N in that region. The model experiment LGMloPO43- performs best in the Central Equatorial Pacific (Figure 11E), which suggests that the additional PO43- stimulates too much N2 fixation there in simulations LGMctl and LGMhiPO43-.

N-Loss Rate Sensitivity

Since one of the main model deficiencies in LGMctl is underestimated δ15N in ODZs (Figures 10, 11C), we performed sensitivity simulations with higher water column N-loss rates (#7 LGMmidWCNl and #8 LGMhiWCNl). Imposing moderately higher water column N-loss rates increases δ15N in existing ODZ zones in the Eastern Pacific and North Indian Ocean and reduces the model-data mismatch in those regions (Figure 11C). However, the enhanced water column N-loss rates also increases δ15N in the Central Equatorial Pacific and Subarctic Pacific, causing those regions to overestimate δ15N compared to the observations, demonstrating how water column N-loss influences δ15N across the entire tropical Pacific and beyond (Figures 8K–N, 9K–N,P).

Higher water column N-loss rates increase global δ15NO3- (Figures 8K–N), but this does not impact sedimentary δ15N everywhere equally (Figures 9K–N). Effects are largest in regions with upwelling of the increased deep ocean δ15NO3-, e.g., tropical Pacific, subarctic North Pacific, and Southern Ocean. Our simulations suggest that if a significant increase in global δ15NO3- occurred during the LGM, this could significantly impact δ15N across the Southern Ocean, which are typically interpreted as changes to surface NO3- utilization alone. Given the sparsity of the global dataset, we cannot exclude a significant increase in global δ15NO3- during the LGM because its large-scale impact affects sedimentary δ15N in regions that already have an increased δ15N trend relative to PI (e.g., Southern Ocean, subarctic Pacific). The model results suggest that the distribution of existing sites is such that calculating a global average from those sites alone would underestimate the true global mean by ~0.8‰ (Figure 11). In fact, those models that fit the mean in the data best (#6, 8) indicate a ~0.8‰ higher global mean in the LGM compared to PI.

Sedimentary N-loss has a weak isotope effect that does not play a prominent direct role determining δ15N patterns in the model. Only areas directly above continental shelves with high sedimentary N-loss (e.g., Bering Sea) rates show a notable δ15N decrease due to its direct impact (Figure 9B). Its small effect on the δ15N makes it a difficult process to isotopically constrain, which creates high uncertainty on N cycling and inventory since it is the largest N-loss process in the global ocean.

The sensitivity simulations testing changes to sedimentary N-loss rates suggest that its most important isotope effect is an indirect one by stimulating changes to N2 fixation. N2 fixation rates are sensitive to N-loss because it enhances N-limitation, which in large part determines diazotroph's ecological niche. In LGMloSedNl, the reduced N2 fixation in response to sedimentary N-loss increased δ15N in the tropics and subtropics, which reproduced observations better in the western Tropical North Atlantic and Pacific (Figure 11E), but worse in Central Equatorial Pacific (Figure 11D) compared to LGMctl.


LGM Model Best Estimates and Uncertainty Ranges

Most LGM model simulations discussed here have similar global δ15N metrics. In general, the sensitivity experiments performed better in some regions but worse in others compared to LGMctl. Therefore, based on these objective global statistical metrics alone, one might conclude that most of our LGM models are similarly consistent with the observations. The only obvious exception is LGMloFe, which has a lower correlation coefficient than the other models, suggesting that enhanced iron fertilization played a crucial role in the glacial nitrogen cycle and its isotopic pattern, consistent with SS16's conclusions.

The simulations that perform among the best with both correlation coefficient (i.e., R > 0.44) and RMS error (i.e., error <1.55) metrics are LGMctl, LGMmidWCNl, and LGMloSedNl that largely determine our uncertainty ranges. However, some of the simulations perform better in regions that are representative of certain processes, which is also considered and discussed below. Overall, the model simulations struggle the most where the observations show small LGM−PI differences between ±1.5‰, whereas the stronger δ15N trends in ODZs and the Southern Ocean are better reproduced (Figure 10).

One of the main model biases in LGMctl is the too low δ15N in ODZs, presumably due to underestimated water column N-loss. The sensitivity simulations that assume higher water column N-loss reproduce those observations more closely. LGMhiWCNl, which predicts a 35% decrease in water column N-loss relative to PIctl, performs best in tropical ODZs but overestimates δ15N in the subarctic Pacific where water column N-loss also occurs in that simulation. The best simulation from SS16 that reproduces tropical ODZ slightly better predicts only a 17% decrease in water column N-loss. Therefore, we choose its global water column N-loss rate (57 Tg N yr−1) as our high-end estimate with our next best performing simulation in ODZs (LGMmidWCNl: 26 Tg N yr−1) as our low-end estimate, which yields a 17–62% decrease relative to PIctl.

Sedimentary N-loss is more difficult to constrain given its weak isotope effect. Our model predicts a lower global rate since sea level was reduced during the LGM. Since the sensitivity simulation LGMloSedNl results in improved model-data statistical metrics compared to LGMhiSedNl (Table 3, Figure 10), we use that simulations as our low-end range, with a value between LGMctl and LGMhiSedNl as our high-end estimate. This yields a 35–69% reduction of sedimentary N-loss during the LGMctl relative to PIctl.

The limited LGM observations interpreted as constraints on N2 fixation are all located in the western tropical North Pacific and Atlantic. Our LGMctl model predicts reduced N2 fixation in these regions relative to PIctl, in general agreement with previous interpretations (Ren et al., 2009, 2012). This occurs in the model in response to lower N-loss rates, which decrease the ecological nice for diazotrophy during the LGM relative to PI by reducing N-limitation. However, model LGMctl predicts an increase of N2 fixation in the tropical and subtropical southern hemisphere due to enhanced atmospheric Fe deposition. Unfortunately no observational constraints exist in these Fe-limited regions in the southern hemisphere, where N2 fixation has been hypothesized to occur increase during the LGM (Falkowski, 1997).

The most important sensitivity for global NO3- in our model is the assumption about changes to global PO43-, which fuels additional N2 fixation that increases global NO3-. The model simulation with a low increase to PO43- (LGMloPO43-) did not stimulate enough productivity to maintain sufficient water column N-loss to reproduce δ15N observations in ODZs, whereas simulations with higher global PO43- performed better with respect to data there. However, many circulation and biogeochemical processes determine the size of ODZs so it remains difficult to confidently conclude that increased PO43- was likely the reason for maintaining ODZs in the LGM, but it remains a possibility that cannot be excluded here. In the Central Equatorial Pacific, LGMloPO43- performed best due to its lower N2 fixation rates (Figure 11D), which suggests that the additional PO43- stimulated too much N2 fixation there when a higher PO43- inventory was assumed.

The best performing model sensitivity simulations predict a large range in changes to the global NO3- inventory (13–22%), which is most sensitive to assumptions for global PO43- inventory that has high uncertainties. If some unresolved physical process, e.g., increased diapycnal mixing (Schmittner et al., 2015), would be responsible for the δ15N bias in ODZs, this could decrease global NO3- similarly to the LGMhiWCNl sensitivity simulation (i.e., +27 Tg N yr−1 water column N-loss caused −2.2 mmol m−3 global NO3-). Applying this idealized linear relationship with our high-end estimate for water column N-loss reduction (57 Tg N yr−1) to LGMctl would reduce its global NO3- increase approximately in half because it underestimates water column N-loss. Thus, for our low-end uncertainty range for global NO3-, we apply this relative decrease to our simulation with lowest PO43- that predicts a +13% global NO3- increase relative to PIctl, which yields a low-end increase of global NO3- of ~6.5% based on our sensitivity simulations performed in this study.

Comparison with Previous Studies

Our model simulations produce lower oxygen levels in the glacial deep ocean below ~1,500 m compared to preindustrial (Figure 4) due to a combination of a weaker Atlantic meridional overturning circulation (Figure 1), which reduces deep ocean ventilation, and enhanced atmospheric Fe deposition over the Southern Ocean (Figure 2), which consumes oxygen during additional organic matter respiration. This large-scale model feature is consistent with proxy reconstructions (Jaccard and Galbraith, 2012). Since atmospheric Fe deposition was mainly increased over the Southern Ocean, enhanced export production stores more nutrients in Antarctic Bottom Water that fills the deep ocean. This causes less preformed nutrients transport to the tropics via Subantarctic Mode Waters, thereby reducing productivity and water column N-loss there, generally consistent with the scenario suggested by Costa et al. (2016).

The glacial upper ocean, on the other hand, was better oxygenated due to enhanced solubility from colder sea surface waters, which reduces the volume of oxygen deficient zones relative to preindustrial. The reduced oxygen deficient zones and exposed continental shelves due to lower sea level decreases water column and sedimentary N-loss in our model by 17–69% and 35–69%, respectively. Our uncertainty range agrees with Galbraith et al. (2013), although noting they examine the period between 15,000 and 8,000 years ago when the strongest deglacial changes occur.

Our simulations suggest that the ratio of sedimentary to water column total N-loss was higher in the LGM assuming no change to global δ15NO3-, in contrast to Galbraith et al. (2013) who suggest a constant ratio. This occurs in our model due to the dilution effect (Deutsch et al., 2004), which decreases the apparent isotope effect of N-loss (per mole N removed) as NO3- utilization increases in N-loss zones. Since in LGMctl N-loss rates were smaller due to diminished upper ocean ODZs and continental shelves, the dilution effect also decreased, thereby increasing the isotope effect of water column and sedimentary N-loss by 5.2 and 1.1‰, respectively, relative to PIctl. Therefore, additional N2 fixation stimulated by sedimentary N-loss would be required to balance global δ15NO3- assuming an oceanic N isotope budget in equilibrium. However, the observed LGM-PI δ15N decrease in ODZs do not support such a large decrease of water column N-loss (Figure 11C). This leads to the possibility that global δ15NO3- was increased during the LGM, which we suggest cannot be excluded with the current sparsity of observations because the true model global sedimentary δ15N average (“full model” in Figure 11A) is not always consistent with limited locations where observations exist (“data-masked model” in Figure 11A). Our simulations suggest some of this potential global δ15NO3- increase would accumulate in the Southern Ocean (Figures 8K,L) and be partially responsible for the sedimentary δ15N increase there (Figures 9K,L).

N2 fixation decreased during the LGM relative to preindustrial in our model in response to reduced N-loss rates following the interpretation by Ren et al. (2009), in contrast to the hypothesis by Falkowski (1997). Enhanced iron deposition during the LGM does not stimulate significant N2 fixation because deposition is estimated to mostly occur in high latitudes (Figure 2; Albani et al., 2014) where N2-fixers are limited by low temperatures in our model. However, model experiments testing the hypothesis of a higher PO43- inventory during the LGM (Wallmann et al., 2016) predict additional N2 fixation that further increases global NO3- inventory, marine productivity, and water column N-loss.

Our model estimated range for the LGM global NO3- increase (6.5–22%) is larger than the more idealized simulations in SS16 (4.5%). The main factors for this difference is that our simulations account for sea level effects on sedimentary N-loss and test higher PO43- inventories, which were neglected in SS16. SS16's best δ15N simulation predicted a 13% decrease to total N-loss, leading to only a 4.5% increase in the LGM global NO3- inventory. Since our model simulations account for more realistic LGM conditions that decrease total N-loss by ~50%, in general agreement with Galbraith et al. (2013), and reproduce the observations with a similar amount of skill, it is likely that SS16's estimate for glacial NO3- inventory (+4.5%) is underestimated.

Previous box modeling studies have estimated LGM N cycling. Deutsch et al. (2004) used a 4-box model constrained by δ15N observations in the eastern tropical North Pacific and account for sea level effects on sedimentary N-loss. They predict conservative changes to the LGM N cycle of ~30% decrease of water column and sedimentary N-loss with a global NO3- inventory that was increased by probably 10% or less than preindustrial. In an expanded 12 box model approach, Eugster et al. (2013) predict a larger global NO3- increase (~30–100%) that is driven by additional N2 fixation stimulated by enhanced atmospheric Fe deposition in the LGM relative to the deglaciation in that model.

While some regions in our model simulations suggest N2 fixation may increase in response to enhanced atmospheric Fe, it was never great enough to significantly increase the global NO3- inventory alone as suggested by Falkowski (1997), Eugster et al. (2013). On the global average, reduced N-loss caused lower N2 fixation that is supported by observations in the western tropical North Atlantic and Pacific. Although, we must note that no δ15N observations exist in the oligotrophic southern tropics and subtropics, where our simulations predict that N2 fixation would be most sensitive to enhanced atmospheric Fe deposition. Therefore, we cannot exclude the scenario by Eugster et al. (2013) of an even larger NO3- inventory and higher N2 fixation rates than predicted by our LGM simulations. However, they acknowledge that their sensitivity simulations with lower N2 fixation during the LGM performed only marginally worse in their model-data comparison, which would be consistent with our model simulations in this study.

Strategies for Future Model Development

Improved simulations may be possible in the future by using a fully interactive iron cycle. Our simplified iron cycle does not account for a higher quota for diazotrophs compared with other phytoplankton that could be another mechanism to further increase N2 fixation in the LGM. On the other hand, reduced sea level effects on dissolved iron concentrations are not accounted for here, which could compensate increase atmospheric deposition to some degree since sedimentary Fe flux is estimated to be a major source of Fe to the ocean (Dale et al., 2015). Our model assumes constant elemental stoichiometry of phytoplankton in the LGM and preindustrial. Some studies suggest that nitrogen to phosphorus quotas may have increased in the glacial ocean (Weber and Deutsch, 2012; Galbraith and Martiny, 2015), potentially enhancing N-limitation that could have stimulated additional N2 fixation, leading to a larger NO3- inventory and stronger biological carbon pump. We suggest future simulations should also test the biogeochemical sensitivity to different physical circulation states as it could be an underlying cause driving biogeochemical biases.


We have presented the first 3D marine nitrogen cycle model with realistic LGM boundary conditions that can be directly constrained by sedimentary δ15N records. Our model sensitivity simulations consistently predict large reductions in N-loss and N2 fixation rates, which together result in higher global NO3- inventories. The models that best reproduce global sedimentary δ15N observations predict that during the LGM N-loss was decreased by 17–62% in the water column, due to higher oxygen concentrations in the thermoclime, and by 35–69% in the sediments, due to lower sea level. We estimate that N2 fixation was decreased by 25–65% with a resulting increase to the NO3- inventory of 6.5–22%. Its uncertainty remains large given the few open ocean δ15N observations that provide constraints on N2 fixation, particularly in the southern hemisphere, and high uncertainties associated with changes to the global PO43- inventory. Our model-data δ15N analysis suggests that changes to circulation, sea surface oxygen solubility, and biological production from elevated marine iron and phosphate inventories all played potentially important roles determining the pattern and intensity of changes to dissolved oxygen and nitrogen cycling during the LGM. We suggest that increases to the oceanic NO3- inventory remains a viable candidate to explain part of the glacial-interglacial variations in atmospheric CO2 and should be considered in future studies.

Author Contributions

CS, AS, and AO designed the study and model experiments. JM calculated the LGM atmospheric iron fluxes and compiled the wind forcing. CS performed the model experiments, analysis, and wrote the paper with contributions from AS, AO, and JM.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.


We appreciate the NICOPP group for compiling the bulk δ15N database and everyone who made their data publicly available. We thank two reviewers for constructive comments that improved the manuscript. CS and AO acknowledge support by the SFB 754 project from the Deutsche Forschungsgemeinshaft and the PalMOD project from the BMBF. AS and JM were supported by National Science Foundation (#1634719). Model data used in this study are available at

Supplementary Material

The Supplementary Material for this article can be found online at:


Adkins, J. F., McIntyre, K., and Schrag, D. P. (2002). The salinity, temperature, and δ18O of the glacial deep Ocean. Science 298, 1769–1773. doi: 10.1126/science.1076252

PubMed Abstract | CrossRef Full Text | Google Scholar

Albani, S., Mahowald, N. M., Perry, A. T., Scanza, R. A., Zender, C. S., Heavens, N. G., et al. (2014). Improved dust representation in the community atmosphere model. J. Adv. Model. Earth Syst. 6, 541–570. doi: 10.1002/2013MS000279

CrossRef Full Text | Google Scholar

Altabet, M. A. (2007). Constraints on oceanic N balance/imbalance from sedimentary 15N records. Biogeosciences 4, 75–86. doi: 10.5194/bg-4-75-2007

CrossRef Full Text | Google Scholar

Altabet, M. A., and Francois, R. (1994). Sedimentary nitrogen isotopic ratio as a recorder for surface ocean nitrate utilization. Global Biogeochem. Cycles 8, 103–116. doi: 10.1029/93GB03396

CrossRef Full Text | Google Scholar

Altabet, M. A., Francois, R., Murray, D. W., and Prell, W. L. (1995). Climate-related variations in denitrification in the Arabian Sea from sediment 15N/14N ratios. Nature 373, 506–509. doi: 10.1038/373506a0

CrossRef Full Text | Google Scholar

Barcelose Ramos, J., Biswas, H., Schulz, K. G., Laroche, J., and Riebesell, U. (2007). Effect of rising atmospheric carbon dioxide on the marine nitrogen fixer Trichodesmium. Global Biogeochem. Cycles 21, GB2028. doi: 10.1029/2006gb002898

CrossRef Full Text | Google Scholar

Bianchi, D., Dunne, J. P., Sarmiento, J. L., and Galbraith, E. D. (2012). Data-based estimates of suboxia, denitrification, and N2O production in the ocean and their sensitivities to dissolved O2. Global Biogeochem. Cycles 26, GB2009. doi: 10.1029/2011GB004209

CrossRef Full Text | Google Scholar

Bohlen, L., Dale, A. W., and Wallmann, K. (2012). Simple transfer functions for calculating benthic fixed nitrogen losses and C:N:P regeneration ratios in global biogeochemical models. Global Biogeochem. Cycles 26, GB3029. doi: 10.1029/2011GB004198

CrossRef Full Text | Google Scholar

Braconnot, P., Harrison, S. P., Kageyama, M., Bartlein, P. J., Masson-Delmotte, V., Abe-Ouchi, A., et al. (2012). Evaluation of climate models using palaeoclimatic data. Nat. Clim. Chang. 2, 417–424. doi: 10.1038/nclimate1456

CrossRef Full Text | Google Scholar

Brandes, J. A., and Devol, A. H. (1997). Isotopic fractionation of oxygen and nitrogen in coastal marine sediments. Geochim. Cosmochim. Acta 61, 1793–1801. doi: 10.1016/S0016-7037(97)00041-0

CrossRef Full Text | Google Scholar

Brandes, J. A., and Devol, A. H. (2002). A global marine-fixed nitrogen isotopic budget: implications for Holocene nitrogen cycling. Global Biogeochem. Cycles 16, 67.1–67.14. doi: 10.1029/2001GB001856

CrossRef Full Text | Google Scholar

Broecker, W. S. (1982). Glacial to interglacial changes in ocean chemistry. Prog. Oceanogr. 11, 151–197. doi: 10.1016/0079-6611(82)90007-6

CrossRef Full Text | Google Scholar

Broecker, W. S., and Henderson, G. M. (1998). The sequence of events surrounding Termination II and their implications for the cause of glacial-interglacial CO2 changes. Paleoceanography 13, 352–364. doi: 10.1029/98PA00920

CrossRef Full Text | Google Scholar

Brovkin, V., Ganopolski, A., Archer, D., and Rahmstorf, S. (2007). Lowering of glacial atmospheric CO2 in response to changes in oceanic circulation and marine biogeochemistry. Paleoceanography 22:PA4202. doi: 10.1029/2006PA001380

CrossRef Full Text | Google Scholar

Casciotti, K. L. (2009). Inverse kinetic isotope fractionation during bacterial nitrite oxidation. Geochim. Cosmochim. Acta 73, 2061–2076. doi: 10.1016/j.gca.2008.12.022

CrossRef Full Text | Google Scholar

Casciotti, K. L. (2016). Nitrite isotopes as tracers of marine N cycle processes. Philos. Trans. R. Soc. A 374:20150295. doi: 10.1098/rsta.2015.0295

CrossRef Full Text | Google Scholar

Christensen, J. P., Murray, J. W., Devol, A. H., and Codispoti, L. A. (1987). Denitrification in continental shelf sediments has major impact on the oceanic nitrogen budget. Global Biogeochem. Cycles 1, 97–116. doi: 10.1029/GB001i002p00097

CrossRef Full Text | Google Scholar

Cline, J. D., and Kaplan, I. R. (1975). Isotopic fractionation of dissolved nitrate during denitrification in the eastern tropical north pacific ocean. Mar. Chem. 3, 271–299. doi: 10.1016/0304-4203(75)90009-2

CrossRef Full Text | Google Scholar

Codispoti, L. A. (2007). An oceanic fixed nitrogen sink exceeding 400 Tg N a−1 vs. the concept of homeostasis in the fixed-nitrogen inventory. Biogeosciences 4, 233–253. doi: 10.5194/bg-4-233-2007

CrossRef Full Text | Google Scholar

Costa, K. M., McManus, J. F., Anderson, R. F., Ren, H., Sigman, D. M., Winckler, G., et al. (2016). No iron fertilization in the equatorial Pacific Ocean during the last ice age. Nature 529, 519–522. doi: 10.1038/nature16453

PubMed Abstract | CrossRef Full Text | Google Scholar

Dale, A. W., Nickelsen, L., Scholz, F., Hensen, C., Oschlies, A., and Wallmann, K. (2015). A revised global estimate of dissolved iron fluxes from marine sediments. Global Biogeochem. Cycles 29, 691–707. doi: 10.1002/2014GB005017

CrossRef Full Text | Google Scholar

Dale, A. W., Sommer, S., Ryabenko, E., Noffke, A., Bohlen, L., Wallmann, K., et al. (2014). Benthic nitrogen fluxes and fractionation of nitrate in the Mauritanian oxygen minimum zone (Eastern Tropical North Atlantic). Geochim. Cosmochim. Acta 134, 234–256. doi: 10.1016/j.gca.2014.02.026

CrossRef Full Text | Google Scholar

Deutsch, C., Gruber, N., Key, R. M., Sarmiento, J. L., and Ganachaud, A. (2001). Decreasing marine biogenic calcification: a negative feedback on rising atmospheric pCO2. Global Biogeochem. Cycles 15, 483–506. doi: 10.1029/2000GB001291

CrossRef Full Text

Deutsch, C., Sigman, D. M., Thunell, R. C., Meckler, A. N., and Haug, G. H. (2004). Isotopic constraints on glacial/interglacial changes in the oceanic nitrogen budget. Global Biogeochem. Cycles 18, GB4012. doi: 10.1029/2003GB002189

CrossRef Full Text | Google Scholar

Devries, T., and Deutsch, C. (2014). Large-scale variations in the stoichiometry of marine organic matter respiration. Nat. Geosci. 7, 890–894. doi: 10.1038/ngeo2300

CrossRef Full Text | Google Scholar

Devries, T., Deutsch, C., Rafter, P. A., and Primeau, F. (2013). Marine denitrification rates determined from a global 3-D inverse model. Biogeosciences 10, 2481–2496. doi: 10.5194/bg-10-2481-2013

CrossRef Full Text | Google Scholar

Eugster, O., and Gruber, N. (2012). A probabilistic estimate of global marine N-fixation and denitrification. Global Biogeochem. Cycles 26, GB4013. doi: 10.1029/2012GB004300

CrossRef Full Text | Google Scholar

Eugster, O., Gruber, N., Deutsch, C., Jaccard, S. L., and Payne, M. R. (2013). The dynamics of the marine nitrogen cycle across the last deglaciation. Paleoceanography 28, 116–129. doi: 10.1002/palo.20020

CrossRef Full Text | Google Scholar

Falkowski, P. G. (1997). Evolution of the nitrogen cycle and its influence on the biological sequestration of CO2 in the ocean. Nature 387, 272–275. doi: 10.1038/387272a0

CrossRef Full Text | Google Scholar

Fischer, T., Banyte, D., Brandt, P., Dengler, M., Krahmann, G., Tanhua, T., et al. (2013). Diapycnal oxygen supply to the tropical North Atlantic oxygen minimum zone. Biogeosciences 10, 5079–5093. doi: 10.5194/bg-10-5079-2013

CrossRef Full Text | Google Scholar

Galbraith, E. D., Gnanadesikan, A., Dunne, J. P., and Hiscock, M. R. (2010). Regional impacts of iron-light colimitation in a global biogeochemical model. Biogeosciences 7, 1043–1064. doi: 10.5194/bg-7-1043-2010

CrossRef Full Text | Google Scholar

Galbraith, E. D., Kienast, M., Galbraith, E. D., Kienast, M., Albuquerque, A. L., Altabet, M. A., et al. (2013). The acceleration of oceanic denitrification during deglacial warming. Nat. Geosci. 6, 579–584. doi: 10.1038/ngeo1832

CrossRef Full Text | Google Scholar

Galbraith, E. D., and Martiny, A. C. (2015). A simple nutrient-dependence mechanism for predicting the stoichiometry of marine ecosystems. Proc. Natl. Acad. Sci. U.S.A. 112, 8199–8204. doi: 10.1073/pnas.1423917112

PubMed Abstract | CrossRef Full Text | Google Scholar

Garcia, H. E., Locarnini, R. A., Boyer, T. P., Antonov, J. I., Baranov, O. K., Zweng, M. M., et al. (2010a). “World Ocean Atlas 2009, Volume 3: dissolved oxygen, apparent oxygen utilization, and oxygen saturation,” in NOAA Atlas NESDIS 70, ed S. Levitus (Washington, DC: U.S. Government Printing Office), 344.

Google Scholar

Garcia, H. E., Locarnini, R. A., Boyer, T. P., Antonov, J. I., Zweng, M. M., Baranov, O. K., et al. (2010b). “World Ocean Atlas 2009, Volume 4: nutrients (phosphate, nitrate, silicate),” in NOAA Atlas NESDIS 71, ed S. Levitus (Washington, DC: U.S. Government Printing Office), 398.

Gent, P. R., and McWilliams, J. C. (1990). Isopycnal mixing in ocean circulation models. J. Phys. Oceanogr. 20, 150–155. doi: 10.1175/1520-0485(1990)020<0150:IMIOCM>2.0.CO;2

CrossRef Full Text | Google Scholar

Getzlaff, J., and Dietze, H. (2013). Effects of increased isopycnal diffusivity mimicking the unresolved equatorial intermediate current system in an earth system climate model. Geophys. Res. Lett. 40, 2166–2170. doi: 10.1002/grl.50419

CrossRef Full Text | Google Scholar

Gradoville, M. R., White, A. E., Böttjer, D., Church, M. J., and Letelier, R. M. (2014). Diversity trumps acidification: lack of evidence for carbon dioxide enhancement ofTrichodesmiumcommunity nitrogen or carbon fixation at Station ALOHA. Limnol. Oceanogr. 59, 645–659. doi: 10.4319/lo.2014.59.3.0645

CrossRef Full Text | Google Scholar

Granger, J., Prokopenko, M. G., Sigman, D. M., Mordy, C. W., Morse, Z. M., Morales, L. V., et al. (2011). Coupled nitrification-denitrification in sediment of the eastern Bering Sea shelf leads to 15N enrichment of fixed N in shelf waters. J. Geophys. Res. 116, C11006. doi: 10.1029/2010JC006751

CrossRef Full Text | Google Scholar

Grosskopf, T., and Laroche, J. (2012). Direct and indirect costs of dinitrogen fixation in Crocosphaera watsonii WH8501 and possible implications for the nitrogen cycle. Front. Microbiol. 3:236. doi: 10.3389/fmicb.2012.00236

PubMed Abstract | CrossRef Full Text | Google Scholar

Großkopf, T., Mohr, W., Baustian, T., Schunck, H., Gill, D., Kuypers, M. M. M., et al. (2012). Doubling of marine dinitrogen-fixation rates based on direct measurements. Nature 488, 361–364. doi: 10.1038/nature11338

PubMed Abstract | CrossRef Full Text | Google Scholar

Gruber, N. (2008). “Chapter 1 - The marine nitrogen cycle: overview and challenges,” in Nitrogen in the Marine Environment, 2nd Edn., eds D. G. Capone, D. A. Bronk, M. R. Mulholland and E. J. Carpenter (San Diego, CA: Academic Press), 1–50.

Google Scholar

Gruber, N., and Galloway, J. N. (2008). An Earth-system perspective of the global nitrogen cycle. Nature 451, 293–296. doi: 10.1038/nature06592

PubMed Abstract | CrossRef Full Text | Google Scholar

Hain, M. P., Sigman, D. M., and Haug, G. H. (2014). “8.18 - The biological pump in the past,” in Treatise on Geochemistry, 2nd Edn., eds K. K. Turekian and H. D. Holland (Oxford: Elsevier), 485–517.

Google Scholar

Hutchins, D. A., Fu, F. X., Zhang, Y., Warner, M. E., Feng, Y., Portune, K., et al. (2007). CO2 control of Trichodesmium N2 fixation, photosynthesis, growth rates, and elemental ratios: implications for past, present, and future ocean biogeochemistry. Limnol. Oceanogr. 52, 1293–1304. doi: 10.4319/lo.2007.52.4.1293

CrossRef Full Text | Google Scholar

Insua, T. L., Spivack, A. J., Graham, D., D'hondt, S., and Moran, K. (2014). Reconstruction of Pacific Ocean bottom water salinity during the Last Glacial Maximum. Geophys. Res. Lett. 41, 2914–2920. doi: 10.1002/2014GL059575

CrossRef Full Text | Google Scholar

Jaccard, S. L., and Galbraith, E. D. (2012). Large climate-driven changes of oceanic oxygen concentrations during the last deglaciation. Nat. Geosci. 5, 151–156. doi: 10.1038/ngeo1352

CrossRef Full Text | Google Scholar

Kalnay, E., Kanamitsu, M., Kistler, R., Collins, W., Deaven, D., Gandin, L., et al. (1996). The NCEP/NCAR 40-year reanalysis project. Bull. Am. Meteorol. Soc. 77, 437–471. doi: 10.1175/1520-0477(1996)077<0437:TNYRP>2.0.CO;2

CrossRef Full Text | Google Scholar

Karl, D., Michaels, A., Bergman, B., Capone, D., Carpenter, E., Letelier, R., et al. (2002). Dinitrogen fixation in the world's oceans. Biogeochemistry 57,58, 47–98. doi: 10.1023/A:1015798105851

CrossRef Full Text | Google Scholar

Key, R. M., Kozyr, A., Sabine, C. L., Lee, K., Wanninkhof, R., Bullister, J. L., et al. (2004). A global ocean carbon climatology: results from Global Data Analysis Project (GLODAP). Global Biogeochem. Cycles 18:GB4031. doi: 10.1029/2004GB002247

CrossRef Full Text | Google Scholar

Knapp, A. N., Sigman, D. M., Lipschultz, F., Kustka, A. B., and Capone, D. G. (2011). Interbasin isotopic correspondence between upper-ocean bulk DON and subsurface nitrate and its implications for marine nitrogen cycling. Global Biogeochem. Cycles 25, GB4004. doi: 10.1029/2010GB003878

CrossRef Full Text | Google Scholar

Large, W. G., Danabasoglu, G., McWilliams, J. C., Gent, P. R., and Bryan, F. O. (2001). Equatorial circulation of a Global Ocean climate model with anisotropic horizontal viscosity. J. Phys. Oceanogr. 31, 518–536. doi: 10.1175/1520-0485(2001)031<0518:ECOAGO>2.0.CO;2

CrossRef Full Text | Google Scholar

Law, C. S., Breitbarth, E., Hoffmann, L. J., McGraw, C. M., Langlois, R. J., Laroche, J., et al. (2012). No stimulation of nitrogen fixation by non-filamentous diazotrophs under elevated CO2 in the South Pacific. Glob. Chang. Biol. 18, 3004–3014. doi: 10.1111/j.1365-2486.2012.02777.x

CrossRef Full Text | Google Scholar

Lehmann, M. F., Sigman, D. M., McCorkle, D. C., Granger, J., Hoffmann, S., Cane, G., et al. (2007). The distribution of nitrate 15N/14N in marine sediments and the impact of benthic nitrogen loss on the isotopic composition of oceanic nitrate. Geochim. Cosmochim. Acta 71, 5384–5404. doi: 10.1016/j.gca.2007.07.025

CrossRef Full Text | Google Scholar

Luo, C., Mahowald, N., Bond, T., Chuang, P. Y., Artaxo, P., Siefert, R., et al. (2008). Combustion iron distribution and deposition. Global Biogeochem. Cycles 22:GB1012. doi: 10.1029/2007gb002964

CrossRef Full Text | Google Scholar

McElroy, M. B. (1983). Marine biological controls on atmospheric CO2 and climate. Nature 302, 328–329. doi: 10.1038/302328a0

CrossRef Full Text | Google Scholar

Meissner, K. J., Galbraith, E. D., and Völker, C. (2005). Denitrification under glacial and interglacial conditions: a physical approach. Paleoceanography 20, PA3001. doi: 10.1029/2004PA001083

CrossRef Full Text | Google Scholar

Middelburg, J. J., Soetaert, K., Herman, P. M. J., and Heip, C. H. R. (1996). Denitrification in marine sediments: a model study. Global Biogeochem. Cycles 10, 661–673. doi: 10.1029/96GB02562

CrossRef Full Text | Google Scholar

Minagawa, M., and Wada, E. (1984). Stepwise enrichment of 15N along food chains: further evidence and the relation between δ15N and animal age. Geochim. Cosmochim. Acta 48, 1135–1140. doi: 10.1016/0016-7037(84)90204-7

CrossRef Full Text | Google Scholar

Minagawa, M., and Wada, E. (1986). Nitrogen isotope ratios of red tide organisms in the East China Sea: a characterization of biological nitrogen fixation. Mar. Chem. 19, 245–259. doi: 10.1016/0304-4203(86)90026-5

CrossRef Full Text | Google Scholar

Mohr, W., Großkopf, T., Wallace, D. W. R., and Laroche, J. (2010). Methodological underestimation of oceanic nitrogen fixation rates. PLoS ONE 5:e12583. doi: 10.1371/journal.pone.0012583

PubMed Abstract | CrossRef Full Text | Google Scholar

Moore, C. M., Mills, M. M., Arrigo, K. R., Berman-Frank, I., Bopp, L., Boyd, P. W., et al. (2013). Processes and patterns of oceanic nutrient limitation. Nat. Geosci. 6, 701–710. doi: 10.1038/ngeo1765

CrossRef Full Text | Google Scholar

Muglia, J., and Schmittner, A. (2015). Wind stress increases glacial atlantic overturning in climate models. Geophys. Res. Lett. 42, 9862–9868. doi: 10.1002/2015GL064583

CrossRef Full Text

Mulholland, M. R., Ohki, K., and Capone, D. G. (2001). Nutrient controls on nitrogen uptake and metabolism by natural populations and cultures of Trichodesmium (Cyanobacteria). J. Phycol. 37, 1001–1009. doi: 10.1046/j.1529-8817.2001.00080.x

CrossRef Full Text | Google Scholar

Nickelsen, L., Keller, D. P., and Oschlies, A. (2015). A dynamic marine iron cycle module coupled to the University of Victoria Earth System Model: the Kiel Marine Biogeochemical Model 2 for UVic 2.9. Geosci. Model Dev. 8, 1357–1381. doi: 10.5194/gmd-8-1357-2015

CrossRef Full Text | Google Scholar

O'Neil, J. M. (1999). Grazer Interactions with Nitrogen-Fixing Marine Cyanobacteria: Adaptation for N-Acquisition? Monaco: Musée océanographique.

Google Scholar

Peltier, W. R. (2004). Global glacial isostacy and the surface of the ice-age Earth: the ICE-5G (VM2) model and GRACE. Annu. Rev. Earth Planet. Sci. 32, 111–149. doi: 10.1146/

CrossRef Full Text | Google Scholar

Ren, H., Sigman, D. M., Chen, M.-T., and Kao, S.-J. (2012). Elevated foraminifera-bound nitrogen isotopic composition during the last ice age in the South China Sea and its global and regional implications. Global Biogeochem. Cycles 26:GB1031. doi: 10.1029/2010GB004020

CrossRef Full Text | Google Scholar

Ren, H., Sigman, D. M., Meckler, A. N., Plessen, B., Robinson, R. S., Rosenthal, Y., et al. (2009). Foraminiferal isotope evidence of reduced nitrogen fixation in the ice age Atlantic Ocean. Science 323, 244–248. doi: 10.1126/science.1165787

PubMed Abstract | CrossRef Full Text | Google Scholar

Ren, H., Studer, A. S., Serno, S., Sigman, D. M., Winckler, G., Anderson, R. F., et al. (2015). Glacial-to-interglacial changes in nitrate supply and consumption in the subarctic North Pacific from microfossil-bound N isotopes at two trophic levels. Paleoceanography 30, 1217–1232. doi: 10.1002/2014PA002765

CrossRef Full Text | Google Scholar

Robinson, R. S., Kienast, M., Luiza Albuquerque, A., Altabet, M., Contreras, S., De Pol Holz, R., et al. (2012). A review of nitrogen isotopic alteration in marine sediments. Paleoceanography 27, PA4203. doi: 10.1029/2012PA002321

CrossRef Full Text | Google Scholar

Sarnthein, M., Schneider, B., and Grootes, P. M. (2013). Peak glacial 14C ventilation ages suggest major draw-down of carbon into the abyssal ocean. Clim. Past 9, 2595–2614. doi: 10.5194/cp-9-2595-2013

CrossRef Full Text | Google Scholar

Schmittner, A., and Egbert, G. D. (2014). An improved parameterization of tidal mixing for ocean models. Geosci. Model Dev. 7, 211–224. doi: 10.5194/gmd-7-211-2014

CrossRef Full Text | Google Scholar

Schmittner, A., Green, J. A. M., and Wilmes, S. B. (2015). Glacial ocean overturning intensified by tidal mixing in a global circulation model. Geophys. Res. Lett. 42, 4014–4022. doi: 10.1002/2015GL063561

CrossRef Full Text | Google Scholar

Schmittner, A., and Somes, C. J. (2016). Complementary constraints from carbon (13C) and nitrogen (15N) isotopes on the glacial ocean's soft-tissue biological pump. Paleoceanography 31, 669–693. doi: 10.1002/2015PA002905

CrossRef Full Text | Google Scholar

Sigman, D. M., Altabet, M. A., McCorkle, D. C., Francois, R., and Fischer, G. (2000). The δ15N of nitrate in the Southern Ocean: nitrogen cycling and circulation in the ocean interior. J. Geophys. Res. 105, 19599–19614. doi: 10.1029/2000JC000265

CrossRef Full Text | Google Scholar

Sigman, D. M., De Boer, A. M., and Haug, G. H. (2007). “Antarctic stratification, atmospheric water vapor, and heinrich events: a hypothesis for late pleistocene deglaciations,” in Ocean Circulation: Mechanisms and Impacts—Past and Future Changes of Meridional Overturning, eds A. Schmittner, J. C. H. Chiang, and S. R. Hemming (Washington, DC: American Geophysical Union), 335–349.

Google Scholar

Somes, C. J., and Oschlies, A. (2015). On the influence of “non-Redfield” dissolved organic nutrient dynamics on the spatial distribution of N2 fixation and the size of the marine fixed nitrogen inventory. Global Biogeochem. Cycles 29, 973–993. doi: 10.1002/2014GB005050

CrossRef Full Text | Google Scholar

Somes, C. J., Oschlies, A., and Schmittner, A. (2013). Isotopic constraints on the pre-industrial oceanic nitrogen budget. Biogeosciences 10, 5889–5910. doi: 10.5194/bg-10-5889-2013

CrossRef Full Text | Google Scholar

Somes, C. J., Schmittner, A., Galbraith, E. D., Lehmann, M. F., Altabet, M. A., Montoya, J. P., et al. (2010). Simulating the global distribution of nitrogen isotopes in the ocean. Global Biogeochem. Cycles 24, GB4019. doi: 10.1029/2009GB003767

CrossRef Full Text | Google Scholar

Tagliabue, A., Aumont, O., Death, R., Dunne, J. P., Dutkiewicz, S., Galbraith, E., et al. (2016). How well do global ocean biogeochemistry models simulate dissolved iron distributions? Global Biogeochem. Cycles 30, 149–174. doi: 10.1002/2015GB005289

CrossRef Full Text | Google Scholar

Tesdal, J. E., Galbraith, E. D., and Kienast, M. (2012). The marine sedimentary nitrogen isotope record. Biogeosci. Discuss. 9, 4067–4097. doi: 10.5194/bgd-9-4067-2012

CrossRef Full Text | Google Scholar

Tesdal, J. E., Galbraith, E. D., and Kienast, M. (2013). Nitrogen isotopes in bulk marine sediment: linking seafloor observations with subseafloor records. Biogeosciences 10, 101–118. doi: 10.5194/bg-10-101-2013

CrossRef Full Text | Google Scholar

Wada, E. (1980). “Nitrogen isotope fractionation and its significance in biogeochemical processes occurring in marine environments,” in Isotope Marine Chemistry, eds E. D. Golderg, Y. Horibe, and K. Saruhashi (Tokyo: Uchida-Rokakuho), 375–398.

Google Scholar

Wada, E., and Hattori, A. (1978). Nitrogen isotope effects in the assimilation of inorganic nitrogenous compounds by marine diatoms. Geomicrobiol. J. 1, 85–101. doi: 10.1080/01490457809377725

CrossRef Full Text | Google Scholar

Wallmann, K. (2010). Phosphorus imbalance in the global ocean? Global Biogeochem. Cycles 24, GB4030. doi: 10.1029/2009GB003643

CrossRef Full Text | Google Scholar

Wallmann, K., Schneider, B., and Sarnthein, M. (2016). Effects of eustatic sea-level change, ocean dynamics, and nutrient utilization on atmospheric pCO2 and seawater composition over the last 130 000 years: a model study. Clim. Past 12, 339–375. doi: 10.5194/cp-12-339-2016

CrossRef Full Text | Google Scholar

Weaver, A. J., Eby, M., Wiebe, E. C., Bitz, C. M., Duffy, P. B., Ewen, T. L., et al. (2001). The UVic earth system climate model: model description, climatology, and applications to past, present and future climates. Atmos. Ocean 39, 361–428. doi: 10.1080/07055900.2001.9649686

CrossRef Full Text | Google Scholar

Weber, T., and Deutsch, C. (2012). Oceanic nitrogen reservoir regulated by plankton diversity and ocean circulation. Nature 489, 419–422. doi: 10.1038/nature11357

PubMed Abstract | CrossRef Full Text | Google Scholar


Model Changes from Previous Versions

Since only minor model formulation changes have been made from SS16, in this appendix we will document only these changes and refer to previous studies for a full description.

Physical Model

We applied a lower vertical mixing rate of 1.5 × 10−5 m2 s−1 in tropical subsurface between 185 and 565 m based on microstructure and tracer release experiments in the tropical Atlantic (Fischer et al., 2013). Our chosen value is near their high-end mixing rate uncertainty range (0.8–1.4 × 10−5 m2 s−1), whereas SS16 assumed a constant value 0.3 × 10−5 m2 s−1 everywhere in the global ocean. Since this lower vertical mixing rate reduced AMOC, a higher tidal mixing e-folding depth (1,000 m compared to 500 m in SS16) was chosen to strengthen AMOC and better reproduce Δ14C (Figure S2).

Marine Ecosystem-Biogeochemical Model

The fast recycling parameter was reduced to 0.001 from 0.015 d−1 (Table A1) to better reproduce global net primary productivity, surface NO3- and surface DON in this study, which were overestimated in SS16. The grazing preference of diazotrophs by zooplankton was reduced (0.04 from 0.1), but the quadratic mortality term, which represents specialist diazotrophs grazers not explicitly modeled, was increased so global N2 fixation rate remained unchanged. With these changes to diazotrophs loss term parameters, the quadratic mortality term is 20% higher than general zooplankton grazing when averaged over the global ocean, whereas in previous versions the loss term from general zooplankton was about an order of magnitude higher than the specialist grazing loss term or not included. We reduced the fractionation factor of water column N-loss to 20‰ (from 25‰ in SS16) since studies have shown that NO2 oxidation reduces the net isotope effect of water column N-loss (Casciotti, 2009, 2016), which is an isotope effect that is not explicitly included in our model. To maintain global preindustrial δ15NO3- at its observed level, the fractionation factor of sedimentary N-loss was increased to 6‰ (from 4‰ in SS16).


Table A1. Marine Ecosystem-Biogeochemistry Parameters.

Keywords: marine nitrogen cycle, nitrogen isotopes, last glacial maximum, oceanic nitrogen budget, biological carbon pump

Citation: Somes CJ, Schmittner A, Muglia J and Oschlies A (2017) A Three-Dimensional Model of the Marine Nitrogen Cycle during the Last Glacial Maximum Constrained by Sedimentary Isotopes. Front. Mar. Sci. 4:108. doi: 10.3389/fmars.2017.00108

Received: 31 October 2016; Accepted: 03 April 2017;
Published: 09 May 2017.

Edited by:

Phyllis Lam, University of Southampton, UK

Reviewed by:

Haojia Ren, National Taiwan University, Taiwan
Moritz Felix Lehmann, University of Basel, Switzerland

Copyright © 2017 Somes, Schmittner, Muglia and Oschlies. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Christopher J. Somes,