Effects of eustatic sea-level change, ocean dynamics, and nutrient utilization on atmospheric p CO 2 and seawater composition over the last 130 000 years: a model study

. We have developed and employed an Earth system model to explore the forcings of atmospheric p CO 2 change and the chemical and isotopic evolution of seawater over the last glacial cycle. Concentrations of dissolved phosphorus (DP), reactive nitrogen, molecular oxygen, dissolved inorganic carbon (DIC), total alkalinity (TA), 13 C-DIC, and 14 C-DIC were calculated for 24 ocean boxes. The bi-directional water ﬂuxes between these model boxes were derived from a 3-D circulation ﬁeld of the modern ocean (Opa 8.2, NEMO) and tuned such that tracer distributions calculated by the box model were consistent with observational data from the modern ocean. To model the last 130 kyr, we employed records of past changes in sea-level, ocean circulation, and dust deposition. According to the model, about half of the glacial p CO 2 drawdown may be attributed to marine regressions. The glacial sea-level low-stands implied steepened ocean margins, a reduced burial of particulate organic carbon, phosphorus, and neritic carbonate at the


Introduction
The discussion of mechanisms that might be responsible for the glacial to interglacial change in the atmosphere's CO 2 content is focused on the ocean (Broecker, 1982b).Everincreasing evidence suggests that CO 2 sequestered in the glacial ocean was rapidly released into the atmosphere at glacial terminations (Schmitt et al., 2012).The fast decline in dust-bound iron deposition in the Southern Ocean (Martin, 1990;Martınez-Garcia et al., 2014) and upwelling pulses in the Southern Ocean (Anderson et al., 2009) and North Pacific (Rae et al., 2014) may have induced the stepwise pCO 2 rise documented in the deglacial ice-core record (Marcott et al., 2014).The preceding CO 2 uptake in the glacial ocean may be attributed to enhanced export production, elevated  (Eakins and Sharman, 2012).The ocean margin at 0-100 m water depth is indicated for the modern ocean (red area) and for the LGM when eustatic sea-level was lowered by 120 m (blue area).The global ocean retreated into steeper terrain during the glacial marine regression.The seafloor areas covered by shallow waters were reduced by this steepening of ocean margins.
seawater alkalinity, and changes in ocean dynamics.The biological pump was probably intensified by iron fertilization (Martin, 1990) and the coeval expansion of nitrate (Deutsch et al., 2004) and phosphate (Broecker, 1982b) stocks in the glacial ocean, while seawater alkalinity may have been enhanced by the demise of neritic carbonate formation (Berger, 1982;Opdyke and Walker, 1992;Kleypas, 1997).The sequestration of atmospheric pCO 2 in the glacial ocean may have been further promoted by the glacial shoaling of the Meridional Overturning Circulation (MOC) in the Atlantic (Duplessy et al., 1988;Sarnthein et al., 1994), a possible increase in Southern Ocean stratification (Toggweiler, 1999), a prolonged residence time of surface waters in the Southern Ocean providing more time for the biota to draw down nutrients and CO 2 (Watson et al., 2015) and a global decline in MOC intensity and deep ocean ventilation (Sarnthein et al., 2013).
In part, a glacial rise in nitrate, phosphate, and alkalinity concentrations, which may have contributed significantly to the drawdown of atmospheric pCO 2 , can be explained by eustatic sea-level fall (Wallmann, 2014).It led to a retreat of ocean margins to steeper terrains that reduced the seafloor area located in shallow waters (Fig. 1).The standing stocks of carbon and nutrients in the glacial ocean may have been significantly enhanced by the marine regression and the decrease in shallow margin area since major removal fluxesnamely accumulation of neritic carbonate, benthic denitrification, burial of particulate organic carbon (POC) and phosphorus (P) -depend on the extent of seafloor located in shallow waters.Various authors and Earth system models con-sidered the glacial decrease in shelf carbonate burial as a major driver of ocean chemistry and atmospheric pCO 2 change (Berger, 1982;Opdyke and Walker, 1992;Brovkin et al., 2012;Ganopolski et al., 1998), since neritic carbonates contribute ≥ 50 % to the carbonate accumulation at the global seafloor (Milliman and Droxler, 1996;Kleypas, 1997;Berelson et al., 2007).However, the effects of sea-level change on POC and nutrient cycling are largely ignored in these stateof-the-art models even though > 50 % of the global benthic denitrification and burial of marine POC and P occur in shelf and upper slope environments (Berner, 1982;Bohlen et al., 2012;Wallmann, 2010).
Against this background, our contribution aims to explore and quantify the effects of sea-level change, ocean dynamics, and nutrient utilization on seawater composition and atmospheric pCO 2 over the last glacial cycle.We use a simple Earth system box model to simulate both chemical and isotopic changes in seawater composition and employ isotope data (δ 13 C) to constrain changes in ocean dynamics and deep ocean ventilation.Atmospheric pCO 2 and the distributions of dissolved oxygen, carbonate, and radiocarbon in the glacial ocean serve as key prognostic model variables.They are compared with independent proxy data to address the following specific questions: What fraction of the glacial pCO 2 drawdown can be ascribed to eustatic sea-level fall (Wallmann, 2014)?To which degree do global 14 C data sets assembled by Sarnthein et al. (2013Sarnthein et al. ( , 2015) ) actually reflect our concepts and ideas on glacial and deglacial ocean circulation and carbon cycling?Do marine 14 C data really form a quantitative proxy of DIC in the glacial deep ocean as proposed by Sarnthein et al. (2013)?To avoid circular reasoning, dissolved oxygen, carbonate, and radiocarbon distributions calculated for the glacial ocean were not used to parametrize our model.These distribution patterns and the atmospheric pCO 2 values calculated in the model are non-trivial consequences of interactions between the various model components and thus are employed to validate the model performance.
Eustatic sea-level change was applied as major model forcing (Fig. 3).Changes in global ocean volume, salinity, depositional area at continental margins, and exposed shelf area were derived from the sea-level record (Stanford et al., (Waelbroeck et al., 2002;Stanford et al., 2011); (b) global ocean volume as calculated from eustatic sea-level and ocean bathymetry data (Eakins and Sharman, 2012); (c) salinity of global mean seawater as calculated from global ocean volume; (d) global burial rate of neritic carbonate as calculated from seafloor area at 0-50 m water depth (Kleypas, 1997;Wallmann, 2014); (ef) seafloor area at 0-100 m and 100-2000 m water depth calculated from sea-level and ocean bathymetry data (Eakins and Sharman, 2012); (g) exposed shelf area calculated from sea-level and ocean bathymetry data (Eakins and Sharman, 2012); (h) global rate of POC weathering calculated from exposed shelf area (Wallmann, 2014).
2011; Waelbroeck et al., 2002) and the hypsographic curve (Eakins and Sharman, 2012).The burial rate of neritic carbonates was reduced during marine regressions in proportion to the decrease in seafloor area available for the growth of tropical reefs and carbonate platforms (Appendix B).We tested the degree to which the decrease in seafloor area at 0-100 and 100-2000 m water depth during glacial sea-level low-stands affected benthic denitrification and the burial of organic carbon and marine phosphorus on the continental shelf and slope (Appendix B) while we assumed that carbonate, P, and POC weathering were promoted by the exposure of shelf sediments (Appendix A).
The comprehensive geological database on benthic foraminiferal δ 13 C (Oliver et al., 2010;Sarnthein et al., 1994) was employed to constrain water fluxes for the Last Glacial Maximum (LGM).Mean δ 13 C-DIC values were calculated for those ocean boxes where sufficient δ 13 C data were available and compared to model results.Water fluxes were varied until the Holocene-LGM differences in δ 13 C generated by the model were consistent with the differences recorded in foraminifera (Table A5).The tuning was done using full transient runs with all forcings applied.The resulting fluxes are shown in Fig. 2. The southward water flux from the Atlantic into the Southern Ocean was relocated from deep (2000-4000 m) to intermediate waters (100-2000 m) to mimic the shoaling of the Atlantic Meridional Overturning Circulation (AMOC) which is inferred not only from δ 13 C data but also from various other proxy records (Curry and Oppo, 2005;Piotrowski et al., 2005;Roberts et al., 2010).Bottom water fluxes from the Southern Ocean into the Atlantic were enhanced during the LGM while the northward flow of surface water was reduced.The overall water exchange between the Atlantic and Southern Ocean was maintained constant at 15.4 Sv.The bottom and deep water exchange between the Southern Ocean and the Indo-Pacific was reduced by 5 Sv to reproduce the δ 13 C data (Table A5).
The eustatic sea-level curve (Fig. 3a) was applied to change ocean dynamics continuously over time so as to define water fluxes over the full model period (Fig. 4a-c).Thus we assumed that AMOC shoaled gradually during the transition from interglacial to full glacial conditions (130-21 ka) while the horizontal exchange flux of intermediate waters between the Southern Ocean and Tropical Indo-Pacific was enhanced over the glacial to mimic the ventilation of tropical oxygen minimum zones (OMZs) observed in various proxy records (Altabet et al., 1995;Jaccard and Galbraith, 2012).Additional rapid changes were implemented for the deglacial period (Fig. 4a and d).NADW formation was strongly reduced during Heinrich Event 1 (H1) and the Younger Dryas (McManus et al., 2004) while upwelling pulses were prescribed in the Southern Ocean during H1 and the Bølling-Allerød (Anderson et al., 2009;Skinner et al., 2010) and in the North Pacific during H1 (Rae et al., 2014).Timing and intensity of these deglacial upwelling/ventilation events were varied until 14 C-DIC values calculated for deep and surface water boxes were consistent with 14 C values recorded in benthic and pelagic foraminifera (Fig. A2).In addition, the deglacial ice-core record of pCO 2 (Marcott et al., 2014) and biogenic opal accumulation rates (Anderson et al., 2009) were employed to constrain the timing and intensity of upwelling pulses in the Southern Ocean.The iron accumulation record from site ODP 1090 was used to constrain changes in nutrient utilization in the Southern Ocean (Martınez-Garcia et al., 2014) assuming that the increase in iron accumulation observed at this site directly translates into an increase in the efficiency of nutrient utilization (Fig. 4e).To calculate realistic marine isotope trends, the changing isotopic com-  (Schmitt et al., 2012) while the solid line defines the values applied in the model.For > 24 kyr BP, where data are not available, the δ 13 C-CO 2 value is set to −6.4 ‰; (g) 14 C value of atmospheric CO 2 , dots indicate values reconstructed from the geological record (Reimer et al., 2013) while the solid line defines the values applied in the model.For > 50 kyr BP, where data are not available, the atmospheric 14 C-CO 2 is assumed to correspond to the pre-anthropogenic modern value (0 ‰). positions of atmospheric CO 2 were set to the values documented in the geological record (Fig. 4f-g).Atmospheric δ 13 C-CO 2 was not calculated as a prognostic model variable because a biased vertical δ 13 C-DIC gradient in the Southern Ocean impeded the simulation of 13 C-CO 2 fluxes across the ocean/atmosphere interface (Appendix A, Sect.A8, Fig. A3).
The major limitations of our simple box model are (i) very low spatial resolution, (ii) water fluxes between model boxes that are not derived from internal model dynamics, (iii) terrestrial inventories of POC in vegetation and soil that are kept constant over the model period.
As a consequence of low resolution, OMZs are not resolved by our model since the entire Indo-Pacific intermediate water at 100-2000 m water depth is pooled in a single ocean box.In our model, we prescribe a constant rate of pelagic denitrification since we are not able to resolve OMZs.Rates of benthic denitrification and P burial are only moderately affected by the lack of OMZs since the area where OMZs impinge the seafloor only amounts to 1 % of the global seafloor (Bohlen et al., 2012).
The modern water fluxes applied to our box model are based on a dynamically consistent circulation field.As explained in Appendix A, these fluxes were modified to obtain tracer distributions that are consistent with observations in the pre-industrial modern ocean (Fig. A1).Glacial and deglacial changes in ocean circulation were not derived from ocean models but from δ 13 C records (Table A5) and additional geochemical observations (Fig. A2).Models with explicit ocean dynamics are superior to any kind of physically unconstrained box model if they generate results that are consistent with observations.However, the physically constrained Earth system models that we are aware of are not yet able to reproduce as many tracer and proxy data as our box model.Our paper shows that shelf and sea-level effects help to explain a wide range of findings (Sect.3) and we think that the outputs of physically better-constrained models may improve in the case these effects are included in the model architecture.It is not our intention to promote box modelling per se as the method of choice.Rather, we hope that the concepts and ideas advanced in our paper may stimulate the community and help to further enhance cutting-edge Earth system models with explicit ocean dynamics (Tschumi et al., 2011;Menviel et al., 2012;Brovkin et al., 2012;Roth et al., 2014;Lambert et al., 2015).
For many decades it was widely assumed that the modern terrestrial carbon pool exceeds the glacial pool by hundreds of Gt because of biosphere regrowth after the glacial termination (Köhler and Fischer, 2004).This concept was initially developed to explain reduced δ 13 C values in glacial seawater (Shackleton, 1977).However, the latest assessment of terrestrial carbon pools indicates that the sum of the modern stocks does not exceed the LGM stock (Brovkin and Ganopolski, 2015).Moreover, this new view on terrestrial carbon cycling suggests a deglacial decline in total carbon stocks since the carbon release from high-latitude areas (melting permafrost soils and soils exposed by the retreat of glacial ice sheets) exceeded the carbon uptake by biosphere regrowth and peat accumulation.Our model explains deglacial and Holocene pCO 2 dynamics and the low glacial δ 13 C values by marine processes and sea-level change only.However, we acknowledge that terrestrial processes that are neglected in our model  (Monnin et al., 2001;Petit et al., 1999;Monnin et al., 2004) may have played a role especially over the Holocene even though more work needs to be done to constrain the sign and magnitude of terrestrial effects on δ 13 C and pCO 2 dynamics.
The major new component included in our Earth system model is a comprehensive formulation of shelf processes and sea-level effects.Appendix B explains in detail how sea-level change affects fluxes at continental margins and how these effects are considered in the model.The full model code is available from the first author on request.

