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).

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 N₂-fixing diazotrophs), zooplankton, sinking detritus, as well as dissolved O₂, dissolved inorganic carbon, alkalinity, and Δ¹⁴C (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 Δ¹⁴C, O₂, and NO3- are reproduced (Figure S2).

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 (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⁻⁴ 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 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).

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.

We perform additional simulations to investigate how different water column and sedimentary N-loss rates affect LGM nitrogen cycling and δ¹⁵N distributions. In the model, water column N-loss 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⁻¹ (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.

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 Δ¹⁴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 Δ¹⁴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).

In LGMctl, oxygen concentrations are higher above ~1,500 m and lower below (Figures 4A,B) relative to PIctl, which is consistent with O₂ 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₂ inventory (Table 2). This is similar to the 20.8% reduction resulting from SS16's best δ¹⁵N simulation, who attribute ~7% of this decrease to physics and ~14% to changes in Fe fertilization.

Higher O₂ 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.

N₂ 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₂ fixation due to enhanced atmospheric Fe deposition (Figure 5A). In the LGMctl model, global N₂ 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 N₂ 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.

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 9–11). It reproduces the largest increase in δ¹⁵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 δ¹⁵N is too strong relative to observations (Figure 11C).

The change to global δ¹⁵NO3- (Table 2) and sedimentary δ¹⁵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 δ¹⁵NO3- (Somes et al., 2013). This relatively small δ¹⁵NO3- 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 δ¹⁵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₂ 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 δ¹⁵NO3- in the LGMctl steady-state solution compared to PIctl.

