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

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 the global response of nitrogen fixation, which decreases by 25‒65% in the model during the LGM, 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.

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 CO 2 during the LGM.
Keywords: marine nitrogen cycle, nitrogen isotopes, last glacial maximum, oceanic nitrogen budget, biological carbon pump

INTRODUCTION
The global nitrate (NO − 3 ) 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, NO − 3 is the main limiting nutrient throughout the tropical and subtropical oceans (Moore et al., 2013) because Nloss processes cause its depletion relative to other macronutrients such as phosphate . Therefore, significant changes to the global NO − 3 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 N 2 fixation by specialized cyanobacteria (diazotrophs) capable of converting dinitrogen gas to ammonium to meet their nitrogen requirements for growth. N 2 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 N 2 -fixing cellular compartment anoxic, as well as additional structural iron requirements of the N 2 -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, O 2 <∼10 mmol m −3 ) where NO − 3 replaces O 2 as the electron acceptor during organic matter respiration. Under these suboxic conditions, forms of fixed nitrogen (NH 4 , NO 2 , NO − 3 ) are converted to N 2 O and N 2 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 N 2 fixation and N-loss mainly determines the preindustrial global NO − 3 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 N 2 fixation, leading to an anthropogenic ocean that could be rapidly losing N (Codispoti, 2007). However, methodological issues with historical N 2 fixation measurements suggest a significant underestimation of N 2 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.
δ 15 N 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). δ 15 N records are typically interpreted qualitatively as changes to nutrient utilization in modern High Nitrate Low Chlorophyll (HNLC) regions, water column N-loss, or N 2 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 δ 15 N 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, threedimensional models have not been used so far to quantify these effects and their impacts on the global nitrogen inventory.
Effects of glacial conditions on N 2 fixation are more uncertain. Falkowski (1997) suggests that increased atmospheric Fe deposition to the glacial ocean stimulated additional N 2 fixation that could have increased the global NO − 3 inventory. Broecker (1982), and more recently (Wallmann et al., 2016) propose that a higher PO 3− 4 inventory, due to reduced loss on continental shelves, could have enhanced N 2 fixation during the LGM since PO 3− 4 is diazotrophs' other main limiting nutrient. On the other hand, Ren et al. (2009Ren et al. ( , 2012 interpret δ 15 N records from the western tropical North Atlantic and Pacific as reduced N 2 fixation, in response to lower N-loss in the glacial ocean. Some specific N 2 -fixing Trichodesmium species may also be limited by CO 2 concentrations of surface waters (Barcelos e Ramos et al., 2007;Hutchins et al., 2007), which could have reduced N 2 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 CO 2 on N 2 fixation.
Previous box modeling studies have produced large differences in their estimates of the global NO − 3 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 CO 2 , whereas the higher estimates could lead to enhanced atmospheric CO 2 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 CO 2 (McElroy, 1983;Christensen et al., 1987;Falkowski, 1997;Broecker and Henderson, 1998), possible effects of changes in the nitrogen inventory on glacialinterglacial CO 2 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 N 2 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 PO 3− 4 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).

Marine Biogeochemical Model
The marine ecosystem-biogeochemical model coupled within the ocean circulation includes 2 nutrients in the inorganic (NO − 3 and PO 3− 4 ) and organic (DON and DOP) phases, 2 phytoplankton (ordinary and N 2 -fixing diazotrophs), zooplankton, sinking detritus, as well as dissolved O 2 , dissolved inorganic carbon, alkalinity, and 14 C (Somes and Oschlies, 2015). Iron limitation is calculated using monthly surface dissolved iron fields prescribed from the BLING model . 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 14 C, O 2, and NO − 3 are reproduced ( Figure S2).

N 2 Fixation
Diazotrophs grow slower than ordinary phytoplankton due to extra energetic demands of N 2 -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 NO − 3 , but still contain sufficient P and Fe (e.g., water with low NO − 3 from denitrification and high iron from atmospheric Fe deposition). They will consume NO − 3 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 PO 3− 4 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).