Results and discussion
The model was run over a period of 130 kyr to simulate the behaviour of the global system over one full glacial cycle.changes in ocean circulation, and nutrient utilization, that is the full model forcing as defined in Figs. 3 and 4. Additional simulations were performed to better understand the controls on atmospheric pCO 2 and the chemical and isotopic composition of seawater.Simulation STD-CC was run with constant circulation -that is, all water fluxes were maintained at the Holocene level (upper panel of Fig. 2, Table A3) over the full model period -whereas simulation STD-CC-CN was performed with the Holocene circulation field and constant nutrient utilization.

Atmospheric pCO 2
The pCO 2 trend recorded in ice cores (Monnin et al., 2001(Monnin et al., , 2004;;Marcott et al., 2014;Petit et al., 1999) was well reproduced by the standard simulation (Figs.5a and 7a).Over the last interglacial, simulated pCO 2 increased from an initial value of 280 ppmv at 130 ka to 285 ppmv at 120 ka (Fig. 5a).This increase was accompanied by a decline in nutrient concentrations (Fig. 5f-g) and export production (Fig. 6d), supported by the high sea-level stand promoting burial of phosphorus and benthic denitrification in continental margin sediments.Over the subsequent glacial period, simulated pCO 2 declined due to sea-level fall, enhanced nutrient utilization in the Southern Ocean, and the decline in deep ocean ventilation until 21 ka, when a pCO 2 minimum of 190 ppmv was reached (Fig. 7b).The simulated glacial pCO 2 drawdown was discontinuous, marked by several steps and turning points (Fig. 5a).Major minima in atmospheric pCO 2 occurred at 90 ka (220 ppmv) and 65 ka (198 ppmv).Both of them are well documented in the ice-core record (Fig. 5a) and accompanied by maxima in nutrient utilization (Fig. 4e) and minima in sea-level (Fig. 3a).Sea-level fall and nutrient uti-lization, thus, may have driven most of the glacial pCO 2 decline.Moreover, they may have induced major turning points in the glacial pCO 2 record.At constant ocean circulation and nutrient utilization (simulation STD-CC-CN), simulated pCO 2 declined to a LGM value of 234 ppmv (Fig. 5a).Additional simulations based on the standard simulation STD helped us to specify the driving forces for this decline (Table 1).To study their effect on pCO 2 we suppressed the temporal changes of individual variables.A first simulation test was based on the assumption of constant modern sea surface temperatures (SSTs).It showed that the glacial decline in global mean SST by ca. 2 • C (Schmittner et al., 2011) induced a pCO 2 decline by 16 ppmv since the solubility of CO 2 in surface waters was enhanced under low temperatures (compare rows 1 and 2 in Table 1).In a second simulation, salinity was set constant, while the other model parameters varied as defined in the STD simulation.Accordingly, the peak glacial increase in salinity induced a relative atmospheric pCO 2 rise by 5 ppmv by lowering the solubility of CO 2 in surface waters (Table 1).In a third test, both salinity and the volume of the ocean boxes were kept constant over time.Changes in these parameters induced an LGM pCO 2 rise by 13 ppmv (Table 1), illustrating that the contraction of the ocean volume during glacial sea-level low-stands reduced the ocean's capacity to sequester atmospheric CO 2 .In summary, the model runs confirmed previous estimates (Broecker, 1982b) showing that the net effect of SST, volume, and salinity changes on glacial pCO 2 is small (decrease by 3 ppmv).Thus, other processes need to be invoked to explain the large glacial drawdown of atmospheric CO 2 simulated by model run STD-CC-CN.Changes in the flux of dissolved phosphorus (DP) exert large effects on pCO 2 , since DP is the ultimate limiting nutrient of the model ocean (Menviel et al., 2012).Neglecting the glacial increase in the weathering of P-bearing solids is raising the LGM pCO 2 value by 50 ppmv (Table 1, rows 5 vs. 1).Most of the P released during chemical weathering originates from apatite, a mineral equally occurring in all rock types (sedimentary, magmatic, and metamorphic).Thus, we assume that the P weathering rate is proportional to the total weathering rate -that is the sum of carbonate, POC, and silicate weathering (Wallmann, 2014).During the glacial, total weathering increased due to the weathering of exposed shelf CaCO 3 and POC (Munhoven, 2002;Wallmann, 2014).This rise led to the increase in P weathering simulated in the model.However, a further simulation shows that the overall pCO 2 change induced by chemical weathering of silicate, POC, CaCO 3 , and P is small (decrease by 3 ppmv) because the glacial CO 2 drawdown induced by P and CaCO 3 weathering was largely compensated by the CO 2 release induced by POC weathering (Table 1, STD run with constant rates of chemical weathering).Applying very high molar C : P ratios for POM in shelf sediments (ca.200), it was previously calculated that shelf weathering resulted in a net increase rather than decrease in atmospheric pCO 2 (Ushie and Matsumoto, 2012).However, most of the phosphorus in shelf sediments and riverine particles is not organic but bound in other reactive, inorganic phases such as carbonate-fluoro-apatite (Berner and Rao, 1994) which release DP when exposed to weathering (Ruttenberg, 1992;Ruttenberg and Berner, 1993).Thus, the global mean atomic ratio of POC over reactive P in shelf sediments is lower than the C : P ratio of marine organic matter (Baturin, 2007;Wallmann, 2010).Hence, the glacial weathering of shelf sediments induced a small drop rather than a rise in LGM pCO 2 (Table 1).
By contrast, a stronger effect results from testing the glacial decrease in depositional areas at continental margins, as revealed by a simulation that ignores the glacial decline in P burial and reveals a glacial pCO 2 rise by 73 ppmv (Table 1) with respect to the standard case due to the decline in DP concentration and export production.An additional simulation with constant depositional area for POC burial resulted in a pronounced drawdown of both atmospheric pCO 2  C4).The contour plots shown here and in Figs. 9 and 11 are based on concentrations calculated in simulation STD for each of the 24 ocean boxes (indicated as grid points).
(by 61 ppmv) and DIC since POC burial at continental margins served as a major sink for CO 2 and DIC in the model system (Table 1).Changes in the burial of neritic carbonates had a less drastic effect on atmospheric pCO 2 (change by 10 ppmv, Table 1) since they were mitigated by carbonate compensation at the deep-sea floor.Thus, the response of the model system to sea-level change was dominated by changes in the burial of P and POC at continental margins.The glacial drop in atmospheric pCO 2 that was induced by a decline in P burial and a corresponding rise in export production was moderated by a coeval drop in POC burial at continental margins.The overall effect was an enhanced sequestration of CO 2 in the glacial deep ocean accompanied by a decrease in sedimentary POC accumulation.Accordingly, most of the glacial pCO 2 decline in simulation STD-CC-CN was driven by the glacial steepening of ocean margins and the resulting expansion of the DP inventory.This conclusion is consistent with the results of previous experiments conducted with more evolved Earth system models showing a strong pCO 2 drawdown in response to an increase in the oceanic phosphate inventory (Tschumi et al., 2011;Menviel et al., 2012).
The glacial pCO 2 value dropped by 31 ppmv (from 234 to 203 ppmv) upon enhanced nutrient utilization (difference between simulations STD-CC-CN and STD-CC).This decrease was amplified by the glacial sea-level fall since the nutrient reservoir that was unlocked by the enhanced utilization in the glacial Southern Ocean was enlarged as a result of glacial marine regression.The remaining portion of the interglacial-to-peak glacial pCO 2 drop by 13 ppmv down to the final LGM value of 190 ppmv was induced by ocean dynamics (difference between simulation STD and STD-CC).Atmospheric pCO 2 rose by 10 ppmv, when all water fluxes between the Atlantic and the Southern Ocean were maintained at their Holocene level over the entire model period (simulation STD with constant AMOC, Table 1).The glacial AMOC shoaling (Fig. 4a-b) thus contributed 10 ppmv to the simulated LGM decline.This effect can be attributed to enhanced CO 2 storage in the deep Atlantic (> 2000 m water depth) which was less ventilated under glacial conditions, since the formation of northern deep waters was greatly diminished and replaced by southern-source waters enriched in DIC.Thus additional DIC was stored in the glacial deep ocean (Ganopolski et al., 2010;Skinner, 2009;Sarnthein et al., 2013).The glacial decrease in water fluxes between the deep Southern Ocean and Tropical Indo-Pacific applied in the model (Fig. 2) likewise supports further sequestration and storage of CO 2 in the deep ocean and the glacial drawdown of atmospheric pCO 2 .The circulation changes employed to simulate LGM conditions (Fig. 2 C4).See legend of Fig. 8 for further information.
by 210 years may be compared to the 600-year increase reconstructed from benthic radiocarbon data (Sarnthein et al., 2013).We suggest that the difference between these two estimates is related to the elevated production rate of radiocarbon in the glacial atmosphere and changes in carbon cycling affecting the marine radiocarbon budget (Sect.3.4).
A stepwise increase in pCO 2 was simulated over the deglaciation (Fig. 7).The first step occurred from 18.5 to 16.3 ka when the simulated pCO 2 rose rapidly from 193 to 220 ppmv.A second step followed at 15.9-14.1 ka with a pCO 2 rise from 222 to 244 ppmv, and a third step at 13.0-10.8ka with a strong increase from 243 to 272 ppmv.The first step was driven by the rapid drop in nutrient utilization at the glacial termination (Fig. 4e) and the ventilation of intermediate and deep water masses in the North Pacific during H1 (Figs. 4d and A2).The second and third steps were driven by the Southern Ocean where CO 2 was released into the atmosphere due to the abrupt decline in stratification and the further decrease in nutrient utilization (Fig. 4d-e).In major parts of the Southern Ocean these steps coincide with maxima in opal accumulation indicating enhanced upwelling (Anderson et al., 2009).The pCO 2 drawdown from 10 to 8 ka reflects a recovery of the ocean system from the antecedent ventilation pulse in the Southern Ocean centred at 11.5 ka (Fig. 4d).According to our model, the ventilation pulse removed CO 2 from the ocean interior, enhanced the O 2 content of the deep ocean, and diminished the vertical DIC and O 2 gradients.The subsequent restoration of the vertical DIC gradient induced the pCO 2 decline observed in the model from 10 to 8 ka.The biological pump needed about 2 kyr to reestablish the vertical DIC gradient due to the large inventory of DIC in the deep ocean.Sea-level change was not uncovered as a major driver for the rapid deglacial pCO 2 rise since the long residence times of DP (13 kyr) and TA (77 kyr) inhibit fast inventory changes in the global ocean (Menviel et al., 2012;Wallmann, 2014).By contrast, the pCO 2 increase over the Holocene (8-0 ka), which is closing the glacial cycle, may have been driven by a high sea-level stand inducing a gradual and slow decline in marine DP and TA inventories.
In our model about equal portions of the extreme deglacial carbon flux are triggered by changes in circulation and nutrient utilization in the Southern Ocean (Fig. 7).However, this model outcome is not well constrained because the lack of both proxy data and physical process understanding impedes an unequivocal determination of the magnitude of these changes in the Southern Ocean.