The processes that lead to slightly higher deep ocean δ¹⁵NO3- in LGMctl are mainly controlled at high latitudes. In the North Atlantic, reduced N₂ fixation elevates surface δ¹⁵NO3- that eventually subducts in the Atlantic meridional overturning circulation and increases deep ocean δ¹⁵NO3- (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 δ¹⁵NO3- 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 δ¹⁵NO3-, not only N₂ fixation and water column N-loss in the tropics.

Lower atmospheric Fe deposition in LGMloFe compared to LGMctl reduces NO3- utilization and sedimentary δ¹⁵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 δ¹⁵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 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 δ¹⁵N in this region. In our high-end simulation for atmospheric Fe deposition (LGMhiFe), enhanced N₂ fixation introducing additional ¹⁵N-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 N₂ fixation was the most important process driving the observed decreased LGM δ¹⁵N trend.

Our sensitivity model experiments in the Southern Hemisphere show a balance between lower δ¹⁵N due to enhanced N₂ fixation in the subtropics and higher δ¹⁵N from enhanced surface NO3- utilization toward high latitudes. Our model predicts that increased Fe could stimulate additional N₂ 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 δ¹⁵N at higher latitudes (Figures 8C–F, 9C–F). As NO3- decreases with enhanced atmospheric Fe deposition, the additional ¹⁵N-depleted nitrogen via N₂ fixation contributes to a larger fraction of total sinking PON, which prevents a higher δ¹⁵N increase in LGMhiFe compared to the δ¹⁵N decrease in LGMloFe (Figure 11B). The absence of open ocean δ¹⁵N observations in the tropical and subtropical southern hemisphere makes it uncertain to validate this increase to N₂ fixation there.

In the model, a larger global PO43- inventory stimulates increases to N₂ fixation, water column N-loss, and the global NO3- inventory (Figure 6). This results in increased sedimentary δ¹⁵N near ODZs and decreased δ¹⁵N outside of them in the more oligotrophic ocean where additional N₂ fixation occurs (Figure 9J). The magnitude of change to δ¹⁵N in the model (Figure 10) is relatively small because N₂ 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 ¹⁵N-depleted nitrogen via fixation roughly compensates the extra ¹⁵N-enriched NO3- via enhanced water column N-loss rates, preventing a larger change to the large scale δ¹⁵N signal (Figures 9G–J).

Despite large changes to the global rates of these processes, changes to simulated δ¹⁵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 LGMhiPO43- better reproduces δ¹⁵N in tropical ODZs (Figure 11C). However, the additional water column N-loss in the subarctic Pacific overestimates δ¹⁵N 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 N₂ fixation there in simulations LGMctl and LGMhiPO43-.

Since one of the main model deficiencies in LGMctl is underestimated δ¹⁵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 δ¹⁵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 δ¹⁵N in the Central Equatorial Pacific and Subarctic Pacific, causing those regions to overestimate δ¹⁵N compared to the observations, demonstrating how water column N-loss influences δ¹⁵N across the entire tropical Pacific and beyond (Figures 8K–N, 9K–N,P).

Higher water column N-loss rates increase global δ¹⁵NO3- (Figures 8K–N), but this does not impact sedimentary δ¹⁵N everywhere equally (Figures 9K–N). Effects are largest in regions with upwelling of the increased deep ocean δ¹⁵NO3-, e.g., tropical Pacific, subarctic North Pacific, and Southern Ocean. Our simulations suggest that if a significant increase in global δ¹⁵NO3- occurred during the LGM, this could significantly impact δ¹⁵N 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 δ¹⁵NO3- during the LGM because its large-scale impact affects sedimentary δ¹⁵N in regions that already have an increased δ¹⁵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 δ¹⁵N patterns in the model. Only areas directly above continental shelves with high sedimentary N-loss (e.g., Bering Sea) rates show a notable δ¹⁵N decrease due to its direct impact (Figure 9B). Its small effect on the δ¹⁵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 N-loss rates suggest that its most important isotope effect is an indirect one by stimulating changes to N₂ fixation. N₂ 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₂ fixation in response to sedimentary N-loss increased δ¹⁵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.

Most LGM model simulations discussed here have similar global δ¹⁵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 δ¹⁵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 δ¹⁵N 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 δ¹⁵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 N-loss rate (57 Tg N yr⁻¹) as our high-end estimate with our next best performing simulation in ODZs (LGMmidWCNl: 26 Tg N yr⁻¹) 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 N₂ fixation are all located in the western tropical North Pacific and Atlantic. Our LGMctl model predicts reduced N₂ 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 N₂ 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₂ 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 N₂ 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 δ¹⁵N 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 N₂ fixation rates (Figure 11D), which suggests that the additional PO43- stimulated too much N₂ 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 δ¹⁵N bias in ODZs, this could decrease global NO3- similarly to the LGMhiWCNl sensitivity simulation (i.e., +27 Tg N yr⁻¹ water column N-loss caused −2.2 mmol m⁻³ global NO3-). Applying this idealized linear relationship with our high-end estimate for water column N-loss reduction (57 Tg N yr⁻¹) 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.

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 δ¹⁵NO3-, 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 N₂ fixation stimulated by sedimentary N-loss would be required to balance global δ¹⁵NO3- assuming an oceanic N isotope budget in equilibrium. However, the observed LGM-PI δ¹⁵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 δ¹⁵NO3- was increased during the LGM, which we suggest cannot be excluded with the current sparsity of observations because the true model global sedimentary δ¹⁵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 δ¹⁵NO3- increase would accumulate in the Southern Ocean (Figures 8K,L) and be partially responsible for the sedimentary δ¹⁵N increase there (Figures 9K,L).

N₂ 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₂ fixation because deposition is estimated to mostly occur in high latitudes (Figure 2; Albani et al., 2014) where N₂-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 N₂ 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 δ¹⁵N 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 δ¹⁵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 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 N₂ 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₂ 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 N₂ fixation that is supported by observations in the western tropical North Atlantic and Pacific. Although, we must note that no δ¹⁵N observations exist in the oligotrophic southern tropics and subtropics, where our simulations predict that N₂ 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 N₂ fixation rates than predicted by our LGM simulations. However, they acknowledge that their sensitivity simulations with lower N₂ fixation during the LGM performed only marginally worse in their model-data comparison, which would be consistent with our model simulations in this study.

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₂ 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₂ 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.

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.