N-Loss
Water column N-loss occurs when dissolved oxygen becomes depleted in poorly ventilated ODZs. NO − 3 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 O 2 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 O 2 . 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 accelSedN−l = [(1 − F SGS ) × Γ SedN−l + 1] × SedN−l, where accelSedN−l is the accelerated sedimentary N-loss rate used in the model, F SGS is the fraction of the model grid covered by the subgrid-scale scheme ( Figure S1), SedN−l is the true sedimentary N-loss function from Bohlen et al. (2012), and Γ SedN−l 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).
The patterns and rates of N 2 fixation and N-loss are generally consistent with its observational biogeochemical indicators N * and sedimentary δ 15 N (Figure 2). Water column N-loss produces low N * and high δ 15 N 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, N 2 fixation conversely generates high N * and low δ 15 N. Low NO − 3 utilization due to iron limitation in HNLC regions of the Southern Ocean, North Pacific, and eastern equatorial Pacific drive low δ 15 N, whereas high utilization in the oligotrophic causes higher δ 15 N in the absence of N 2 fixation.
Global mean δ 15 NO − 3 is determined by global rates and isotope effects of N 2 fixation, and water column and sedimentary N-loss (Brandes and Devol, 2002;Altabet, 2007). Sedimentary Nloss rates and its isotope effect represent the largest uncertainties for determining global mean δ 15 NO − 3 (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 δ 15 NO − 3 near its observed ∼5‰ value given the model predicted rates of N 2 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.
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 , 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). 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 . 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 dFe LGM = (atmFeflx LGM /atmFeflx PI ) B × dFe PI . This pragmatic approach is used due to high uncertainties associated with 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). 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).