Dissolved nutrients and oxygen
The standing stock of DP in the global ocean rose under glacial conditions since P burial was diminished by the decrease in depositional area located at shallow water depths (Fig. 6h) while chemical weathering was promoted by the exposure of shelf sediments (Fig. 6g).Most of the glacial DP rise found in our simulations was induced by the glacial steepening of ocean margins reducing the burial of P in margin sediments (Table 1, Fig. 1, Appendix B).Vice versa, enhanced utilization in the Southern Ocean induced a strong decrease in glacial DP stocks since more DP was taken up by phytoplankton to be drawn down and finally buried in marine sediments (simulations STD-CC vs. STD-CC-CN, Figs.5f, 6d, h), while the glacial DP stock was largely restored by changes in ocean dynamics (simulations STD vs. STC-CC, Fig. 5f) separating the large nutrient pool in the deep ocean from the surface layer.The spatial distribution of DP in the global ocean reflects the export of POM by the biological pump and ocean circulation.The overall pattern, that is a strong vertical gradient between depleted surface waters and enriched deep water masses and a significant horizontal gradient between the deep North Atlantic and North Pacific, was maintained over the glacial cycle (Fig. 8).However, the vertical DP gradient was amplified over the LGM due to enhanced utilization and the decrease in deep ocean ventilation.Reactive P accumulation rates in marine sediments can be used to validate our model results.A global compilation of these data confirmed that P accumulation in shelf sediments decreased drastically under glacial conditions (Tamburini and Föllmi, 2009).The resulting decline in global P burial induced an increase in the glacial DP inventory by 17-40 % (Tamburini and Föllmi, 2009) as predicted by our model.Cd / Ca ratios in LGM sediments from the Atlantic Ocean (Boyle and Keigwin, 1982) and δ 13 C records (Duplessy et al., 1988;Sarnthein et al., 1994;Oliver et al., 2010) suggest a steepening of the vertical DP gradients broadly consistent with our model results.
The dissolved oxygen (DO) content of the global ocean decreased under glacial conditions due to the decline in deep ocean ventilation and increase in export production (Fig. 5h).In contrast, it recovered and peaked over the deglaciation since ocean ventilation was enhanced in the Southern Ocean and the North Pacific.The spatial distribution of DO changed significantly under LGM conditions (Fig. 8).Concentrations declined at > 2000 m water depth due to the decrease in ocean ventilation and increase in export production while glacial cooling induced a small DO rise in surface waters.The DO minimum in intermediate waters of the Indo-Pacific expanded and spread into the deep ocean under glacial conditions.The lowest value was calculated for the intermediate water box of the North Pacific where the DO concentration declined to 69 µM at 21 ka.The glacial DO decrease in the deep ocean is consistent with a large data set showing that deep waters below 1500 m water depth were significantly depleted in the LGM over all major ocean basins (Jaccard and Galbraith, 2012).The glacial oxygen depletion in the intermediate Indo-Pacific (76 µM during the LGM vs. 96 µM in the modern ocean) seems to be at odds with the geological record which shows that OMZs located in the tropical ocean were better ventilated under glacial conditions (Altabet et al., 1995).This discrepancy may probably arise from the spatial resolution of the box model that is too coarse to resolve OMZs.Moreover, all sediment cores that have been used to reconstruct the oxygen conditions in glacial intermediate waters were taken at continental margins (Jaccard and Gal- The glacial oxygen decline simulated in the model had no significant effect on other model parameters since the oxygen level stayed above the threshold values for diminished phosphorus burial (20 µM) in all ocean boxes (Wallmann, 2010).It was only benthic denitrification at the deep-sea floor that was enhanced by the glacial oxygen depletion in bottom waters and rising export production (Bohlen et al., 2012).Nevertheless, the global rate of benthic denitrification calculated in simulation STD decreased under glacial conditions (Fig. 6j) since ocean margins retreated into steeper terrain such that less nitrate was consumed on the continental shelf.Nitrogen fixation was assumed to increase when the DN / DP ratio in the surface ocean fell below the N / P ratio in exported POM (Eq.A8).Nitrogen fixation thus traced the temporal evolution of benthic denitrification (Fig. 6i-j).This negative feedback mechanism (Tyrrell, 1999;Redfield, 1958) maintained the simulated DN / DP ratio close to its modern value over the entire glacial cycle.The DN inventory peaked during the LGM where it exceeded the modern value by 16 %.A similar increase in the LGM nitrate inventory (10-30 %) and glacial decline in denitrification and nitrogen fixation was simulated with a box model constrained by the marine δ 15 N record (Deutsch et al., 2004).

Dissolved inorganic carbon, carbonate ion
concentrations, and δ 13 C of dissolved inorganic carbon In the standard case (STD) the global mean seawater concentrations of dissolved inorganic carbon (DIC) and total alkalinity (TA) decreased over the last interglacial, attained a minimum at its end (2267 and 2394 µM at 118.5 ka, respectively), increased over the glacial up to a maximum prior to the glacial termination (2467 and 2631 µM at 19.5 ka), and decreased again over the Holocene (Fig. 5b and e).On the basis of our model runs these trends were mainly driven by sea-level change that controlled the burial of neritic carbonate (Fig. 3d) and POC (Fig. 6e), and the rates of POC and carbonate weathering by shelf exposure (Figs.3h and 6c).The glacial DIC and TA rise was mitigated by nutrient utilization enhancing marine export production and carbon burial (compare simulations STD-CC-CN and STD-CC in Figs. 5 and 6).
In turn, it was amplified by the glacial decrease in deep ocean ventilation, a reduced turnover rate that also implied a decrease in marine export production, POC burial, and pelagic carbonate accumulation (STD-CC vs. STD).Considering changes in ocean volume by about 3 % (Fig. 3b), the increase in DIC over the last glacial (118.5-19.5 ka) translates into a mean DIC accumulation rate of 1.70 Tmol yr −1 .By comparison, the CO 2 uptake from the atmosphere as calculated from the glacial rate of pCO 2 decline amounts to 0.17 Tmol yr −1 .Thus, only 10 % of the glacial DIC rise was induced by CO 2 uptake from the atmosphere.According to our model, the glacial demise of neritic carbon pools (carbonate and POC) was the major forcing of the DIC rise, while the sequestration of atmospheric CO 2 only was of minor importance for the glacial change in seawater composition.Most of the excess DIC accumulating in the glacial ocean originated from exposed shelf carbonate and POC and from riverine DIC which was not buried due to the contraction of depositional areas at ocean margins.The accumulation of TA and DIC in the deep ocean (Fig. 9) was corroborated by a change in Atlantic deep water chemistry.As outlined above, this LGM ocean basin was filled with corrosive southernsource waters compromising the preservation of carbonates at the deep-sea floor and diminishing the rate of pelagic carbonate burial.
The global mean concentration of carbonate ions (CO 2− 3 ) in the deep ocean (> 2000 m) rose over glacial times in simulation STD due to the decline in neritic carbonate burial and dropped over the deglaciation, at least in part, due to the recovery of neritic carbonate deposition (Fig. 10a).The glacial CO 2− 3 rise was mitigated by the decline in deep ocean ventilation and increase in ocean productivity promoting the sequestration of CO 2 in the deep ocean.Interestingly, deep ocean pH and CO 2− 3 trends diverged during the transition into the LGM (Fig. 10a-b) -that is, pH dropped while CO 2− 3 was maintained at a constant level over this period (30-20 ka).Due to this divergence, the late glacial pH was lower (that is, more acidic) than the modern value, while the CO 2− 3 concentration exceeded the modern concentration in ocean deep waters.This apparent discrepancy may be explained by the fact that alkalinity and DIC were strongly elevated in late glacial seawater (Fig. 5b and e) thereby enhancing the concentrations of both H + and CO 2− 3 ions with respect to the preindustrial modern ocean.The deglacial CO 2− 3 minimum is related to ventilation pulses in the Southern Ocean and North Pacific employed in the model.Export production of CaCO 3 and pelagic carbonate burial were enhanced by these upwelling events and removed dissolved CO 2− 3 from the global ocean.The Holocene was marked by a continuous CO 2− 3 decline probably induced by the high sea-level stand promoting neritic carbonate burial.
The strong enrichment of dissolved CO 2− 3 in glacial surface waters was induced by the decline in atmospheric pCO 2 (Fig. 9).According to our standard simulation, these CO 2− 3 concentrations exceeded Holocene values down to water depths of 1000 m, likewise at northern high latitudes where deep-water formation transmitted the signature of glacial surface waters into the ocean's interior.The carbonate ion concentration was almost constant over the entire Indo-Pacific at > 1000 m water depth since the strong increase in DIC (Fig. 9) was balanced by a corresponding TA rise.These model results well compare with those of B / Ca ratios in benthic foraminifera which probably record CO 2− 3 changes in ambient bottom waters (Yu et al., 2008(Yu et al., , 2013)).The model is consistent with glacial to interglacial changes in deepsea CO 2− 3 reconstructed from this proxy (Table C1 in Appendix C).The only deviation occurs at > 4 km water depth in the Atlantic where the model predicts elevated LGM values while the data show a glacial CO 2− 3 depletion (Table C1), possibly the result of a strong east-west gradient in bottom water chemistry not resolved yet by the B / Ca data that accordingly may not be fully representative for the Atlantic at large.
The glacial distribution pattern of δ 13 C-DIC values calculated in the standard simulation (Fig. 11) is consistent with observations (Oliver et al., 2010) because glacial δ 13 C-DIC data were employed to define the LGM circulation pattern (Sect.2, Table A5 in Appendix A).In all simulations global mean δ 13 C-DIC values mirror inversely DIC concentrations (Fig. 5b-c) getting depleted with rising DIC and vice versa.This anti-correlation is linked to the turnover of POC being strongly depleted in 13 C as compared to average seawater.The glacial demise of the sedimentary POC pool, induced by the weathering of exposed shelf sediments and the decline in depositional areas along continental margins, contributed significantly to the glacial DIC rise and affected the isotopic evolution of seawater (Broecker, 1982b;Wallmann, 2014).The glacial δ 13 C-DIC depletion was widely ascribed to a glacial loss of terrestrial biomass (Shackleton, 1977;Köhler and Fischer, 2004).However, our model can reproduce almost the entire glacial shift to depleted δ 13 C-DIC values recorded in benthic foraminifera (0.34 ± 0.19 ‰; Peterson et al., 2014) without invoking any net changes in terrestrial biomass (Table A5).This outcome is consistent with results of a new model study suggesting that the rise in carbon buried in permafrost and under ice largely compensated for the decline in peat, soil, and biomass carbon over the LGM (Brovkin and Ganopolski, 2015).C4). 14 C-DIC values represent the difference between seawater 14 C-DIC and the atmospheric value (0 ‰ for the pre-anthropogenic modern, +446 ‰ for the LGM).See legend of Fig. 8 for further information.

Radiocarbon
Atmospheric 14 C-CO 2 was forced to follow the IntCal13 values derived from the geological record by varying the 14 Cproduction rate in the atmosphere (Appendix A, Sect.A8).The 14 C production rate was calculated for each time step considering the 14 C transfer from the atmosphere to the oceans and changes in the inventories of atmospheric CO 2 and 14 C-CO 2 . 14C-CO 2 and the model-derived difference between global mean 14 C-DIC and atmospheric 14 C-CO 2 ( 14 C-DIC = 14 C-DIC -14 C-CO 2 ) were anticorrelated, since the 14 C uptake from the atmosphere was insufficient to raise the radiocarbon content of the model ocean up to the level of elevated 14 C-CO 2 values attained during periods of strong atmospheric radiocarbon production (Figs. 5d,4g,and 12).This anti-correlation was also observed in an additional simulation (Fig. 12b) where the carbon cycle and ocean circulation were maintained at a steady state representing Holocene boundary conditions, whereas 14 C-CO 2 values were forced to follow the IntCal13 record.In further steady-state simulations the radiocarbon production rate in the atmosphere was held constant over time while the carbon cycle operated in a Holocene steady-state mode with a constant pCO 2 value of 280 µatm.These model runs show that the steady-state 14 C-DIC values attained after about 100 kyr simulation time decreased with increasing production rate (Appendix C, Table C2).As previously shown, changing production rates of radiocarbon in the atmosphere may affect the difference between 14 C values in planktonic and benthic foraminifera (Adkins and Boyle, 1997) and the contrast between atmospheric and marine 14 C values (Franke et al., 2008).However, these effects were regarded as transient features induced by a slow 14 C transfer from the atmosphere into the ocean.In contrast, the results of our steady-state model suggest that 14 C depletion of the ocean during periods of elevated atmospheric 14 C production can be a permanent steady-state feature (Table C2), a conclusion further substantiated by a simple steady-state model presented in Appendix D.
The standard simulation STD yielded a mean ocean 14 C-DIC of −270 ‰ for the LGM (21 ka) and −152 ‰ at 0 ka corresponding to a glacial 14 C-DIC decline by 118 ‰.The simulations depicted in Fig. 12b show that multiple processes contributed to the 14 C glacial depletion of the ocean with respect to the atmosphere.These include the glacial changes in deep ocean ventilation, sea-level, nutrient utilization, and atmospheric radiocarbon production.The 14 C-DIC decline observed upon sea-level fall in part was induced by the glacial decline in sedimentary carbon pools adding fossil carbon to the global ocean.Glacial changes in ocean circulation contributed to the 14 C-DIC decline since the glacial demise of ventilation across the 2000 m depth horizon isolated the deep ocean from the atmosphere.However, the simulations suggest that changes in ocean ventilation possibly were responsible for less than one-third of 14 C-DIC values calculated as difference between radiocarbon in seawater DIC and atmospheric CO 2 .The results of simulations STD-CC-CN, STD-CC, and STD are indicated as blue, red, and black lines, respectively.The green line indicates the results obtained in a steady-state simulation under Holocene boundary conditions where all variables except atmospheric 14 C-CO 2 and radiocarbon production rate were kept constant over time.(c) Production rates of radiocarbon in the atmosphere calculated in the model runs and normalized to the pre-anthropogenic modern value (1.64 atoms cm 2 s −1 ).(d) 14 C production rates calculated from the geo-magnetic record (Laj et al., 2002); 10 Be production rate as reconstructed from Greenland ice-core data (Muscheler et al., 2005;Reimer et al., 2013), 10 Be production rate as reconstructed from sediment data (Frank et al., 1997); all rates are normalized to their preanthropogenic modern values.
the glacial rise in the radiocarbon contrast between global ocean and atmosphere.
Atmospheric 14 C production rates calculated in the model showed the same trends for all simulations (Fig. 12c).They attained very high values at 25 ka and declined over time.The only significant difference between the model runs occurred during the deglaciation.The standard simulation yielded elevated production rates for this period since the rapid ventilation of the deep ocean considered in simulation STD drew radiocarbon from the atmosphere and released 14 Cdepleted CO 2 into the atmosphere such that the production rate was enhanced to maintain atmospheric 14 C-CO 2 at the deglacial level documented by IntCal13 (Fig. 12a).Our model results may support the hypothesis that a significant fraction of the 14 C-CO 2 record is controlled by changes in atmospheric radiocarbon production (Köhler et al., 2006;Broecker et al., 2004).However, in contrast to our model approach various authors have proposed that most of the 14 C-CO 2 record can be explained by changes in glacial ocean dynamics and carbon cycling without invoking significantly elevated rates of atmospheric radiocarbon production (Muscheler et al., 2005;Robinson et al., 2005).The Holocene trends calculated in the model are similar to those observed in the 10 Be ice-core record (Muscheler et al., 2005) and derived from Holocene geo-magnetic data (Laj et al., 2002) while the glacial values are closer to the stacked sedimentary 10 Be record (Frank et al., 1997).Various reasons have been evoked to explain the deviations between different records of atmospheric radionuclide production (Köhler et al., 2006).The controversy suggests a clear need to develop a better constrained record of atmospheric 14 C production suitable for model validations (Fig. 12d).
The spatial distribution of radiocarbon in the global ocean changed significantly during the LGM (Fig. 11).According to the standard simulation  et al., 2013).If this correlation also holds for the glacial ocean, glacial 14 C-DIC values may be used as a proxy for DIC concentrations in glacial seawater (Sarnthein et al., 2013).The model results show that the correlation was indeed maintained in the glacial ocean and the slope of the correlation was similar for all model runs and time slices (Fig. 13).However, the regression line for glacial conditions was shifted to lower DIC and 14 C-DIC values due to changes in ocean carbon cycling and possibly elevated radiocarbon production rates in the glacial atmosphere.Thus 14 C-DIC values may serve as a new proxy for DIC concentrations in ancient seawater, if suitable methods are found to correct for the glacial shift observed in the simulations (Fig. 13).The overall LGM pattern calculated in the standard simulation (Fig. 11) compares well with trends derived from the radiocarbon contents of planktonic and benthic foraminifera, even though radiocarbon data indicate strong gradients within ocean basins, which were not resolved by the box model (Appendix C, Table C3).A recent review of glacial 14 C-DIC data (Sarnthein et al., 2013) revealed radiocarbon depletions in the deep Atlantic, Southern Ocean, and Indo-Pacific broadly consistent with those calculated by the model.However, the model was not able to reproduce very strong radiocarbon depletions measured at some deep water sites due to its coarse spatial resolution (Table C3).Moreover, it predicts significant 14 C-depletions in the Atlantic thermocline which are inconsistent with coral 14 C data (Robinson et al., 2005). 14C measurements in foraminiferal shells and corals from the glacial ocean feature strong spatial and temporal variability (Broecker et al., 2004;Sarnthein et al., 2013).More data will help to resolve this variability and constrain the radiocarbon distribution and dynamics of the glacial ocean.

Conclusions
For a first time we show model results that are consistent with both the atmospheric pCO 2 record (Figs. 5 and 7) and data on past distribution changes of dissolved oxygen, carbonate, and radiocarbon in the glacial ocean (Figs. 8,9,11 and Tables C1 and C3).Atmospheric pCO 2 and the glacial distribution of seawater tracers were not prescribed but calculated as prognostic model variables.Only marine δ 13 C data were used to parametrize the glacial circulation model.A comprehensive formulation of shelf processes and sea-level effects is a major new component included in our Earth system model.Thus, the conformity between independent proxy data and key model results (atmospheric pCO 2 change over the last 130 kyr, distribution of dissolved oxygen, carbonate, and radiocarbon in the LGM ocean) supports our hypothesis that the glacial sea-level drop induced a decline in atmospheric pCO 2 and a rise in the inventories of nutrients, DIC, and alkalinity in the glacial ocean (Wallmann, 2014).Also, we  Monnin et al., 2001Monnin et al., , 2004;;Petit et al., 1999;Waelbroeck et al., 2002;Stanford et al., 2011) reflect the internal non-linear dynamics of the Earth system and its response to external insolation forcing.
first show that the slope of DIC vs. radiocarbon observed in the modern deep ocean (Sarnthein et al., 2013) was probably maintained in the glacial ocean (Fig. 13).However, a glacial shift in the intercept now complicates the use of 14 C as DIC proxy.
The shelf hypothesis was originally developed to explain the deglacial rise in atmospheric pCO 2 (Broecker, 1982a).In contrast, our model analysis reveals that shelf and sea-level effects were not responsible for this rapid rise but account for a major portion of the slow glacial decline of atmospheric pCO 2 (Figs. 5 and 7).The deglacial sea-level rise induced a decline in nutrient and carbon stocks in the global ocean.However, these stocks changed only slowly due to their large size (Menviel et al., 2012).In contrast to the "early anthropogenic hypothesis" (Ruddiman et al., 2016), we attribute the gradual pCO 2 rise over the Holocene to the slow relaxation of nutrient and carbon stocks promoting CO 2 transfer from the ocean into the atmosphere.The slow relaxation may also be responsible for the imbalance in phosphate and TA sources and sinks observed in the modern ocean (Wallmann, 2010(Wallmann, , 2014)).Stocks of these chemical species may decline until today since tens of thousands of years may be needed to draw down the dissolved P and TA inventories from their peak values attained over the last glacial maximum.
According to standard Milankovitch theory (Milankovitch, 1941), variations in summer insolation at high latitudes (Berger and Loutre, 1991) cause waxing and waning of northern ice sheets (Fig. 14).Most of the global climate change over a glacial cycle is thus believed to be driven by northern summer insolation and ice sheet dynamics (Denton et al., 2010).However, it has always been difficult to explain why atmospheric pCO 2 declined over glacial periods and how this drop was connected to the buildup of large continental ice sheets.The sea-level effects explored in this paper provide the missing link between glacial ice sheet and pCO 2 dynamics.The sea-level-driven pCO 2 decline was amplified by a decrease in deep ocean ventilation, a decline in sea surface temperature, and enhanced nutrient utilization.These additional changes were driven by a combination of greenhouse gas, albedo, and insolation forcing (Fig. 14).Glacial terminations occurred when summer insolation increased at northern latitudes (Raymo et al., 1997), ice sheets reached a critical size (Denton et al., 2010), and carbonate compensation at the deep-sea floor reversed the declining pCO 2 trend (Wallmann, 2014).The deglacial warming was again driven by greenhouse gas, albedo, and insolation forcing promoting the retreat of continental ice sheets, sea-level rise, ocean ventilation, and the decline in nutrient utilization in a positive feedback mode.
Due to their internal non-linear dynamics, continental ice sheets are able to generate 100 kyr cycles with a slow glacial expansion and rapid deglacial contraction of ice volume under Milankovitch forcing even though insolation oscillates on much shorter timescales (Imbrie and Imbrie, 1980;Abe-Ouchi et al., 2013;Ganopolski and Calov, 2011;Pollard, 1983).Positive feedbacks embedded in the global carbon cycle are able to generate a 100 kyr cycle without any form of external forcing when surface temperature, ice volume, sealevel, and ocean circulation are assumed to be controlled by pCO 2 (Wallmann, 2014).Thus, both, continental ice sheets and the global carbon system have the inherent tendency to generate cycles with a length of 100 kyr.They interact via sea-level and pCO 2 change, respond to insolation forcing, control changes in the climate system (surface temperature, ocean and atmospheric circulation) and may generate the 100 kyr cycle dominating late Quaternary climate change.