PO 3− 4 Inventory
We test a hypothesis that the global PO 3− 4 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 PO 3− 4 inventory is conserved in each model simulation, we implement this hypothesis by imposing the estimated + ∼15% increase in the global PO 3− 4 inventory from Wallmann et al. (2016) in LGMctl (

N-Loss Rates
We perform additional simulations to investigate how different water column and sedimentary N-loss rates affect LGM nitrogen cycling and δ 15 N distributions. In the model, water column Nloss can be increased pragmatically by increasing the elemental ratio of oxygen consumed per organic nitrogen remineralized (R O: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 R O : 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 R O : 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.

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 14 C 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 14 C 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 O 2 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 O 2 inventory ( Table 2). This is similar to the 20.8% reduction resulting from SS16's best δ 15 N simulation, who attribute ∼7% of this decrease to physics and ∼14% to changes in Fe fertilization.
Higher O 2 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 NO − 3 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.
N 2 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 N 2 fixation due to enhanced atmospheric Fe deposition ( Figure 5A). In the LGMctl model, global N 2 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 NO − 3 inventory (+22%). Note that DOP utilization by diazotrophs does not play a significant role stimulating N 2 fixation in LGMctl because nitrogen limitation is reduced due to the larger NO − 3 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 NO − 3 is thus consumed at the surface and stored in the deep ocean (Figures 7A,B), where it can remain for millennia.
Enhanced export production at high latitudes reduces the amount of preformed surface NO − 3 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 N 2 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 NO − 3 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 N 2 fixation in our model. However, on the global average, these areas of additional N 2 fixation are much smaller than areas in which N 2 fixation decreased in model LGMctl relative to PIctl due to reduced N-loss ( Figure 5A).
Water column δ 15 NO − 3 is strongly influenced by nitrogen cycling changes in response to enhanced atmospheric Fe deposition (Figure 8). Enhanced (reduced) surface NO − 3 utilization increases (decreases) surface δ 15 NO − 3 in the Southern Ocean and Subarctic Pacific in LGMhiFe (LGMloFe). The minor increase (decrease) to N 2 fixation rates in the southern tropics and subtropics in LGMhiFe (LGMloFe) substantially decreases (increases) δ 15 NO − 3 . This large isotopic signal occurs because these systems switch from being driven by high NO − 3 utilization, which causes extremely 15 N-enriched surface NO − 3 , to a N 2 fixation isotopic regime that introduces 15 N-depleted nitrogen to the surface ocean. This model prediction suggests that all relevant processes should be considered when interpreting δ 15 N and the change to the isotopic signal may not be linearly representative of the rate change to one specific nitrogen cycling process alone.

Global PO 3− 4 Inventory Sensitivity
The model simulations testing a hypothesized increase to the global PO 3− 4 inventory predict larger NO − 3 inventories due to enhanced N 2 fixation (Table 2, Figures 6, 7G-J). PO 3− 4 is the major limiting macronutrient for N 2 -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 PO 3− 4 stimulates additional growth of diazotrophs via N 2 fixation that is most pronounced in the tropical ocean. The additional organic matter production via enhanced N 2 fixation and remineralization at depth stimulates additional water column N-loss, which balances the enhanced N 2 fixation at a new steady-state with a larger global NO − 3 (Figure 6) and reduced O 2 inventory (Figures 4I,J). A 10% increase in the PO 3− 4 inventory thus causes an 8% increase in the NO − 3 inventory (LGMctl-LGMloPO 3− 4 ). However, a further 10% increase in the PO 3− 4 inventory (LGMhiPO 3− 4 -LGMctl) leads to diminishing returns in terms of NO − 3 (4%) due to the further expansion of ODZs and N-loss. The additional N 2 fixation decreases surface δ 15 NO − 3 in the tropics and subtropics, whereas the additional water column N-loss increases subsurface δ 15 NO − 3 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 NO − 3 in the eastern tropical Pacific, North Indian, and subarctic Pacific (Figures 7K-N). This significantly increases subsurface δ 15 NO − 3 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 NO − 3 (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) N 2 fixation to balance additional (less) sedimentary N-loss in LGMhiSedNl (LGMloSedNl). This decreases (increases) surface δ 15 NO − 3 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  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 δ 15 N 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, δ 15 N 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 N 2 fixation, and higher values in high latitude HNLC regions (i.e., Southern Ocean, western subarctic Pacific) indicating enhanced surface NO − 3 utilization (Figures 9-11). These sedimentary δ 15 N records are compared to the δ 15 N values of sinking particulate organic nitrogen (PON) in bottom waters above the seafloor, which we refer to as sedimentary δ 15 N in the model below. Our model does not account for terrestrial sources of δ 15 N, 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.

LGM Control Simulation
Model simulation LGMctl, which includes enhanced atmospheric Fe deposition and assumes a +15% increase to the global PO 3− 4 inventory, reproduces the major trends in the LGM reconstructions (R = 0.45, Table 3; Figures 9-11). It reproduces the largest increase in δ 15 N 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 δ 15 N is too strong relative to observations ( Figure 11C).
The change to global δ 15 NO − 3 ( Table 2) and sedimentary δ 15 N ( 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 δ 15 NO − 3 (Somes et al., 2013). This relatively small δ 15 NO − 3 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 NO − 3 utilization decreases in the smaller ODZs. This effect decreases the isotopic signature of δ 15 N 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 N 2 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 δ 15 NO − 3 in the LGMctl steady-state solution compared to PIctl.
The processes that lead to slightly higher deep ocean δ 15 NO −

in
LGMctl are mainly controlled at high latitudes. In the North Atlantic, reduced N 2 fixation elevates surface δ 15 NO − 3 that eventually subducts in the Atlantic meridional overturning circulation and increases deep ocean δ 15 NO − 3 (Figure 8A). In the Southern Ocean, enhanced atmospheric Fe deposition increases NO − 3 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 δ 15 NO − 3 as well. Since most of the increased NO − 3 inventory in LGMctl accumulates in these old, sluggish deep waters, our modeling suggests that high latitude process also play a role determining global mean δ 15 NO − 3 , not only N 2 fixation and water column N-loss in the tropics.

Lower atmospheric Fe deposition in LGMloFe compared to
LGMctl reduces NO − 3 utilization and sedimentary δ 15 N 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 δ 15 N 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 NO − 3 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 NO − 3 utilization and hence sedimentary δ 15 N in this region. In our high-end simulation for atmospheric Fe deposition (LGMhiFe), enhanced N 2 fixation introducing additional 15 Ndepleted 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 NO − 3 utilization or increased N 2 fixation was the most important process driving the observed decreased LGM δ 15 N trend.
Our sensitivity model experiments in the Southern Hemisphere show a balance between lower δ 15 N due to enhanced N 2 fixation in the subtropics and higher δ 15 N from enhanced surface NO − 3 utilization toward high latitudes. Our model predicts that increased Fe could stimulate additional N 2 fixation at the poleward edge of the southern subtropics as far south as ∼40 • S, which counteracts with the isotope  effect of enhanced NO − 3 utilization that increases δ 15 N at higher latitudes (Figures 8C-F, 9C-F). As NO − 3 decreases with enhanced atmospheric Fe deposition, the additional 15 N-depleted nitrogen via N 2 fixation contributes to a larger fraction of total sinking PON, which prevents a higher δ 15 N increase in LGMhiFe compared to the δ 15 N decrease in LGMloFe (Figure 11B). The absence of open ocean δ 15 N observations in the tropical and subtropical southern hemisphere makes it uncertain to validate this increase to N 2 fixation there.

Global PO 3− 4 Inventory Sensitivity
In the model, a larger global PO 3− 4 inventory stimulates increases to N 2 fixation, water column N-loss, and the global NO − 3 inventory (Figure 6). This results in increased sedimentary δ 15 N near ODZs and decreased δ 15 N outside of them in the more oligotrophic ocean where additional N 2 fixation occurs ( Figure 9J). The magnitude of change to δ 15 N in the model (Figure 10) is relatively small because N 2 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 15 N-depleted nitrogen via fixation roughly compensates the extra 15 N-enriched NO − 3 via enhanced water column N-loss rates, preventing a larger change to the large scale δ 15 N signal (Figures 9G-J).
Despite large changes to the global rates of these processes, changes to simulated δ 15 N at the location of observations only moderately affect the model-data statistical metrics (Table 3, Figure 10). The higher water column N-loss in LGMhiPO 3− 4 better reproduces δ 15 N in tropical ODZs ( Figure 11C). However, the additional water column N-loss in the subarctic Pacific overestimates δ 15 N in that region. The model experiment LGMloPO 3− 4 performs best in the Central Equatorial Pacific (Figure 11E), which suggests that the additional PO 3− 4 stimulates too much N 2 fixation there in simulations LGMctl and LGMhiPO 3− 4 .

N-Loss Rate Sensitivity
Since one of the main model deficiencies in LGMctl is underestimated δ 15 N 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 δ 15 N 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 δ 15 N in the Central Equatorial Pacific and Subarctic Pacific, causing those regions to overestimate δ 15 N compared to the observations, demonstrating how water column N-loss influences δ 15 N across the entire tropical Pacific and beyond (Figures 8K-N, 9K-N,P).
Higher water column N-loss rates increase global δ 15 NO − 3 (Figures 8K-N), but this does not impact sedimentary δ 15 N everywhere equally (Figures 9K-N). Effects are largest in regions with upwelling of the increased deep ocean δ 15 NO − 3 , e.g., tropical Pacific, subarctic North Pacific, and Southern Ocean. Our simulations suggest that if a significant increase in global δ 15 NO − 3 occurred during the LGM, this could significantly impact δ 15 N across the Southern Ocean, which are typically interpreted as changes to surface NO − 3 utilization alone. Given the sparsity of the global dataset, we cannot exclude a significant increase in global δ 15 NO − 3 during the LGM because its largescale impact affects sedimentary δ 15 N in regions that already have an increased δ 15 N 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 δ 15 N patterns in the model. Only areas directly above continental shelves with high sedimentary N-loss (e.g., Bering Sea) rates show a notable δ 15 N decrease due to its direct impact ( Figure 9B). Its small effect on the δ 15 N 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 Nloss rates suggest that its most important isotope effect is an indirect one by stimulating changes to N 2 fixation. N 2 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 N 2 fixation in response to sedimentary N-loss increased δ 15 N 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 δ 15 N 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 δ 15 N trends in ODZs and the Southern Ocean are better reproduced (Figure 10).

One of the main model biases in
LGMctl is the too low δ 15 N in ODZs, presumably due to underestimated water column Nloss. 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 δ 15 N 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 Nloss 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 modeldata 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 N 2 fixation are all located in the western tropical North Pacific and Atlantic. Our LGMctl model predicts reduced N 2 fixation in these regions relative to PIctl, in general agreement with previous interpretations (Ren et al., 2009(Ren et al., , 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 N 2 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 N 2 fixation has been hypothesized to occur increase during the LGM (Falkowski, 1997).
The most important sensitivity for global NO − 3 in our model is the assumption about changes to global PO 3− 4 , which fuels additional N 2 fixation that increases global NO − 3 . The model simulation with a low increase to PO 3− 4 (LGMloPO 3− 4 ) did not stimulate enough productivity to maintain sufficient water column N-loss to reproduce δ 15 N observations in ODZs, whereas simulations with higher global PO 3− 4 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 PO 3− 4 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, LGMloPO 3− 4 performed best due to its lower N 2 fixation rates (Figure 11D), which suggests that the additional PO 3− 4 stimulated too much N 2 fixation there when a higher PO 3− 4 inventory was assumed. The best performing model sensitivity simulations predict a large range in changes to the global NO − 3 inventory (13-22%), which is most sensitive to assumptions for global PO 3− 4 inventory that has high uncertainties. If some unresolved physical process, e.g., increased diapycnal mixing , would be responsible for the δ 15 N bias in ODZs, this could decrease global NO − 3 similarly to the LGMhiWCNl sensitivity simulation (i.e., +27 Tg N yr −1 water column N-loss caused −2.2 mmol m −3 global NO − 3 ). 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 NO − 3 increase approximately in half because it underestimates water column Nloss. Thus, for our low-end uncertainty range for global NO − 3 , we apply this relative decrease to our simulation with lowest PO 3− 4 that predicts a +13% global NO − 3 increase relative to PIctl, which yields a low-end increase of global NO − 3 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 δ 15 NO − 3 , 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 NO − 3 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 N 2 fixation stimulated by sedimentary N-loss would be required to balance global δ 15 NO − 3 assuming an oceanic N isotope budget in equilibrium. However, the observed LGM-PI δ 15 N decrease in ODZs do not support such a large decrease of water column N-loss ( Figure 11C). This leads to the possibility that global δ 15 NO − 3 was increased during the LGM, which we suggest cannot be excluded with the current sparsity of observations because the true model global sedimentary δ 15 N 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 δ 15 NO − 3 increase would accumulate in the Southern Ocean (Figures 8K,L) and be partially responsible for the sedimentary δ 15 N increase there (Figures 9K,L).
N 2 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 N 2 fixation because deposition is estimated to mostly occur in high latitudes (Figure 2; Albani et al., 2014) where N 2 -fixers are limited by low temperatures in our model. However, model experiments testing the hypothesis of a higher PO 3− 4 inventory during the LGM (Wallmann et al., 2016) predict additional N 2 fixation that further increases global NO − 3 inventory, marine productivity, and water column N-loss. Our model estimated range for the LGM global NO − 3 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 PO 3− 4 inventories, which were neglected in SS16. SS16's best δ 15 N simulation predicted a 13% decrease to total N-loss, leading to only a 4.5% increase in the LGM global NO − 3 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 NO − 3 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 δ 15 N 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 NO − 3 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 NO − 3 increase (∼30-100%) that is driven by additional N 2 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 N 2 fixation may increase in response to enhanced atmospheric Fe, it was never great enough to significantly increase the global NO − 3 inventory alone as suggested by Falkowski (1997), Eugster et al. (2013). On the global average, reduced N-loss caused lower N 2 fixation that is supported by observations in the western tropical North Atlantic and Pacific. Although, we must note that no δ 15 N observations exist in the oligotrophic southern tropics and subtropics, where our simulations predict that N 2 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 NO − 3 inventory and higher N 2 fixation rates than predicted by our LGM simulations. However, they acknowledge that their sensitivity simulations with lower N 2 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 N 2 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 N 2 fixation, leading to a larger NO − 3 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.

CONCLUSIONS
We have presented the first 3D marine nitrogen cycle model with realistic LGM boundary conditions that can be directly constrained by sedimentary δ 15 N records. Our model sensitivity simulations consistently predict large reductions in N-loss and N 2 fixation rates, which together result in higher global NO − 3 inventories. The models that best reproduce global sedimentary δ 15 N 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 N 2 fixation was decreased by 25-65% with a resulting increase to the NO − 3 inventory of 6.5-22%. Its uncertainty remains large given the few open ocean δ 15 N observations that provide constraints on N 2 fixation, particularly in the southern hemisphere, and high uncertainties associated with changes to the global PO 3− 4 inventory. Our model-data δ 15 N 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 NO − 3 inventory remains a viable candidate to explain part of the glacial-interglacial variations in atmospheric CO 2 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.

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 m 2 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 m 2 s −1 ), whereas SS16 assumed a constant value 0.3 × 10 −5 m 2 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 14 C ( 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 NO − 3 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 N 2 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 NO 2 oxidation reduces the net isotope effect of water column N-loss (Casciotti, 2009(Casciotti, , 2016, which is an isotope effect that is not explicitly included in our model. To maintain global preindustrial δ 15 NO − 3 at its observed level, the fractionation factor of sedimentary N-loss was increased to 6‰ (from 4‰ in SS16).