A1 Data and procedures for model calibration
Mean tracer concentrations were calculated for each of the model boxes using the GLODAP database for total alkalinity, DIC, 13 C-DIC, and 14 C-DIC (Key et al., 2004) and the World Ocean Atlas (WOA01) for temperature, salinity, PO 4 , NO 3 , and O 2 (Conkright et al., 2002). 14C-DIC data were corrected by subtracting the bomb-14 C signal and DIC data were corrected for the intrusion of anthropogenic CO 2 (Key et al., 2004), whereas 13 C-DIC data were not corrected and are thus affected by 13 C-depleted anthropogenic CO 2 .The model was run into steady state under pre-anthropogenic boundary conditions and resulting tracer concentrations were compared to data to validate the model output and calibrate the model (Tables A1 and A2, Fig. A1).The pCO 2 value and global export production of particulate organic carbon were used as additional constraints, i.e. the calculated values had to comply with the corresponding observations (ca.280 µatm and ca.700-900 Tmol yr −1 , respectively, Sarmiento and Gruber, 2006).For these initial model runs, the riverine fluxes to the ocean were enhanced to compensate for the removal fluxes observed in the modern ocean.The atmospheric pCO 2 value was calculated by applying a constant continental CO 2 uptake rate balancing the CO 2 being produced in the modern ocean by carbonate burial and degassing processes (Wallmann, 2014).The isotopic composition of atmospheric CO 2 was maintained at a constant level representative for the pre-human atmosphere (δ 13 C = −6.5 ‰, 14 C = 0‰).Water fluxes (Table A3) and parameter values for key biogeochemical processes (Table A4) were varied until pCO 2 , global export production, and the tracer distribution fields generated by the steady-state box model were consistent with data (Tables A1 and A2, Fig. A1).A good fit was obtained for all tracers except 13 C-DIC.Some of the 13 C mismatch was induced by anthropogenic 13 C which has a strong effect on the observations in the modern ocean but was not considered in the model simulations.

A2 Ocean circulation and tracer transport
Water fluxes between adjacent boxes were calculated by using output of the Opa 8.2 ocean circulation model in the framework of NEMO (Nucleus for European Modelling of the Ocean) to start with a configuration that is dynamically consistent with a 3-D forward ocean model (Madec et al., 1998).NEMO was forced by atmospheric reanalysis data as described in Aumont and Bopp (2006).A more detailed analysis of the resulting large-scale circulation pattern is given in Bordelon- Katrynski and Schneider (2012).The calculation of water exchange was based on horizontal and vertical velocities on the box model grid.In order to consider two-way exchange between the boxes, not only the net transports, but fluxes in both directions (northward/southward, up/down) were taken into account.Water fluxes at the atmosphereocean boundary were calculated as residual fluxes balancing the water exchange for each vertical column.Test runs with the box model revealed, however, that tracer distributions calculated with the NEMO-derived water fluxes were inconsistent with tracer data when the NEMO circulation field was applied in the box model.This mismatch was induced by the coarse spatial resolution of the box model and by errors inherent to the NEMO simulations.The NEMOderived circulation field was, thus, modified to allow for a better fit to the independent observations listed in Tables A1  and A2 and to bring the circulation field in line with other observations and GCM modelling results.Thus, NEMO features a North Atlantic Deep Water (NADW) formation rate of only ca. 10 Sv whereas tracer data (radiocarbon, phosphate, oxygen) constrain this rate at ca. 15 Sv (Broecker et al., 1998).Moreover, NEMO predicts that Antarctic Bottom Water (AABW) upwells in the Indo-Pacific all the way to the thermocline and surface ocean.This pattern is consistent with the "great ocean conveyer" (Broecker, 1991) but in conflict with other more recent ocean models suggesting that deep water ascent occurs in the Southern Ocean rather than in the Indo-Pacific (Sarmiento and Gruber, 2006).These models show that AABW flowing into the Indo-Pacific returns as deep water to the Southern Ocean where it ascends to form intermediate water masses flowing northwards into the major ocean basins (Imbrie et al., 1993;Gnanadesekian and Hallberg, 2002;Marinov et al., 2006).NEMO also predicts an extremely high rate of vertical water exchange in the Southern Ocean across 2000 m water depth of more than 200 Sv.These strong upward and downward water fluxes are inconsistent with tracer observations showing strong vertical gradients between the deep ocean (> 2000 m water depth) and the overlying water masses.The water fluxes derived from NEMO were modified to remove these biases and to provide more realistic water fluxes for the box model (Table A3).The corresponding best-fit water fluxes are bidirectional, i.e. water flows in both directions between each of the adjacent model boxes (Table A3).The net fluxes (Table A3) were calculated as the difference between these opposing fluxes.They represent the meridional overturning circulation (MOC) as implemented in the box model (Fig. 2): NADW is formed in the North Atlantic and Arctic basins at an overall rate of ca. 15 Sv.It flows towards the Southern Ocean where Antarctic Bottom Water (AABW) is formed at a rate of ca.18 Sv.A minor AABW fraction flows northwards into the Atlantic while most of the AABW is filling the deep basins of the Indo-Pacific at a rate of ca.16 Sv where it upwells and returns into the Southern Ocean as deep water.Intermediate water formed by deep water ascent in the Southern Ocean flows into the Indo-Pacific at a rate of ca. 14 Sv where it upwells to form surface water flowing back towards the Southern Ocean.Surface waters flowing northward into the Atlantic and returning as NADW to the Southern Ocean are closing the loop.This overall MOC pattern and the corre- sponding flow rates are consistent with tracer data and other observations (Sarmiento and Gruber, 2006).The exchange fluxes between adjacent boxes reflect, both, opposing water flows across the box boundaries and eddy diffusive mixing (Table A3).The large vertical exchange flux between tropical surface and intermediate waters in the Indo-Pacific is thus supported by intense Ekman-driven upwelling and downwelling while wind-driven eddy diffusive mixing explains most of the vertical exchange between surface and intermediate water boxes in the North Atlantic, Southern Ocean, and North Pacific.The overall vertical water exchange across the 100 m water depth horizon (ca.297 Sv, Table A3) is sufficiently high to ventilate the global thermocline and to support a global rate of new and export production in the order of 700-900 Tmol yr −1 .About 37 % of this vertical exchange flux occurs in the Southern Ocean (> 30 • S).The global bidirectional water flux across the 2000 m water depth level is much lower (only 45 Sv, 44 % in the Southern Ocean) while the flux across the 4000 m line amounts to ca. 179 Sv with a 74 % contribution by the Southern Ocean.The box model's major internal boundary for vertical exchange is thus located at 2000 m water depth between the thermocline and the underlying deep ocean.
In box modelling, water fluxes (F Wab ) are multiplied by tracer concentrations (C j ) to calculate tracer fluxes (F Tab ) between adjacent boxes: where F Wab and F Tab are the water and tracer fluxes from box a to box b while C a is the concentration of the considered tracer in box a.The back fluxes from box b to box a are defined correspondingly: Tracer fluxes arising from the water exchange fluxes listed in Table A3 (F Wex ) are thus proportional to the concentration difference between adjacent boxes: These fluxes can be regarded as diffusion-analogue mass transfer processes since their magnitude is proportional to concentration differences rather than concentrations.In contrast, the tracer fluxes arising from the net water fluxes in   A1 and A2).
where F TIS is the tracer flux from the intermediate water box to the overlying surface water box, F WIS is the corresponding water flux, and C TH is the tracer concentration in the upwelling thermocline water calculated as  DIC, and 13 C-DIC.A smaller value was applied for total alkalinity (f I = 0.3) to mimic the deeper regeneration of this tracer while larger values were applied for other tracers featuring steeper thermocline gradients (DP and DN: f I = 0.55; DO: f I = 0.9; 14 C-DIC: 0.7).

A3 Salinity and surface temperatures
Salinity (Sal) was treated as an inert tracer.The mass balance equations for Sal in each of the 24 ocean boxes were thus simply defined as Freshwater fluxes were varied in the initial steady-state simulations until the calculated salinity values were consistent with observations.These simulations showed, however, that unrealistically high freshwater fluxes were needed to reproduce the low salinity values observed in the Arctic Ocean surface water box.This problem arises since the standard box model procedure demands that water masses leaving the Arctic carry a chemical signature corresponding to the mean salinity value integrated over the entire Arctic surface ocean.Observations show, however, that Arctic surface waters sinking into the abyss are more salty than mean Arctic surface water.The standard model procedure was, hence, modified to consider this characteristic feature of deep water formation and to avoid unrealistically high freshwater fluxes to the Arctic surface ocean -that is an enhanced salinity (+ 0.9 PSU with respect to the mean salinity of Arctic surface water) was employed for the waters sinking into the underlying intermediate water box.A1 and A2).Over a glacial cycle, surface temperatures are regulated by changes in, both, albedo and the partial pressure of greenhouse gases.In the model it was assumed that 50 % of the temperature change is proportional to the prescribed sealevel change (e.g.continental ice sheet formation and albedo change) while the remaining 50 % is proportional to the log-arithm of atmospheric pCO 2 calculated as prognostic model variable.The global mean atmospheric surface temperature was assumed to fall by 3 • C during the glacial while the average SST was allowed to drop by ca. 2 • C (Schmittner et al., 2011).The SST drop was assumed to be twice as high as the global mean at high latitudes and only half as high in the Clim.Past, 12, 339-375, 2016 www.clim-past.net/12/339/2016/low-latitude surface water boxes.Temperatures in intermediate and deep waters were maintained at their modern values, for simplicity.

A4 Phosphorus
The model includes a comprehensive phosphorus cycle.Rivers transport dissolved phosphorus (DP) into the ocean where it is taken up by phytoplankton, gets exported, degraded, buried in marine sediments, and removed via hydrothermal activity (Wallmann, 2014).Export production (F EPOP ) of particulate organic P (POP) from the individual surface water boxes across 100 m water depth was calculated applying Liebig's law: (A7) POP export was thus limited either by dissolved reactive nitrogen (DN) or DP where K DP is a Monod constant (K DP = 0.01 µM), r NP is the atomic N to P ratio in exported particulate organic matter (r NP = 17, Körtzinger et al., 2001), Vol S is the volume of the considered surface water box, and k EXP is a site-specific kinetic constant defined by fitting the model to DP concentrations observed in the modern surface ocean (Table A4).Most of the exported POP was degraded in the water column while a small but significant fraction was permanently buried in marine sediments (Appendix B).POP degradation in the water column (including the bioturbated surface layer of marine sediments) was distributed between intermediate (93 %), deep (6 %), and bottom water boxes (1 %).Export and degradation of particulate organic carbon (POC) and nitrogen (PON) and oxygen respiration were derived from the corresponding POP turnover applying constant Redfield ratios (PON/POP = 17, POC/PON = 123/17, O 2 /POC = 1.34,Körtzinger et al., 2001).

A5 Nitrogen
Nitrogen cycling was simulated considering export, degradation, and burial of PON, nitrogen fixation, benthic and pelagic denitrification, and riverine fluxes of dissolved reactive nitrogen (DN).Nitrogen fixation in surface water boxes (F NF ) was calculated as It was controlled by the ambient DP concentration and modulated by the DN / DP ratio in surface water (r DNDPS ) such that nitrogen fixation decreased when r DNDPS exceeded the N / P ratio in exported biomass (r NP ) and vice versa.The biomass of nitrogen-fixing organisms was completely degraded in the surface ocean boxes and the organic nitrogen compounds were transformed into DN.Nitrogen fixation thus enhanced the DN pool in the surface ocean but did not contribute to export production.The kinetic constant for nitrogen fixation (k NF ) was determined by fitting the model to the nitrate concentrations observed in the modern surface ocean (Tables A1 and A4).Denitrification in the water column was limited to the intermediate water box of the Tropical Indo-Pacific where the major oxygen minimum zones (OMZs) are located.It proceeded at a constant rate of 5 Tmol yr −1 (Deutsch et al., 2001).No attempt was made to simulate the temporal evolution of water column denitrification since the box model did not resolve the dynamics and spatial extent of OMZs.Benthic denitrification was calculated as a function of bottom water chemistry (nitrate and oxygen) and the rain rate of particulate organic carbon (POC) to the seafloor (Appendix B) using an empirical transfer function (Bohlen et al., 2012).

A6 Oxygen
Oxygen was produced in the surface ocean via export production and consumed in the ocean's interior by degradation of particulate organic matter.The oxygen exchange between surface ocean and atmosphere was calculated as where k W is piston velocity, A SUR is the ice-free surface area of the considered box (Köhler et al., 2005), DO S is the concentration of DO in the considered surface water box, while DO SEQ is the temperature-dependent equilibrium concentration of DO in surface water (García and Gordon, 1992).The piston velocity was determined for each surface water box by fitting the model to the observed 14 C values (Tables A2   and A4).The deep ocean is ventilated by cold and oxygenenriched surface waters sinking into the ocean's interior at high latitudes with a temperature of ca.−1.5 • C. To mimic this process in the box model, the oxygen concentration in downward-flowing water masses was calculated by applying a temperature of −1.5 • C rather than the significantly higher mean SSTs of North Atlantic and Southern Ocean surface waters.

A7 Carbon
POC burial depended on export production and depositional area while neritic carbonate burial was proportional to the shelf area at 0-50 m water depths (Appendix B).Export of pelagic PIC (particulate inorganic carbon) was calculated from POC export production and the PIC / POC export ratios which were derived by fitting the TA values in the surface ocean to observations (Tables A2 and A4).Exported PIC dissolved in intermediate water boxes until the PIC / POC export ratio at 2000 m water depth reached a value of unity as observed in sediment trap studies (Berelson et al., 2007;Honjo et al., 2008).Carbonate compensation was implemented at > 2000 m water depth where PIC dissolution was controlled by the carbonate ion concentrations calculated from TA and DIC values in deep and bottom water boxes (Wallmann, 2014).The remaining PIC was buried in pelagic sediments.The TA mass balance equations considered al-kalinity production via denitrification and PON export production and alkalinity consumption via nitrogen fixation and PON degradation.The CO 2 gas flux from the surface ocean into the atmosphere across the seawater/atmosphere boundary layer (F CO 2 ) was calculated as (Sarmiento and Gruber, 2006) where CO 2S is the concentration of CO 2 in the considered surface water box (as calculated from ambient DIC and TA), while CO 2SEQ is the equilibrium concentration of CO 2 in surface water (as calculated from atmospheric pCO 2 ).The thermodynamic equations included in the box model considered the effects of SST and salinity on CO 2SEQ and CO 2S (Zeebe and Wolf-Gladrow, 2001).Constant rates were applied for on-shore volcanic and metamorphic degassing, degassing at mid-ocean ridges, alteration of oceanic crust, and silicate weathering (Wallmann, 2014).The rate of silicate weathering was set to a constant value since the weathering of exposed shelf sediments was assumed to compensate for the glacial decrease in silicate weathering in the continental hinterland (Munhoven, 2002).The rate of carbonate weathering was assumed to depend on surface temperature, run-off, and the size of the exposed shelf area (Wallmann, 2014).Riverine POC fluxes are ignored in the model.However, POC weathering is considered.It has two components: (i) weathering of POC in exposed shelf sediments and (ii) weathering of fossil POC in continental hinterland (Wallmann, 2014).Both components produce atmospheric CO 2 depleted in 13 C and 14 C.

A8 Carbon isotopes
The model includes 13 C-DIC and 14 C-DIC as tracers in addition to total DIC.Isotope ratios as well as δ 13 C and 14 C values of DIC were calculated from 13 C-DIC/DIC and 14 C-DIC/DIC mole fractions.The gas exchange of 13 C-CO 2 across the seawater/atmosphere boundary layer (F 13CO2 ) was calculated as (Schmittner et al., 2013;Zhang et al., 1995) where α aq−g is the equilibrium fractionation factor for CO 2 gas exchange between seawater and air, α k is the corresponding kinetic fractionation factor, α DIC−g is the equilibrium fractionation factor defining the 13 C fractionation between DIC and gaseous CO 2 , R 13DIC is the 13 C / 12 C ratio in DIC, and R 13CO2A is the 13 C / 12 C ratio in atmospheric CO 2 .The isotopic composition of DIC species (CO 2 , HCO − 3 , CO 2− 3 ) was calculated using equilibrium fractionation factors given in Zeebe and Wolf-Gladrow (2001).These values were applied to calculate the isotopic composition of exported POC and neritic and pelagic carbonates applying isotopic fractionation factors according to Ridgewell (2001) and Romanek et al. (1992), respectively.
Figure A3 shows the atmospheric δ 13 C-CO 2 record as calculated in our standard simulation.According to this simulation and our previous studies (Sarnthein et al., 2013), the strong negative δ 13 C-excursion observed in the ice-core record was largely caused by deglacial upwelling pulses in the Southern Ocean, though the amplitude of the simulated δ 13 C-CO 2 decline is much larger than that observed in the data set (Schmitt et al., 2012).The trends of our results and the empiric data are similar.The different extent of shift of the two records is probably related to the poor representation of the 13 C-DIC turnover in the Southern Ocean in our coarse-resolution model.We were not able to reproduce more closely the observed δ 13 C-DIC distribution in our model calibration even though a good fit was attained for all other tracers (salinity, DIC, 14 C-DIC, alkalinity, phosphate, nitrate, oxygen, Fig. A1).Hence, we do not conclude that this deviation for δ 13 C-CO 2 implies erroneous model results for all remaining model variables.It rather reflects a specific weakness in the model set-up with respect to the simulation of δ 13 C-DIC in the Southern Ocean.Our model predicts a negative value for intermediate waters in the modern Southern Ocean (δ 13 C-DIC = −0.05‰) while observations yield a positive value of +0.72 ‰ for this ocean box (SO I in Table A2).Due to this deviation the vertical gradient in the model exceeds the observed δ 13 C gradient between surface and intermediate waters by more than a factor of two.The glacial rise and deglacial drop in δ 13 C-CO 2 are amplified by this model artefact.The deglacial intermediate water, overly depleted in 13 C, upwells into the surface ocean where it induces a far too strong atmospheric δ 13 C-CO 2 decline.Since the biased vertical δ 13 C-DIC gradient in the Southern Ocean impedes a meaningful simulation of atmospheric δ 13 C-CO 2 , we tuned the 13 C-CO 2 fluxes between the surface ocean and the atmosphere such that the resulting atmospheric δ 13 C-CO 2 values were consistent with the ice-core record.By this way we effectively employed the ice-core data to force the δ 13 C-DIC model (Fig. 4f).The δ 13 C values of intermediate, deep and bottom water boxes employed to derive the glacial circulation field (Table A5) were not significantly affected by this tuning since the inventory of 13 C residing in the global ocean exceeds the atmospheric inventory by almost two orders of magnitude.
The 14 C / 12 C fractionation between DIC, CO 2 , POC, and CaCO 3 was calculated using the squared 13 C equilibrium fractionation factors since the mass difference between 14 C and 12 C exceeds the 13 C-12 C difference by a factor of 2. The gas exchange of 14 C-CO 2 across the seawater/atmosphere boundary layer (F 14CO2 ) was thus calculated as where R 14DIC is the 14 C / 12 C ratio in DIC of the considered surface water box and R 14CO2A the 14 C / 12 C ratio in atmospheric CO 2 .Moreover, 14 C-DIC was subject to radioactive decay with a decay constant of λ = 1 8267 yr −1 . 14C-DIC values were calculated as where 14DIC is the 14 C mole fraction ( 14DIC = 14 C-DIC/DIC), abs is the 14 C mole fraction of the standard (1.175 × 10 −12 ; Mook and Plicht, 1999) pre-human atmosphere with a δ 13 C value of −25 ‰, f N is the normalization factor defined as and δ 13 C is the δ 13 C value of DIC in the considered box in ‰ PDB.The isotopic fractionation experienced by the considered DIC pool was thus taken into account in the calculation of marine 14 C-DIC values by applying the δ 13 C-DIC calculated for the considered ocean box (Stuiver and Polach, 1977).The time-dependent radiocarbon production rate in the atmosphere (R 14 ) was calculated by applying R 14 = k 14 ( 14 C − CO 2 (data) where k 14 is a constant (≥ 10 5 mmol yr −1 ‰), 14 C-CO 2 (model) is the atmospheric value calculated for each time step of the model, and 14 C-CO 2 (data) is the data trend reconstructed from the geological record (Reimer et al., 2013).R 14 thus increased when 14 C-CO 2 (model) was smaller than 14 C-CO 2 (data) and vice versa.With this approach, the production rate was varied such that the model always complied with the atmospheric 14 C record.The atmospheric radiocarbon model considered production and the decay of radiocarbon in the atmosphere as well as exchange processes with the continents and the surface ocean (Eq.A12).The major output of the atmospheric 14 C model was the timedependent 14 C production rate (Eq.A15).

B1 Particulate organic carbon (POC) turnover
The overwhelming portion of POC produced in the euphotic zone is degraded in the water column before it can reach the seabed.Hence, in the open ocean, only ca. 1 % of the primary production is deposited at the deep-sea floor (Suess, 1980;Jahnke, 1996;Seiter et al., 2005;Dunne et al., 2007).However, the fraction reaching the seabed increases drastically at continental margins where shallow water depths limit the transit time of POC sinking through the water column.
Global models and observations thus indicate that ca.30 % of ambient primary production reaches the shallow seafloor at 0-50 m water depth (Dunne et al., 2007).Due to this effect and the high productivity of continental margins, the margin seabed located at < 2 km water depths receives ca.85 ± 15 % of the global POC rain rate (Dunne et al., 2007;Burdige, 2007) even though only 16 % of the global seabed is located at < 2 km water depth (Eakins and Sharman, 2012).Continental margins are even more dominant in terms of POC burial because burial is promoted by the deposition of riverine particles (Berner, 1982(Berner, , 2004) ) accumulating mostly on the continental shelf during interglacial sea-level high-stands (Burwicz et al., 2011).Thus, 90 ± 10 % of the global POC burial takes place at < 2 km water depth (Dunne et al., 2007;Burdige, 2007;Wallmann et al., 2012).POC rain and burial rates at continental margins declined during glacial sea-level low-stands since the oceans retreated into steeper terrains.During the LGM when the sea-level was 120 m lower than today, the shelf seafloor area at 0-100 m contracted by 73 % while the outer shelf and upper slope area located at 100-2000 m water depth was reduced by 13 % (Eakins and Sharman, 2012) neglecting isostatic adjustment.Considering the high rain rates at shallow water depths (Dunne et al., 2007), the glacial margin contraction diminished the global POC rain rate by possibly up to 50 %.The burial rate of marine POC may have been reduced by a similar proportion since POC burial is ultimately limited by the amount of POC reaching the seabed.However, there are a number of additional factors that affect the rate of POC burial.These include bulk sedimentation rate (Berner, 1982), surface area of sediment particles (Mayer et al., 2004), oxygen exposure time (Hartnett et al., 1998), and the resuspension and down-slope transport of shelf POC promoting POC burial at the upper slope (Walsh et al., 1981;Dale et al., 2015).These secondary processes control the burial efficiency of POC, that is the ratio between POC burial and POC rain rate.
POC burial efficiency is to a large degree controlled by sedimentation processes on the shelf that are strongly affected by sea-level change.At high sea-level most of the riverine particle load is deposited on the shelf (Burwicz et al., 2011) while low sea-level stands promote down-slope trans-port (Hay and Southam, 1977).Hence, data on Quaternary shelf and deep-sea fan sedimentation clearly show that the riverine particle flux was discharged over the shelf edge onto deep-sea fans and abyssal plains by turbidity currents over most of the glacial period (Hay and Southam, 1977;Hay, 1994;Schlünz et al., 1999).The corresponding increase in sedimentation rate probably led to a rise in burial efficiency and POC burial at > 2 km water depths (Burwicz et al., 2011;Wallmann, 2014).The efficiency of POC burial at the continental rise and deep-sea floor may have been further amplified by the glacial decline in dissolved oxygen concentrations in the deep ocean (Jaccard and Galbraith, 2012) favouring the preservation of POC in marine sediments (Hartnett et al., 1998;Dale et al., 2015).It is difficult to validate glacial changes in burial efficiency at < 2 km water depths.Here, POC preservation was possibly reduced by the intense ventilation of the glacial thermocline (Jaccard and Galbraith, 2012) while preservation might have been enhanced if sedimentation rates on the outer shelf and upper slope were significantly elevated by the glacial loss of inner shelf regions.Considering the available evidence it can be concluded that the glacial marine regression induced a strong decline in POC burial on the continental shelf while POC burial was enhanced at the continental rise and deep-sea floor.The overall effect was an increase in water column degradation and decline in marine POC burial since the focus of POC burial was shifted to > 2 km water depth where rain and burial rates are limited by the almost complete degradation of marine POC in the water column.

B2 Nutrient turnover
Continental margins are also major sinks for nitrate and phosphate since > 50 % of the global benthic denitrification and burial of marine phosphorus occur in sediments deposited at < 2 km water depth (Bohlen et al., 2012;Archer et al., 2002;Froelich et al., 1982;Baturin and Savenko, 1997;Wallmann, 2010;Middelburg et al., 1996).These fluxes are driven by the rain of marine POM to the seabed which is focused on shallow water environments (Dunne et al., 2007).The strong decrease in shallow seafloor area during glacial marine regressions thus induced a decline in nitrate and phosphate removal contributing to the expansion of the nutrient inventory in the glacial ocean (Broecker, 1982b;Deutsch et al., 2004;Eugster et al., 2013;Wallmann, 2014).A negative feedback was probably established whereby the expansion of the nutrient inventory induced a rise in export production and rain rate to the seabed which in turn promoted the burial of POC and removal of nutrients from the ocean.The glacial decline in POC burial and the glacial rise in the standing stocks of macronutrients may have been mitigated by this negative feedback mechanism (Middelburg et al., 1996).
Phosphate cycling in marine sediments is affected by oxygen conditions in ambient bottom waters and sediments (Krom and Berner, 1981;Van Cappellen and Ingall, 1994; Wallmann, 2003).Phosphate is released from sediments under suboxic and anoxic conditions due to the reduction of iron and manganese oxides and the preferential degradation of P-bearing organic matter (POP).However, a large fraction of the released phosphate is precipitated and retained in the sediment as authigenic carbonate fluorapatite (CFA).Hence, OMZ sediments are depleted in Fe / Mn-bound P, enriched in CFA and characterized by high POC / POP ratios exceeding the Redfield ratio by a factor of 2-8 (Schenau and De Lange, 2001;Lomnitz et al., 2015).Ratios between POC and reactive P (P reac : sum of POP, CFA, and Fe / Mn-bound P) amount to POC/P reac = 100-300 in OMZ sediments and 20-70 in continental margin sediments underlying oxygenated bottom waters (Schenau and De Lange, 2001;Noffke et al., 2012).Hence, the burial efficiency of reactive P and total P is reduced under low-oxygen conditions (Ingall and Jahnke, 1994;Schenau and De Lange, 2001) whereas POC is more efficiently buried in OMZ sediments covered by oxygendepleted bottom waters (Dale et al., 2015).

B3 Equations and parameter values employed to simulate carbon and nutrient turnover
The parametrization of margin processes applied in the model is summarized in Tables B1 and B2.A very simple approach was chosen to calculate the burial of neritic carbonates.It was assumed that the burial rate is proportional to the seafloor area at 0-50 m water depth (A NM ) which is controlled by sea-level change, only (see Fig. 3).Burial of marine POC and phosphorus (P) at continental margins was assumed to be proportional to POC export production (F EPOC ) and the depositional areas at 0-100 m (A S ) and 100-2000 m (A I ) that were controlled by sea-level change (Fig. 3).The decline in seafloor area at 0-2000 m water depths was applied to parametrize the rise in POC burial at the deep-sea floor induced by the glacial marine regression.We thus effectively assumed that the burial efficiency of POC remained constant at < 2 km water depth but increased at > 2km during glacial sea-level low-stands.
Phosphorus considered in the model refers to the sum of organic and reactive inorganic P phases (CFA and Fe / Mnbound P).The latter fractions contribute strongly to P burial since most of the marine particulate P deposited at the seafloor is degraded and transformed into inorganic authigenic phases during early diagenesis (Ruttenberg and Berner, 1993).P burial decreases when the dissolved oxygen content of ambient bottom waters (DO) falls below a threshold value of about 20 µM (Wallmann, 2010).This effect was taken into account by introducing corresponding Monod terms in the P burial flux definitions (Table B2).Moreover, the model formulation ensured that the molar POC / P reac burial ratio did not exceed the maximum value of about 400 observed in Quaternary sediments (Anderson et al., 2001).
Benthic denitrification was calculated from the POC rain rate and ambient dissolved oxygen and nitrate concentrations (DN) applying an empirical transfer function calibrated by in-situ benthic flux data (Bohlen et al., 2012).The rain rates needed for this function were derived from POC burial rates using the corresponding burial efficiencies.Marine PON burial was calculated from POC burial applying a molar PON / POC ratio of 17/123 (Körtzinger et al., 2001).
In contrast to previous sediment models (Heinze et al., 1999;Gehlen et al., 2006), our model does not resolve transport processes and reactions within surface sediments.We prefer to employ observational data on POC and P burial and empirical transfer functions to constrain benthic turnover rates (Bohlen et al., 2012;Wallmann, 2010) because most depth-resolving transport-reaction models yield results that are not yet consistent with key data such as benthic oxygen and nitrate fluxes and POC burial rates (Stolpovsky et al., 2015).
The applied parameter values were constrained by field data from the modern ocean (Table B1).The available POC and P burial data suggest that burial is rather evenly distributed between 0-100 and 100-2000 m water depths (Ta-ble B1) whereas POM rain rates at 0-100 m clearly exceed the corresponding rates at 100-2000 m (Dunne et al., 2007).This difference is caused by bottom currents transporting marine POM from the inner shelf towards outer shelf and upper slope environments (Walsh et al., 1981).Note that the POC burial rates applied in our model (Table B1) are conservative -that is, lower than most previous estimates (Burdige, 2007;Dunne et al., 2007).The burial efficiency applied in the model is low at shallow water depth since winnowing by bottom currents affects large parts of the shallow seafloor such that about 70 % of the modern shelf sediments are non-accumulating, relict sands with very low POC contents (Burdige, 2007).POM exported laterally from these shallow areas provides POM for slope deposits (Walsh et al., 1981).The highest burial efficiency is applied at 100-2000 m water depth where the deposition of fine-grained riverine particles and low oxygen values in ambient bottom waters favour the preservation of marine POC.Parameter values and fluxes listed in Tables B1 and B2 refer to global fluxes.These fluxes were distributed among the six ocean basins defined in the box model considering the respective seafloor areas and export productions.Neritic carbonate burial was distributed between the Tropical Indo-Pacific and Atlantic (Kleypas, 1997).
Our model predicts that the global burial rates of marine POC and P would decline to 7.9 and 0.14 Tmol yr −1 during the LGM, respectively, if export production and oxygen concentrations were maintained at their modern value.The export production was, however, promoted by the decline in P burial such that the best-fit simulation STD produced LGM burial rates of 9.6 and 0.165 Tmol yr −1 compared to modern global rates of 11.5 and 0.18 Tmol yr −1 , respectively (see Table B1).Our standard model run suggests that the shelf (0-100 m water depths) trapped a total of 4650 Gt POC over the last glacial cycle (130-0 ka) while 7870 Gt POC accumulated on the continental slope (100-2000 m water depth).Shelf weathering released a total of 1940 GtC over the last 130 ka in simulation STD -namely, less than 50 % of the POC accumulating on the shelf over the last glacial cycle.According to the model, the glacial marine regression affected the chemical and isotopic composition of seawater and the CO 2 content of the atmosphere via a chain of interconnected processes: ocean margins retreated into steeper terrain and shelf areas were exposed by the marine regression; the burial of phosphorus and neritic carbonate and benthic denitrification declined due to the steepening of ocean margins; carbonate, POC, and P weathering rates increased due to the exposure of shelf sediments; atmospheric CO 2 was consumed and converted into dissolved alkalinity by enhanced carbonate weathering while isotopically depleted CO 2 was released into the atmosphere by POC weathering; standing stocks of dissolved nutrients and alkalinity in the ocean expanded due to the decrease in burial and denitrification and the increase in weathering; export production rose due to the increase in the dissolved nutrient stocks and CO 2 was transferred from www.clim-past.net/12/339/2016/Clim.Past, 12, 339-375, 2016 the atmosphere into the ocean interior by the intensified biological pump while CO 2 sequestration was supported by enhanced seawater alkalinity.

B4 Model limitations
Model parametrizations were chosen to the best of our knowledge.It should, however, be noted that key processes such as glacial changes in POC and P burial efficiency are only poorly constrained by available data.Moreover, we assumed that the global mean morphology of continental margins was retained over the glacial cycle and that the average global change in relative sea-level was equal to eustatic sealevel change.This approach neglects the glacial isostatic adjustment, i.e. the glacial subsidence of northern land masses loaded by large ice sheets, the uplift in flanking regions, and the numerous far-field effects (Daly, 1934;Milne and Mitrovica, 2008).It also neglects changes in margin morphology induced by the erosion and down-slope transport of shelf sediments during glacial sea-level low-stands (Hay, 1994) and the deglacial tilting of continental margins (Clark et al., 1978).Moreover, we assumed that during the LGM eustatic sea-level was 120 m lower than today while growing evidence supports the view that LGM sea-level fall was in fact larger than this consensus value (Austermann et al., 2013;Lambeck et al., 2014).The changes in the size of depositional and exposed areas at continental margins applied in the model should thus be regarded as rough estimates.Clearly, more work needs to be done to improve these estimates.However, there is no doubt that ocean margins retreated into steeper terrain while large shelf areas were exposed during glacial marine regressions and that these changes had a profound effect on glacial seawater composition and atmospheric pCO 2 .
The accumulation of terrestrial POC in marine sediments was not considered in our model.Moreover, our model did not consider the growth of land plants and soil formation on emerged shelf regions during glacial sea-level low-stands.Trees and other plants may use sediment nutrients after shelf exposure and accumulate terrestrial POC on the emerged shelf.However, we think that the POC accumulation associated with these processes is small compared to the sedimentary POC turnover considered in the model.Modern continental margins (shelf and rise) accumulate sedimentary POC at a rate of about 100-200 Gt kyr −1 (Hedges and Keil, 1995;Burdige, 2007;Wallmann et al., 2012;Dunne et al., 2007).This enormous flux is induced by the high marine productivity of the region and the rapid accumulation of sediments facilitating the burial of marine POC.Trees and soils growing on the emerged shelf would have to accumulate POC in the order of 10 000 Gt C to maintain this high carbon flux over the glacial period (ca.80 kyr), an unlikely scenario since the global terrestrial carbon stock is ≤ 5000 GtC.The standing stock of POC in margin sediments exceeds the global terrestrial stock since POC is buried and preserved more efficiently in sediments than in most soils and plants.Sedimentary POC burial and preservation are promoted by high sedimentation rates and the lack of oxygen in these water-saturated deposits.
Table C1.Deep-sea CO 2− 3 concentrations (in µmol kg −1 ): differences between 21 and 5 ka calculated in the standard simulation STD are compared to corresponding observational data (derived from B / Ca ratios in foraminifera, LGM-Holocene differences).

Figure 1 .
Figure 1.Morphology of global ocean margin.The black line is the cumulated seafloor area as derived from the high-resolution ETOPO 1 grid(Eakins and Sharman, 2012).The ocean margin at 0-100 m water depth is indicated for the modern ocean (red area) and for the LGM when eustatic sea-level was lowered by 120 m (blue area).The global ocean retreated into steeper terrain during the glacial marine regression.The seafloor areas covered by shallow waters were reduced by this steepening of ocean margins.

Figure 2 .
Figure 2. Set-up of the box model: the global ocean is separated into 24 boxes representing surface (0-100 m), intermediate (100-2000 m), deep (2000-4000 m), and bottom (> 4000 m) waters in the Arctic (AR), North Atlantic (NA), Tropical Atlantic (TA), Southern Ocean (SO), Tropical Indo-Pacific (TIP), and North Pacific (NP).Arrows with numbers indicate net water fluxes between boxes in Sv; major fluxes (> 5 Sv) are represented by large arrows, minor fluxes (< 5 Sv) by small arrows.Arrows crossing the top boundary of the surface water boxes (seawater-atmosphere interface) indicate net freshwater fluxes (precipitation + river water fluxes -evaporation).The upper panel shows the circulation field applied for the modern ocean and the previous interglacial; the lower panel shows the circulation applied over the LGM.

Figure 3 .
Figure 3. Model forcing related to sea-level change.(a) Eustatic sea-level (Waelbroeck et al., 2002; Stanford et al., 2011); (b) global ocean volume as calculated from eustatic sea-level and ocean bathymetry data (Eakins and Sharman, 2012); (c) salinity of global mean seawater as calculated from global ocean volume;(d) global burial rate of neritic carbonate as calculated from seafloor area at 0-50 m water depth(Kleypas, 1997;Wallmann, 2014); (ef) seafloor area at 0-100 m and 100-2000 m water depth calculated from sea-level and ocean bathymetry data(Eakins and Sharman, 2012); (g) exposed shelf area calculated from sea-level and ocean bathymetry data(Eakins and Sharman, 2012); (h) global rate of POC weathering calculated from exposed shelf area(Wallmann, 2014).

Figure 4 .
Figure 4. Model forcing applied to define ocean circulation, nutrient utilization in the Southern Ocean, and the isotopic composition of atmospheric CO 2 .(a-b) Net water fluxes between the Southern Ocean (SO) and Tropical Atlantic (TA).The horizontal flows are given for (from top to bottom) surface water (subscript S), intermediate water (subscript I), deep water (subscript D), and bottom water (subscript B); (c) horizontal exchange flux between the Southern Ocean and Tropical Indo-Pacific intermediate waters; (d) vertical water exchange fluxes in the Southern Ocean and North Pacific across 100 m water depth (solid lines) and 2000 m water depth (broken line); (e) nutrient utilization in the Southern Ocean (Martınez-Garcia et al., 2014); (f) δ 13 C value of atmospheric CO 2 .Dots indicate ice-core data(Schmitt et al., 2012) while the solid line defines the values applied in the model.For > 24 kyr BP, where data are not available, the δ 13 C-CO 2 value is set to −6.4 ‰; (g) 14 C value of atmospheric CO 2 , dots indicate values reconstructed from the geological record(Reimer et al., 2013) while the solid line defines the values applied in the model.For > 50 kyr BP, where data are not available, the atmospheric 14 C-CO 2 is assumed to correspond to the pre-anthropogenic modern value (0 ‰).

Figure 7 .
Figure 7. Atmospheric pCO 2 over the last 25 kyr.(a) Model results for simulations STD (black line), STD-CC (red line), and STD-CC-CN (blue line); solid dots indicate ice-core data by Monnin et al. (2001, 2004) while ice-core data reported in Marcott et al. (2014) are shown as open circles; (b) relative contribution of sealevel change (blue), nutrient utilization (red), and ocean circulation changes (black) to pCO 2 model results; the left-hand column shows contributions to the glacial pCO 2 drawdown (51 % induced by sealevel fall, 34 % by enhanced nutrient utilization, 15 % by changes in ocean circulation); the right-hand column indicates the driving forces for the deglacial pCO 2 rise (23 % induced by sea-level rise, 43 % by decrease in nutrient utilization, 34 % by changes in ocean circulation).

Figure 8 .
Figure 8. Concentrations (in µM) of dissolved phosphorus (DP) and oxygen (DO) in the pre-anthropogenic modern ocean (PRE, model results for 0 ka, a: DP, c: DO) and during the LGM (model results for 21 ka, b: DP, d: DO, TableC4).The contour plots shown here and in Figs. 9 and 11 are based on concentrations calculated in simulation STD for each of the 24 ocean boxes (indicated as grid points).
) led to a reduction in the global water exchange across the 2000 m depth horizon from a modern amount of 45 Sv down to 31 Sv at 21 ka.This corresponds to an increase in the average residence time of water in the deep ocean (> 2000 m) from 470 years in the modern ocean to 680 years during the LGM where the residence time is calculated as the ratio of the deep ocean volume (6.65 × 10 17 m 3 at > 2000 m) and the global vertical water fluxes across 2000 m.The glacial increase in residence time

Figure 10 .
Figure 10.Global mean carbonate ion concentrations (a) and pH values (b) below 2000 m water depth for simulations STD (black line), STD-CC (red line), and STD-CC-CN (blue line).

Figure 11 .
Figure 11.Isotopic composition of dissolved inorganic carbon (DIC) in the global ocean.δ 13 C of dissolved inorganic carbon (δ 13 C-DIC in ‰) and radiocarbon composition of dissolved inorganic carbon ( 14 C-DIC in ‰) in the pre-anthropogenic modern ocean (PRE, model results for 0 ka, a: δ 13 C-DIC, c: 14 C-DIC) and during the LGM (model results for 21 ka, b: δ 13 C-DIC, d: 14 C-DIC, TableC4). 14C-DIC values represent the difference between seawater 14 C-DIC and the atmospheric value (0 ‰ for the pre-anthropogenic modern,

Figure 12 .
Figure 12.Radiocarbon values and production rates.(a) Atmospheric 14 C-CO 2 ; dots are IntCal13 data (Reimer et al., 2013) while the black line shows the values applied in the model runs.(b) Marine14 C-DIC values calculated as difference between radiocarbon in seawater DIC and atmospheric CO 2 .The results of simulations STD-CC-CN, STD-CC, and STD are indicated as blue, red, and black lines, respectively.The green line indicates the results obtained in a steady-state simulation under Holocene boundary conditions where all variables except atmospheric 14 C-CO 2 and radiocarbon production rate were kept constant over time.(c) Production rates of radiocarbon in the atmosphere calculated in the model runs and normalized to the pre-anthropogenic modern value (1.64 atoms cm 2 s −1 ).(d)14 C production rates calculated from the geo-magnetic record(Laj et al., 2002); 10 Be production rate as reconstructed from Greenland ice-core data(Muscheler et al., 2005;Reimer et al., 2013), 10 Be production rate as reconstructed from sediment data(Frank et al., 1997); all rates are normalized to their preanthropogenic modern values.

Figure 13 .
Figure 13.DIC vs. difference between radiocarbon in seawater DIC and atmospheric CO 2 ( 14 C-DIC) at > 2000 m water depth.Data are mean values for deep water and bottom water boxes derived from water column measurements.Model results are shown for the standard case (STD), for constant circulation (STD-CC), and constant values for circulation and nutrient utilization (STD-CC-CN).LGM refers to model results at 21 ka.
14 C-DIC values were strongly depleted over the entire ocean and reached a minimum of −356 ‰ in North Pacific deep water.Both, vertical and horizontal gradients were strengthened during the LGM.In the modern ocean, marine 14 C-DIC values are correlated with www.clim-past.net/12/339/2016/Clim.Past, 12, 339-375, 2016 DIC concentrations at water depths below 2000 m (Sarnthein

Figure 14 .
Figure 14.Key elements of the 100 kyr cycle.Summer insolation at high northern latitudes (June insolation at 60 • N, diagram in the upper right corner, Berger and Loutre, 1991) affects the growth and melting of continental ice sheets and thereby eustatic sea-level change.The glacial drawdown of atmospheric pCO 2 and its deglacial rise are supported by sea-level change.The cycle is closed by atmospheric pCO 2 affecting global climate and thereby the volume of continental ice sheets.It is accelerated and further strengthened by additional positive feedbacks: ice sheets affect the Earth's albedo and climate while changes in ocean and atmosphere circulation support the glacial pCO 2 drawdown and are largely responsible for the rapid deglacial rise in atmospheric pCO 2 .The records of atmospheric pCO 2 and eustatic sea-level change (diagrams in the upper left corner;Monnin et al., 2001Monnin et al., , 2004;;Petit et al., 1999;Waelbroeck et al., 2002;Stanford et al., 2011) reflect the internal non-linear dynamics of the Earth system and its response to external insolation forcing.

Figure A1 .
Figure A1.Tracer concentrations in ocean boxes: model vs. data.Open circles indicate concentrations applied as initial values at 130 ka (TablesA1 and A2).Crosses are concentrations obtained at the end of simulation STD at 0 ka after the completion of a full glacial cycle.Lines indicate the 1 : 1 relationship, e.g. the best fit to the data.δ 13 C data at < 2000 m water depth are excluded from the model-data comparison since they are affected by anthropogenic CO 2 .
Figure A1.Tracer concentrations in ocean boxes: model vs. data.Open circles indicate concentrations applied as initial values at 130 ka (TablesA1 and A2).Crosses are concentrations obtained at the end of simulation STD at 0 ka after the completion of a full glacial cycle.Lines indicate the 1 : 1 relationship, e.g. the best fit to the data.δ 13 C data at < 2000 m water depth are excluded from the model-data comparison since they are affected by anthropogenic CO 2 .

Figure A2 .
Figure A2.Deglacial benthic -pelagic radiocarbon record in the North Pacific: model (line) vs. data (squares).The 14 C-DIC difference between deep water and surface water boxes in the NorthPacific is compared to data ( 14 C difference between benthic and pelagic foraminifera, B-P 14 C) from core MD02-2489 taken at the Alaskan Margin at 3.6 km water depth(Rae et al., 2014).The vertical mixing across 100 and 2000 m water depth was enhanced by a factor of 10 at 16.5 ka (Fig.4d) to reproduce the radiocarbon data.
A5)    where C I is the concentration in the intermediate water box, C S the concentration in the overlying surface water box, and f I defines the fraction of C I in the ascending two-component mixture.The weighing factor f I was set to 0.5 for salinity, A6)where F Sal-in are the salinity fluxes from the neighbouring boxes into the considered box while F Sal-out gives the corresponding fluxes from the considered box into the adjacent boxes.The volume of the considered box (Vol) was allowed to change over time to mimic the contraction of the ocean volume during glacial sea-level low-stands.Fluxes were calculated as products of water flux and salinity in the source box (Eqs.A1 and A2) with the exception of intermediate to surface water fluxes where thermocline concentrations were applied (Eqs.A4 and A5).Surface water boxes were sub-ject to evaporation and precipitation and received river input from the continents.The freshwater fluxes listed in Table A3 represent the overall budget of evapotranspiration and river water input.Negative fluxes thus indicate that evaporation exceeds the sum of precipitation and river water input.

Figure A3 .
Figure A3.Stable carbon isotopic composition of atmospheric CO 2 .Model results vs. ice-core data(Schmitt et al., 2012).The model results were obtained in the standard simulation STD without any tuning of the 13 C-CO 2 fluxes across the seawater-atmosphere interface.
i = (S, I, D, B) F BDENi = F BPOCi BE i • (0.06 + 0.19 • 0.99 (DO i −DN i ) ) * Subscripts indicate modern values (M) and the following environments: shelf (S, 0-100 m water depth), outer shelf and slope (I, 100-2000 m), continental rise and deep-sea floor (D, 2000-4000 m), deep-sea floor and abyssal plain (B, > 4000 m).The equations define global fluxes.These were distributed among the ocean basins according to their export production and the seafloor areas of individual boxes.

Table 1 .
Controls on atmospheric pCO 2 and mean dissolved carbon and phosphorus concentrations in the global ocean.
* Corrected for anthropogenic CO 2 ; + not corrected for anthropogenic CO 2 .

Table A3 .
Water fluxes (in Sv) derived from the NEMO model run and fluxes applied in the box model to reproduce observed tracer distributions (Tables

Table A4 .
Biogeochemical parameter values determined by fitting the model to observations.

Table B1 .
Parameter values applied in the simulation of margin processes.

Table B2 .
Flux parametrizations applied in the simulation of margin processes.

Table C4 .
LGM tracer concentrations in model boxes at 21 ka in simulation STD (see TablesA1 and A2for further information).