Evaluating the Biological Pump Efficiency of the Last Glacial Maximum Ocean using δ 13 C

. Although both physical and biological marine changes are required to explain the 100 ppm lower atmospheric p CO 2 of the Last Glacial Maximum (LGM, ~21 ka) as compared to pre-industrial (PI) times, their exact contributions are debated. Proxies of past marine carbon cycling (such as δ 13 C) document these changes, and thus provide constraints for quantifying the drivers of long-term carbon cycle variability. This modelling study 15 explores the relative rolespresents a realization of the of physical and biological changes in the ocean needed to simulate an LGM ocean in satisfactory agreement with proxy data, and here especially δ 13 C. We prepared a PI and LGM equilibrium simulation using the ocean model state (NorESM-OC) with full biogeochemistry (including the carbon isotopes δ 13 C and radiocarbon) and dynamic sea ice. The modelled LGM-PI differences are evaluated against a wide range of physical and biogeochemical proxy data, and show agreement for key aspects of the 20 physical ocean state within the data uncertainties. However, the lack of a simulated increase of regenerated nutrients for the LGM indicates that additional biogeochemical changes are required to simulate an LGM ocean in agreement with proxy data. In order to examine these changes, we explore the potential (offline) theoretical effects of different global mean biological pump efficiencies on the simulated marine biogeochemical tracer distributions. We estimate that (besides changes in ocean an approximate doubling of the global mean pump efficiency from 38 % (PI) to 75 % (LGM) model-proxy biases the remaining absolute model-proxy error in δ 13 C is 0.07 ‰ larger than the 0.19 ‰ data uncertainty) indicates that additional changes in ocean dynamics are needed to simulate an LGM ocean in agreement with proxy data, on which we include a detailed discussion.such as increased aging or volume of Southern Source Waters. Besides that, our theoretical approach of increasing the biological pump efficiency may be too simplified to capture the vertical redistribution of regenerated nutrients – also suggested by a too weak chemocline. Our results underline that only those coupled climate models that contain the processes and/or components that realistically change both ocean circulation and biogeochemistry will be able to simulate an LGM ocean in satisfactory agreement with proxy data – and hence be reliable for use in climate projections. Therefore, future research should aim to identify the exact physical and biogeochemical processes that could have doubled the global mean biological pump efficiency (i.e., the interior regenerated signature) between the PI and LGM, with a likely central role for Southern Source Waters.


Introduction
Model and proxy reconstructions of the Last Glacial Maximum (LGM) suggest major redistributions of marine biogeochemical tracers and water masses as compared to pre-industrial (PI) times, as well as lower carbon storage in both the land biosphere and atmosphere. The culmination of these changes into a ~100 ppm lower LGM atmospheric pCO2 concentration (pCO2 atm ; EPICA Project Members, 2004) has driven extensive research to 5 identify, understand, and quantify the processes contributing to these major atmospheric pCO2 atm variations (e.g., Broecker, 1982;Broecker and Peng, 1986;Heinze and Hasselmann, 1993;Heinze et al., 2016;Sigman et al., 2010;Adkins, 2013;Jeltsch-Thömmes et al., 2019;Rae et al., 2018). The oceans are of particular interest as they form the largest carbon reservoir available for atmospheric exchange on millennial timescales, and in addition need to have stored the extra carbon coming from the land biosphere and atmosphere during the LGM. Both physical 10 (circulation, solubility) and biological processes (biological pump efficiency) likely played a role in the differences between the LGM and PI oceans, although their relative importance is under debate.: Between ~25 and ~60 % is attributed to biological processes and the remainder to physical changes (Bouttes et al., 2011;Buchanan et al., 2016;Khatiwala et al., 2019;Yamamoto et al., 2019;Kobayashi et al., 2015;Stein et al., 2020;Marzocchi and Jansen, 2019). 15 Here, we explore thethe relative roles of marine physical and biological changes needed to simulate an LGM ocean in optimal agreement with proxy data. We use the concept of the biological pump efficiency (defined as the ability of marine organisms to consume surface ocean phosphate, or more specifically the ratio of global mean regenerated to total phosphate, Sect. 2.4) to examine its effect on LGM marine biogeochemical tracer distributions in concert with physical changes. The global mean efficiency of the biological pump is strongly and nearly linearly correlated 20 with atmospheric CO2 concentrations (Ito and Follows, 2005) and is considered a key concept to understand the atmospheric CO2 drawdown potential of the ocean (Ödalen et al., 2018) through its influence on the vertical gradient of marine dissolved inorganic carbon (DIC). In our evaluation and discussion, we pay particular attention to the role of Southern Source Waters (SSW, waters originating in the Southern Ocean), which are thought to be a key component in altering ocean interior tracer distributions and glacial atmospheric pCO2 atm drawdown (e.g., 25 Lynch-Stieglitz et al., 2016;Schmitt et al., 2012;Moy et al., 2019;Sigman et al., 2010;Ferrari et al., 2014;Morée et al., 2018;Khatiwala et al., 2019;Rae et al., 2018;Kobayashi et al., 2015;Marzocchi and Jansen, 2019;Stein et al., 2020).
water mass ventilation and circulation rates, and for comparison with reconstructed ∆ 14 C (Skinner et al., 2017;Gebbie and Huybers, 2012). We furthermore evaluate the LGM simulation against proxy and/or model reconstructions of water mass distributions, sea surface temperature, salinity, sea ice extent, export production, vertical nutrient redistribution, atmospheric pCO2 atm , the change in marine dissolved inorganic carbon, and O2 (Sect. 3.2). We apply the concept of True Oxygen Utilization (TOU; Ito et al., 2004) instead of Apparent Oxygen 5 Utilization (AOU) and make use of the explicit simulation of preformed biogeochemical tracers in our model (Tjiputra et al., 2020). This approach makes it possible to separate physical and biogeochemical drivers of the tracer distributions more thoroughly, and accounts for the role of the air-sea carbon disequilibrium pump (Khatiwala et al., 2019). We acknowledge that without a land source of C in our simulated LGM ocean (of ~ 850 Gt C, Jeltsch-Thömmes et al., 2019), nor sediments that could alter CaCO3 cycling and long-term organic matter 10 burial (Sigman et al., 2010), we do not expect to simulate the full range of processes contributing to glacialinterglacial pCO2 atm changes. Rather, we include estimates of these carbon reservoir changes in our evaluation of the LGM biological pump efficiency (Sect. 3.3).
The evaluation against proxy data allows us to evaluate both the physical and biological changes needed for shown that changes in model (biogeochemical) parameterizations are needed to simulate glacial-interglacial cycles in agreement with proxy records (e.g., Jeltsch-Thömmes et al., 2019;Ganopolski and Brovkin, 2017;Buchanan et 20 al., 2016;Heinze et al., 1991;Heinze and Hasselmann, 1993;Heinze et al., 2016). In our forced ocean model setup, we are able to reveal aspects important for modelling the LGM and relevant for improving the agreement between fully coupled paleo modelling and proxy data. Moreover, our work will help to gain insight in the changes (i.e. physical and biological) needed to simulate a different climate state (such as the LGM) -which also applies to Earth System Model-based climate projections. 25

Model description
We apply the ocean carbon-cycle stand-alone configuration of the Norwegian Earth System Model (NorESM) as described by Schwinger et al. (2016), but with several modifications for the next generation NorESM version 2 already included. The physical ocean component MICOM (Miami Isopycnic Coordinate Ocean Model; Bentsen 30 et al., 2013) has been updated as described in Guo et al. (2019). The biogeochemistry component HAMOCC (HAMburg Ocean Carbon Cycle model;Maier-Reimer, 1993;Maier-Reimer et al., 2005) adopted for use with the isopycnic MICOM (Assmann et al., 2010;Tjiputra et al., 2013;Schwinger et al., 2016) has undergone a few minor technical improvements (e.g. updated initialisation based on latest data products, additional diagnostic tracers) and employs a new tuning of the ecosystem parameterization as described in Tjiputra et al. (2020). 35 In addition to these changes, the carbon isotopes ( 13 C and 14 C) are implemented in HAMOCC (Tjiputra et al. (2020)), a prognostic box atmosphere is made available for atmospheric CO2 (including 13 CO2 and 14 CO2; Tjiputra et al., 2020), and an LGM setup is made (Sect. 2.2). This is an ocean-only modelling study, where the atmospheric forcing is prescribed from a data set (except atmospheric CO2, δ 13 C and ∆ 14 C, which evolve freely; Sect. 2.2). All simulations in this study are done without the sediment module of HAMOCC (this is done in order to avoid prohibitively long spin-up times, especially for the carbon isotopes; an acceleration method for the model spin-up including interactive water column-sediment interaction is work in progress for a separate manuscript). Applying the current set-up, detritus arriving at the sediment-water interface is evenly redistributed over the entire water 5 column, while opal and CaCO3 are dissolved immediately in the bottom-most mass containing layer. Riverine input of carbon and nutrients is also turned off. Furthermore, nitrogen deposition, denitrification and nitrogen fixation are excluded from our simulations, as these processes cause a long-term drift in the alkalinity inventory of the ocean (and thereby the pCO2 of the prognostic atmosphere).
The two main isotopes of carbon, 13 C and 14 C, are newly implemented in HAMOCC (Tjiputra et al., 2020). The 10 model includes fractionation during air-sea gas exchange and photosynthesis, as well as radiocarbon decay.
Fractionation during CaCO3 formation is small as compared to the effects of air-sea gas exchange and photosynthesis, as well as uncertain (Zeebe and Wolf-Gladrow, 2001) and is therefore omitted (e.g., Schmittner et al., 2013;Lynch-Stieglitz et al., 1995). Air-sea gas exchange fractionation (~8-11 ‰) is a function of surface ocean temperature and the CO3 2fraction of total DIC such that fractionation increases with decreasing temperatures 15 (Zhang et al., 1995;Mook, 1986). Biological fractionation (~19 ‰) increases surface water δ 13 C of DIC while producing low-δ 13 C organic matter. In the interior ocean, this light isotope signal from organic matter is released back into the water column during remineralization and respiration, thereby creating a vertical gradient. HAMOCC applies the parameterization by Laws et al. (1997), where the biological fractionation εbio depends on the ratio between phytoplankton growth rate and the aqueous CO2 concentration. For 14 C, each fractionation factor is set to 20 the quadratic of the respective 13 C value (i.e., α14C = α 2 13C). In addition, 14 C is radioactive and decays with a halflife of 5730 years to 14 N.
In order to evaluate the carbon isotopes against observations, we derive δ 13 C and ∆ 14 C. δ 13 C is calculated using the standard equation δ 13 C = (  (Tjiputra et al., 2020). This approach facilitates comparison with the radiocarbon disequilibrium data by Skinner et al. (2017). 30

Last Glacial Maximum set-up
Several adjustments were made to the model in order to obtain an LGM circulation field. First, the land-sea mask and ocean bathymetry were adjusted for the ~120 m lower sea level in the LGM caused by the increased land ice volume as compared to the PI. Following the PMIP4 guidelines in Kageyama et al. (2017) the Bering Strait is closed, and the Canadian Archipelago (including Borrow Strait and Nares Strait), Barents Sea, Hudson Bay, Black 35 Sea, Red Sea, as well as the Baltic and North Seas are defined as land in the LGM. The PI land-sea mask formed the basis for the LGM land-sea mask, through shifting the PI bathymetry 116 m upwards. If the resulting depth in a grid cell was between 0-25 meters, the depth was set to 25 m and negative depths were set to land grid points.
After this, any channels with a width of only one grid cell were closed off as well, as these inhibit sea ice movement in the sea ice model causing unrealistic sea ice build-up. LGM freshwater runoff is routed to the nearest ocean grid cell but otherwise unadjusted.
Changes in isopycnal densities and sea surface salinity restoring are applied in the LGM model setup in order to ensure an adequate vertical model resolution and ocean circulation. A net LGM increase in density due to 5 decreased ocean temperatures and increased ocean salinity required increasing all 53 isopycnal layer densities by 1.3 kg m -3 in the LGM setup as compared to the PI model setup. NorESM-OC uses salinity restoring to avoid longterm drift away from a predefined SSS state. Here, this predefined state is chosen, consistent with the atmospheric forcing (see below) as the mean of the LGM minus PI SSS anomaly modelled by PMIP3 models added to a PI SSS climatology. However, the unadjusted application of the PMIP3-based SSS anomaly caused an Atlantic water 10 mass distribution and overturning strength in poor agreement with proxy reconstructions (SM 12). Earlier studies have shown a high sensitivity of models to SSS restoring, especially in the North Atlantic (Rahmstorf, 1996;Spence et al., 2008;Bopp et al., 2017). Indeed, the density contrast between Northern and Southern source waters drives the simulated Atlantic Meridional Overturning Circulation (AMOC) strength in many of the PMIP2 models (Weber et al., 2007), and is therefore important for the simulation of overturning strength in agreement with proxy 15 records. Therefore, we adjust the SSS restoring present in NorESM-OC to obtain a circulation field in better agreement with proxy reconstructions: In addition to the PMIP3-based SSS anomaly, we apply a region of -0.5 psu in the North Atlantic and +0.5 psu in the Southern Ocean (for specifics, see SM 1), as done similarly by Winguth et al. (1999) or through freshwater fluxes by Menviel et al. (2017) and Bopp et al. (2017).
An atmospheric LGM forcing for NorESM-OC was created by adding anomalies (relative to the pre-industrial 20 state) derived from PMIP3 models (Morée and Schwinger, 2019; version 1, SM 2) to the CORE Normal Year Forcing (NYF; Large and Yeager, 2004). The use of mean PMIP/CMIP anomalies to force stand-alone models is a standard approach that has been tested before (Mitchell et al., 2017;Chowdhury and Behera, 2019;Muglia et al., 2015;Muglia et al., 2018). Through this approach, the effect of the presence of sea ice on the atmospheric state is included in the forcing, but the sea ice model handles the actual formation/melt of sea ice. Compared to the PI 25 CORE-NYF, the LGM forcing over the ocean has a lower specific humidity (especially in the tropics), decreased downwelling longwave radiation, precipitation and air temperature, and a heterogeneous change in downwelling shortwave radiation and zonal and meridional winds. In addition to the adjustments to the NYF, the dust fluxes of Lambert et al. (2015) are used in the LGM model setup, following PMIP4 guidelines . The change to the dust forcing alters the amount of iron available for biology during photosynthesis, and is considered 30 an important component of glacial ocean biogeochemistry and pCO2 atm drawdown (e.g., Yamamoto et al., 2019).
The PI setup uses the Mahowald et al. (2006) dust dataset.

Initialization and tuning
All marine biogeochemical tracers were initialized in the LGM as done for the PI spin-up using the WOA and GLODAPv2 data sets. For the LGM, consistent with the decreased ocean volume, all biogeochemical tracer 35 concentration are increased at initialization by 3.26 %. Similarly, ocean salinity is uniformly increased by 1 psu, following PMIP recommendations. The carbon isotopes (Sect. 2.1) are only enabled after an initial spin-up of the model in order to first obtain reasonably stabilized total carbon tracer distributions. DI 13 C is initialized after 1000 years using the correlation between δ 13 C and apparent oxygen utilization (AOU) in combination with the model's −0.0075 • AOU + 1.72) and converted to absolute model 13 C using model DIC and AOU. As this approach uses the model's 'native' AOU and DIC, the equilibration time of δ 13 C was reduced as compared to initialisation with a δ 13 C data product such as that of Eide et al. (2017). Model DI 14 C is initialized after 4000 years by first calculating δ 14 C using a combination of pre-industrial δ 13 C (Eide et al., 2017) (with the missing upper 200m copied from 5 200m depth to all empty surface layers) and the observational-based estimate of pre-industrial ∆ 14 C (Key et al., 2004). Then, model DI 14 C is derived from the δ 14 C by rewriting and solving the standardization equation (δ 14 C = ( 14 ⁄ ( 14 ⁄ ) − 1) * 1000 ‰, with model DIC as C). Subsequently, isotopic DOC, POC, phytoplankton C, and zooplankton C are initialized as done for the corresponding total carbon variable, but multiplied with 0.98 (as an estimate of the photosynthetic fractionation effect) and the respective DI 13 C/DI 12 C or DI 14 C/DI 12 C ratio. Isotopic 10 CaCO3 is initialized as for total carbon, multiplied with DI 13 C/DI 12 C or DI 14 C/DI 12 C, as we do not consider fractionation during CaCO3 formation.
The prognostic atmospheric pCO2 atm is initialized at 278 ppm for both spin-ups. At initialization of the carbon isotopes, atmospheric δ 13 C is set to -6.5 ‰ and atmospheric ∆ 14 C is set to 0 ‰. Atmospheric pCO2 atm at the time of initialization is then used to calculate the absolute 13 C and 14 C model concentrations ( 13 C atm and 14 C atm , in ppm). 15 Two main spin-ups have been made with NorESM-OC: One for the LGM and one for the PI, designed as described in Sect. 2.1-2.3. Both the PI and LGM simulations are run for a total of 5600 years.

Analysis of the biological pump efficiency
Here, we explore the effect of an increase in the global mean biological pump efficiency ( ̅̅̅̅̅̅̅̅ , Eq. 1), which we define, following Ito and Follows (2005), as the ratio between global mean regenerated phosphate ( 4 ̅̅̅̅̅̅̅̅ ) and 20 global mean total phosphate ( 4 ̅̅̅̅̅ ).
Regenerated phosphate is calculated as the difference between total phosphate and preformed phosphate (PO4 pref ).
PO4 pref is explicitly simulated in the model (Tjiputra et al., 2020), and represents phosphate that leaves the mixed layer in inorganic form (unutilized by biology). 25 We work with the global mean value of ̅̅̅̅̅̅̅̅ as this governs pCO2 atm (Ito and Follows, 2005;Ödalen et al., 2018).
However, we note that major local differences in the ratio of regenerated to total phosphate exist in the ocean, for example between North Atlantic Deep Water (high-ratio) and Antarctic Bottom Water (low-ratio) (Ito and Follows, 2005;DeVries et al., 2012), which thus indicate the differences in potential to sequester carbon and nutrients in the ocean interior. Here, changes in ̅̅̅̅̅̅̅̅ are calculated in a theoretical framework (off-line) to better 30 understand the LGM redistribution of carbon between the land, atmosphere and ocean, and its effects on marine biogeochemistry (and corresponding proxy data). Our approach also allows us to give an upper estimate of the ̅̅̅̅̅̅̅̅ of the LGM ocean.
The spatial distribution of ∆ 4 is an important consideration. We therefore explore the effect of changes in ̅̅̅̅̅̅̅̅ on local ∆ 4 by applying three different methods (visualized in Fig. S3): The first method (method 'add') equally distributes the mean change in ∆ 4 over the entire ocean ( 4 = 4 + ∆ 4 ̅̅̅̅̅̅̅̅̅̅ ). 10 . The second method (method 'factor') takes into account the original distribution of 4 (by first calculating the global value ∆ 4 ̅̅̅̅̅̅̅̅̅̅ = 4 ̅̅̅̅̅̅̅̅ / 4 ̅̅̅̅̅̅̅̅ and calculating for every grid-cell The third method is as the first , additive, method but only adding the extra regenerated tracers to SSW, the location of which isas determined from the conservative PO tracer (method 'SSW', see Sect.

and Fig. 1b for the LGM PO tracer distribution). 15
It is important to note that ̅̅̅̅̅̅̅̅ can be changed by several processes: through the soft-and hard tissue biological pumps, the solubility pump (Heinze et al., 1991;Volk and Hoffert, 1985) and by changes in the physical carbon pump (circulation/stratification of the water column).

The Bern3D model
In order to put the theoretical (offline) approach applied to NorESM-OC (i.e., Sect. 2.4) in perspective, we 20 performed an additional analysis using results from the Bern3D v2.0s intermediate complexity model. The Bern3D model includes a 3-D geostrophic-frictional balance ocean (Edwards et al., 1998, Müller et al., 2006 coupled to a thermodynamic sea-ice component, a single-layer energy-moisture balance atmosphere (Ritz et al., 2011), and a 10-layer ocean sediment module (Tschumi et al., 2011). The Bern3D model components share the same horizontal resolution of 41×40 grid cells and the ocean has 32 logarithmically spaced depth layers (Roth et al., 2014). A four-25 box representation of the land biosphere (Siegenthaler and Oeschger, 1987) is used to calculate the dilution of carbon isotopic perturbations. Further information on the Bern3D model is provided in the supplementary materials (SM 3). We used the Bern3D model to estimate the LGM-PI change in the marine DIC inventory (∆DIC) in regard to changes in four observational targets (pCO2 atm , δ 13 C atm , marine δ 13 C of DIC, and deep equatorial Pacific CO3 2-), following the approach of Jeltsch-Thömmes et al. (2019). This is done by evaluating forcing-response 30 relationships in seven generic carbon cycle mechanisms obtained from factorial deglacial simulations with the Bern3D model in a reduced form emulator that is constrained with the observational targets. The seven generic forcings comprise changes in wind stress, air-sea gas exchange, the export rain ratio between organic and inorganic carbon, the remineralization depth of organic particles, coral growth, weathering of organic material, and carbon stocks in the land biosphere. More details on the experimental setup are given in Appendix A. ∆DIC is estimated 35 with this approach to be ~3900 Gt C (±1σ range of [3350,4480] Gt C), based on the four given constraints. If only a single constraint of the four is applied, ∆DIC estimates are shifted by about 1000 Gt C towards lower values (see Fig. 7). Contributions to ∆DIC from different carbon reservoirs that changed over the deglacial and eventually drew carbon from the ocean (atmosphere, land biosphere, coral reef regrowth, sedimentation-weathering imbalances) are discussed in Appendix A.

Results and discussion
The results presented in Sect. 3.1 and Sect. 3.2 are the annual mean climatologies over the last 30 years of the 5 5600-year PI and LGM spin-ups. The total integration time is long enough to eliminate most model drift at this point in time. Over the last 1000 years of the LGM spin-up simulation, ocean temperature does not show a significant trend, rather small oscillations with an amplitude of 0.015℃ (Fig. S4). SSS is stable for all practical purposes with a small positive trend of about 0.004 psu per 1000 years (Fig. S4). Likewise, AMOC strength is stable, while Drake Passage transport shows a small negative trend superimposed by variations in its strength of 10 about 2 Sv (Fig. S5). We present an evaluation of the PI (Sect. 3.1) and LGM (Sect. 3.2) spin-ups and compare the latter to proxy reconstructions, and discuss the LGM-PI changes in a theoretical framework (offline) exploring the efficiency of the biological pump (Sect. 3.3).

The simulated pre-industrial ocean
The simulated pre-industrial ocean state has a maximum AMOC strength of ~18 ± 0.5 Sv north of 20° N, which 15 compares favourably to the mean observational estimates of 17.2-18.7 Sv (Srokosz and Bryden, 2015;McCarthy et al., 2015), especially when noting the wide range of modelled AMOC strengths in similar forced ocean setups (Danabasoglu et al., 2014). The interannual variability of the simulated AMOC is small compared to observations (about ±4 Sv; Srokosz and Bryden, 2015), due to the annually repeating forcing. Drake Passage transport is simulated at ~114 Sv, lower than recent observational estimate of 173.3 ± 10.7 Sv (Donohue et al., 2016). The 20 depth of the transition between the upper and lower Atlantic overturning cells at 30° S lies at ~2700 m, comparable to other model estimates (Weber et al., 2007;Fig. S6). Temperature biases (Figs. S7a,b and S8a,b) are generally modest (smaller than ± 1.5℃) for most of the ocean above 3000m, except for a warm bias related to a too deep tropical and subtropical thermocline and a warm bias related to a too strong Mediterranean outflow at mid-depth in the Atlantic. At depth (>3000 m) there is a widespread cold bias that originates from the Southern Ocean (too 25 much deep mixing and associated heat loss to the atmosphere). Salinity biases (Figs. S7c,d and S8c,d) are generally small, except for a positive bias related to athe too strong Mediterranean outflow at mid-depth in the Atlantic. where PP is generally too low. Because of too high PP and export in the equatorial Pacific, a far too large oxygen minimum zone (OMZ) with elevated concentrations of regenerated phosphate develops in the model (Fig. S10).
Otherwise, the global nutrient concentrations are in reasonable agreement with modern observations (WOA, Glodapv2). ̅̅̅̅̅̅̅̅ is 38 % for the simulated PI ocean, in good agreement with observational estimates of 32-46 % (Ito and Follows, 2005;Primeau et al., 2013).

The physical ocean state
Proxy-based reconstructions describe an LGM circulation that includes a shoaling of the upper circulation cell in 5 the Atlantic (Glacial NADW) and expansion and slow-down of a cooler and more saline lower circulation cell (Glacial AABW; Adkins, 2013;Sigman et al., 2010;Ferrari et al., 2014). In this study, we assume these aspects of the LGM ocean to be qualitatively correct, and therefore aim for a model simulation in agreement with these features. We note that discussion continues as to the magnitude and veracity of these changes (e.g., Gebbie, 2014).
Most reconstructions estimate a weakened AMOC for the LGM as compared to the PI state, although estimates 10 vary between a 50 % weakening and an invigoration of AMOC (McManus et al., 2004;Kurahashi-Nakamura et al., 2017;Böhm et al., 2014;Lynch-Stieglitz et al., 2007;Muglia et al., 2018). The maximum overturning strength north of 20° N simulated by NorESM-OC is 15.6 Sv (~7 % weaker than simulated for our PI ocean, which we attribute to our adjustments and tuning of the salinity restoring). Higher uncertainties are involved with reconstructions of the strength of the Antarctic Circumpolar Current (ACC), with consensus leaning towards a 15 slight invigoration (Lynch-Stieglitz et al., 2016;Lamy et al., 2015;McCave et al., 2013;Mazaud et al., 2010;Buchanan et al., 2016). We simulate a Drake Passage transport of 129 Sv (varying between ~128 and ~130 Sv, Fig. S5), which is about 13 % stronger than simulated for the PI ocean. The simulated transition between the Atlantic overturning cells shoals by ~350 m (Fig. S6), which falls within the uncertainty of reconstructions (Gebbie, 2014;Adkins, 2013;Oppo et al., 2018), (Fig. S5)and is not captured in most PMIP simulations (Muglia 20 and Schmittner, 2015). Specifically, tThe transition between the overturning cells lies well above the main bathymetric features of the Atlantic Ocean in our LGM simulation (as visible in Fig. 1 from the transition line in the conservative PO water mass tracer; Fig. 1). This could have been an important feature of the glacial Atlantic water mass configuration due to reduced mixing along topography (Adkins, 2013;Ferrari et al., 2014) -i.e., shifting water mass boundaries away from the regions of intense internal mixing increases chemical and tracer 25 stratification. The changes in water mass circulation cause an increased SSW volume contribution to the Atlantic and Pacific basins, as visible from the conservative PO tracer (Fig. 1a,b, and Fig. S114 for the Pacific) (Broecker, 1974), and in agreement with proxy reconstructions. Furthermore, the Atlantic abyssal cell is slower in the LGM simulation as compared to the PI, which is important for the effects of the Southern Ocean on LGM pCO2 atm drawdown (e.g., Kobayashi et al., 2015). Indeed, Rradiocarbon age increases at depth ( Fig. 1c and Fig. S114 for 30 the Pacific), with a global volume-weighted mean increase of 269 years. This is low compared to the estimate of Skinner et al. (2017) of 689±53 years, and we find the majority of our radiocarbon age bias to lie at depth in the Atlantic (not shown) indicating too strong ventilation and/or biased equilibration of these waters (which have a southern source , Fig 1a-b). Our simulation shows a radiocarbon age increase particularly at mid-depth in the South Atlantic and South Pacific (at 2-3 km; Fig. 1c and Fig. S11). Burke et al. (2015) describe this 'mid-depth 35 radiocarbon bulge' to relate to the reorganization of the Atlantic overturning cell structure driven by Southern Ocean sea ice expansion. We simulate an increase in Southern Ocean sea ice cover for both summer and winter ( Fig. S912), but less than is inferred from proxy-based reconstructions for the LGM (Gersonde et al., 2005)similar to PMIP models (Roche et al., 2012). If the driving mechanism of the mid-depth radiocarbon bulge is indeed expanded sea ice as argued by Burke et al. (2015), improving our simulated LGM sea ice extent (Fig. S12) may increase (mid-depth) radiocarbon ages and improve agreement with proxy records like those of Skinner et al. (2017). FurthermoreBesides radiocarbon aging of in particular SSW, the SSW salinity increases (Fig. S613)alsowhich is in good agreement with proxy reconstructions as well (Adkins, 2013). The simulated increased SSW salinity, radiocarbon age and larger volume are important for simulation of LGM biogeochemistry and not a 5 straightforward combination in glacial simulations (e.g., Kobayashi et al., 2015, PMIP). Nevertheless, our simulation of the LGM physical ocean state is most likely still not fully alike the LGM ocean (see also our discussion in Sect. 4). Also relevant for the water mass circulation as well as marine biogeochemistry (Marzocchi and Jansen, 2019;Jansen, 2017), are the low LGM atmospheric temperatures that cause a mean ocean temperature decrease of 1.9 °C in the model. This is less than the 2.57 ± 0.24 °C estimated from proxy reconstructions of mean 10 ocean temperature (Bereiter et al., 2018), likely because the SSW may not carry a strong enough temperature decrease from the atmosphere into the interior ocean in our simulation (Fig. S147) -implying an underestimation of negative buoyancy fluxes. While the differences between our simulated LGM-PI changes in SST (Fig. S158) do not exhibit the same amount of heterogeneity as observed in proxy reconstructions (MARGO Project Members, 2009), the simulated mean SST change (-1.97 °C) seems reasonable when taking into account the uncertainty of 15 (2009) estimate -1.9 ± 1.8 °C; Ho and Laepple, 2015). Further, the general pattern of stronger SST cooling outside of the polar regions is captured which is important for the airsea disequilibrium pump (Khatiwala et al., 2019). We simulate an increase in Southern Ocean sea ice cover for both summer and winter (Fig. S9), but less than is inferred from proxy-based reconstructions for the LGM (Gersonde et al., 2005) -similar to PMIP models (Roche et al., 2012). Southern Ocean sea ice may have played a 20 major role in LGM marine biogeochemistry and interior ocean carbon storage (Ferrari et al., 2014;Marzocchi and Jansen, 2017;Stephens and Keeling, 2000), and could therefore explain some of the model biases. Examples are the lack of stratification in the Southern Ocean which is thought to be driven by brine rejection from sea ice (Jansen, 2017), and affect air-sea equilibration of biogeochemical tracers such as DIC and O2 (Gottschalk et al., 2016a). 25

The biogeochemical ocean state
Proxies for the past biogeochemical state of the ocean (such as export production, oxygen concentrations, δ 13 C) allow us to make a further evaluation of our simulated LGM ocean (Fig. 2). The global mean efficiency of the biological pump ̅̅̅̅̅̅̅̅ decreases from 38 % in the PI simulation to 33 % in the LGM simulation, as opposed to 30 reconstructions which infer an increased regenerated signature in the interior ocean (Jaccard et al., 2009;Umling et al., 2018;Freeman et al., 2016;Yamamoto et al., 2019). The simulated increase in preformed phosphate ( Fig. 2) represents an increased (but unused) potential for the ocean to draw down atmospheric pCO2 atm (Ödalen et al., 2018). We can attribute our simulated decrease in ̅̅̅̅̅̅̅̅ to the increase in SSW volume (Fig. 1), as SSW has a low regenerated signature (Ito and Follows, 2005). Despite the decrease in ̅̅̅̅̅̅̅̅ , simulated pCO2 atm is 21 ppm 35 lower in our LGM setup as compared to the PI setup. We attribute this change to the net effect of the i) smaller ocean volume, causing the concentration of alkalinity, DIC and salinity and a reportedly ~16 ppm pCO2 atm increase (Sarmiento and Gruber, 2006), and ii) the decrease in water temperature, which drives a pCO2 atm drawdown of ~30 ppm (Sigman and Boyle, 2000). As we made no additional changes to the marine biogeochemical model (except for thea LGM dust input field), and have no sediment or land model included in our simulation, the ~20 ppm pCO2 atm drawdown as well as limited changes in regenerated nutrient inventories is expected and found in earlier studies (e.g., Buchanan et al., 2016). Similarly, the atmospheric δ 13 C change due to glacial land-vegetation loss is not simulated because we only simulate the ocean. Physical processes thus mostly drive the simulated changes in pCO2 atm and marine δ 13 C between our LGM-PI simulations, also evidenced by the LGM increase in preformed DIC (Fig. S16). 5 Simulated changes in Atlantic total phosphate (Fig. 2) agree well qualitatively with reconstructed nutrient redistributions, which describe a deep ocean nutrient increase and mid-to surface decrease (Buchanan et al., 2016;Gebbie, 2014;Marchitto and Broecker, 2006;Oppo et al., 2018). North Pacific waters >2.5 km depth exhibit a lower LGM phosphate (and DIC) as compared to the PI, due to the lack of accumulated regenerated phosphate (Fig. S107). In agreement with the expectation of increased interior carbon storage, simulated interior DIC 10 increases -especially in SSW (Fig. 2). As for phosphate, this increase is driven by the physical carbon pump only, through higher saturation of surface DIC in the Southern Ocean driven by lower T and increased alkalinity (see also not shownSect. 4 for a discussion on the contribution from the different DIC components).
However, the biases in simulated O2 and δ 13 C LGM-PI changes and their respective proxy reconstructions are large and in disagreement with proxy data (Fig. 2). Any mismatch in the absolute values of δ 13 C is not shown here 15 because we compare LGM-PI differences in both the sediment cores and model. In line with decreased remineralisation and increased O2 solubility due to lower temperatures in the model, O2 concentrations increase throughout the ocean (Fig. 2). There is a notable difference between Northern and Southern end-members in the Atlantic: Northern-source deep waters have increased O2 concentrations due to physical O2 pumping (colder waters have higher O2 solubility) as visible in preformed O2, while SSW O2 increases due to a lack of remineralisation at 20 depth due to low oxygen utilization (Fig. 3). The general increase in O2 concentrations, mostly due to the lack of biological O2 consumption at depth (Fig. 3), is in disagreement with proxy reconstructions (Jaccard and Galbraith, 2011) and shows the missing regenerated nutrients should mostly come from lacking biological processes (remineralization). δ 13 C of DIC is governed by both ocean circulation (ventilation rate) and the efficiency of the biological pump (respiration rate), and their relative importance depends on location (Gruber et al., 1999;25 Schmittner et al., 2013;Eide et al., 2017).
As for O2, the overall increase in simulated δ 13 C of DIC contradicts δ 13 C records from sediment cores, in which the strengthening of the vertical gradient is a main feature (Fig. 2). Deep δ 13 C and its vertical gradient is for an important part governed by biological processes (Morée et al., 2018). As a last comparison, we evaluate our modelled changes in export production against proxy data, even though such data have poor global coverage and 30 large spatial heterogeneity, and are largely qualitative (Kohfeld et al., 2005). LGM export production generally decreases outside of upwelling zones in our model (Fig. 4) and increases in upwelling zones (model and proxy data, Fig. 4). We conclude that the simulated export production increase may be too weak in the sub-Antarctic and is lacking in the tropical Atlantic (Fig. 4). Especially Southern Ocean export production has a large potential to affect interior and lower latitude nutrient and DIC concentrations (Sarmiento et al., 2004;Primeau et al., 2013), 35 and likely contributes to the simulated low ̅̅̅̅̅̅̅̅ and low-regenerated signature in the biogeochemical tracer distributions. Considering the large influence of SSW nutrient supply on lower latitude productivity we anticipate a key role for SSW here -as supported by proxy data (e.g., Winckler et al., 2016;Costa et al., 2016), although local changes in iron fertilization may play an additional role depending on local iron limitation (e.g., Oka et al., 2011).

The potential of the biological pump
A decrease in pCO2 atm and increase in regenerated nutrients, despite an increase in low regenerated nutrient SSW volume, is likely to occur through the increase of (and thus regenerated nutrients) of SSW (Jaccard et al., 5 2009). In addition, an increase in the (Southern Ocean) air-sea disequilibrium of DIC may have kept more carbon sequestered in the deep ocean, through increased stratification and inhibition of air-sea gas exchange (Khatiwala et al., 2019) by for example sea ice (Marzocchi and Jansen, 20197). An increase in the regenerated signature of northern source water would have further contributed to a global increase in regenerated carbon and nutrients in the interior ocean (Yu et al., 2019), although occupying a smaller volume. As natural preformed concentrations in 10 SSW are high (Ito and Follows, 2005), there is a high potential for these waters to obtain a stronger regenerated signature, and thereby facilitate pCO2 atm drawdown (Ödalen et al., 2018). Our model results for the LGM-PI change in O2 and δ 13 C show a strong miss-match with proxy records (Sect. 3.2.2 and Fig. 2). Here, we explore the effect of a potential theoretical ('offline') increase of regenerated nutrients in the ocean, by increasing regenerated phosphate and updating O2, DIC and δ 13 C accordingly (Sect. 2.4). The increase is projected on the same simulated 15 physical ocean state presumed to represent LGM conditions (i.e. Sect 3.2.1). To the extent that this state represents true glacial conditions (see also the discussion in Sect. 4.1), it allows for an assessment of the magnitude and nature of marine biogeochemical changes needed for lowering LGM pCO2 atm . In our approach, the additional regenerated phosphate is taken from preformed phosphate, thus leaving the total phosphate inventory unchanged. Proxy reconstructions of global LGM phosphate, however, show that LGM phosphate was redistributed as well as 20 elevated (Tamburini and Föllmi, 2009;Filippelli et al., 2007). As we consider a closed system (no sediments or land input of phosphate or other elements), only redistributions of phosphate can be captured in our model setup.
We compare the mean error between the model and the δ 13 C proxy data across a wide range of ̅̅̅̅̅̅̅̅ (20-100 %, Using NorESM-OC, tThe best fit between the modelled and sediment core LGM-PI changes in δ 13 C is found for a ̅̅̅̅̅̅̅̅ of 75 % (Fig. 5). A ̅̅̅̅̅̅̅̅ of 75 % would lead to the adjusted tracer distributions shown in Fig. 6 (applying 35 Eq. 2-4). This is true for the approach 'factor' (described in Sect. 2.4), indicating that taking the distribution of the original simulated LGM PO4 reg and strengthening that regenerated signal gives the best agreement with sediment core data. The ̅̅̅̅̅̅̅̅ of 75 % corresponds to an LGM-PI ∆DIC of ~1850 Gt C. This ∆DIC estimate falls below the overall range (2400 to 5500 Gt C) given by the Bern3D model and its constraints (Sect. 2.5). If using only one of the Bern3D model above constraints (pCO2 atm , δ 13 C atm , marine δ 13 C of DIC, and deep equatorial Pacific CO3 2-), however (pCO2 atm , δ 13 C atm , marine δ 13 C of DIC, and deep equatorial Pacific CO3 2-), the ∆DIC estimate of ~1850 Gt C lies in the range of Bern3D results (see SM 3Appendix A for details). The estimated ∆DIC of ~1850 Gt C allows for full (i.e. ~80 ppm more than simulated, which is ~170 Gt C) LGM pCO2 atm drawdown, a profound decrease in land carbon (~850 Gt C as estimated in Jeltsch-Thömmes et al., 2019) as well as a source of DIC from 5 the deep ocean sediments/CaCO3. We can however not determine these sources and sinks in our model setup. % and up to 37 (i.e. 75-38) % higher than in the PI. We note that the total marine ∆DIC (~1850 Gt C) estimate points towards a central role for the Pacific basin to store glacial carbon, if we consider Atlantic LGM storage >3 km depth to be in the order of 50 Gt C . 15 In addition to evaluating the model-data fit of δ 13 C, we evaluate the effect of the tracer adjustments on O2. The decrease in O2 for an adjusted ̅̅̅̅̅̅̅̅ (Fig. 6) shows better agreement with qualitative proxy data in the Atlantic (Compare Figs. 2 and 6;Jaccard and Galbraith, 2011) as well as the estimated 175±20 μmol kg −1 LGM-PI decrease in Sub-Antarctic Atlantic bottom-water (Gottschalk et al., 2016a). Absolute LGM O2 (Fig. S181) would decrease 20 to sub-zero O2 concentrations in the North Pacific (ca. -100 µmol kg -1 ). This, which is of a similar magnitude as the size of the PI model-observation bias in this area (Tjiputra et al., 2020), but may be too extreme as indicated by qualitative proxy data that describe an LGM-PI O2 increase for the North Pacific above 3 km depth ( Fig. S170; Jaccard and Galbraith, 2011). An increase in denitrification could have played a role here, but cannot be evaluated in our model setup. Additional (quantitative) estimates of LGM-PI O2 for major water masses would thus help to 25 provide further constraints on LGM ̅̅̅̅̅̅̅̅ .
Our mean absolute error in δ 13 C decreases from 0.67 ‰ for the original LGM state estimate (Fig. 2) to 0.26 ‰ for the 75% ̅̅̅̅̅̅̅̅ ocean (Fig. 5 and 6). The remaining absolute δ 13 C error is therefore 0.07 ‰ larger than the proxy data uncertainty. After adjustment of ̅̅̅̅̅̅̅̅ to 75 %, the remaining mean model-data δ 13 C mismatch of 0.07 ‰ (Fig. 5) as well as the possibly too low Pacific O2 (Fig. S18) indicate that projecting changes in ̅̅̅̅̅̅̅̅ onto the 30 simulated glacial circulation field still does not fully align with the actual LGM state -despite exploring different approaches for the redistribution of the regenerated nutrients (Sect. 2.4). In Sect. 4.2 we discuss possible mechanisms that could contribute to understand this remaining mismatch. This remaining model-data δ 13 C mismatch of 0.07 ‰ (Fig. 5) as well as the possibly too low Pacific O2 (Fig. S10) indicate that projecting changes in ̅̅̅̅̅̅̅̅ onto the simulated glacial circulation field still does not fully align with the actual LGM state -despite 35 exploring different approaches for the redistribution of the regenerated nutrients (Sect. 2.4). Specifically, even though a 75 % ̅̅̅̅̅̅̅̅ and the reorganised circulation captures most of the magnitude of the LGM-PI δ 13 C change, the strength of the glacial chemocline (the vertical δ 13 C gradient) remains too weak (Fig. 6). Other modelling studies of the glacial ocean show similar biases (e.g., Heinze et al., 2016;Schmittner and Somes, 2016). This suggests additional processes, which would allow stronger chemical stratification between intermediate and deep waters, are missing in the model(s) and are not explained in a simple way by intensification of the biological pump or our simulated changes in circulation. We recognize that optimizing the model ̅̅̅̅̅̅̅̅ to additional (quantitative) proxies such as the nitrogen isotopes could provide more constraints (Schmittner and Somes, 2016).

Model-data bias of the LGM simulation 5
Earth System Models are generally found to incompletely capture the biogeochemistry and strengthening of the biological pump for the LGM ocean, and identification of the exact processes that are missing in these models is a major challenge in modelling the LGM ocean (e.g., Galbraith and Skinner, 2020). Previous studies show that both physical and biological marine changes must have contributed to the LGM pCO2 atm drawdown (Bouttes et al., 2011;Buchanan et al., 2016;Khatiwala et al., 2019;Yamamoto et al., 2019;Kobayashi et al., 2015;Stein et al., 10 2020;Marzocchi and Jansen, 2019). The physical ocean state of our LGM simulation compares quite well to proxy data (Sect. 3.2.1), while the simulated biogeochemical state (Sect. 3.2.2) reveals major disagreements with proxy data and a too weak LGM ̅̅̅̅̅̅̅̅ . Nevertheless, our modelled ̅̅̅̅̅̅̅̅ can be increased through both physical and biological mechanisms.
To start, our results may indicate that the water mass production processes in the model are not (yet) fully adequate 15 and contribute to the too low simulated LGM ̅̅̅̅̅̅̅̅ and too weak chemocline. Interior δ 13 C is influenced by the source of deep waters as well as intermediate waters, the extent of deep convection as well as mixing processes between interior water masses (e.g., Duplessy et al., 1988), partly on spatial and temporal scales possibly not resolved by our model. An example is the role of shelf processes and the related localized sea ice formation and brine excretion in the Southern Ocean, which are highly relevant for SSW characteristics as well as the transition 20 between the overturning cells (Klockmann et al., 2016;Mariotti et al., 2016;Bouttes et al., 2016;Jansen, 2017;Marzocchi and Jansen, 2017;Stein et al., 2020). The PI simulation of a generally too deep Southern Ocean mixed layer depth as well as Southern Ocean-attributed model biases in the PI biogeochemical tracers (Sect. 3.1; Tjiputra et al., 2020) suggest that deep water formation processes indeed are not simulated in full agreement with observational data -which may lead to biases in the LGM simulation as well. Furthermore, the lack of a reliable 25 glacial freshwater forcing is likely to be partly responsible for biases in the LGM simulation. The weak LGM-PI reduction of SSW temperature (Fig. S14) is in disagreement with proxy records and implies an underestimation of negative buoyancy fluxes. Improving this bias may however lead to more vigorous sinking of SSW in the Southern Ocean and possible speed-up of the abyssal overturning cell, equally in disagreement with proxy data.

Additionally, too strong exchange (mixing) in the LGM ocean between the deep and intermediate waters could 30
contribute to maintain a too weak glacial chemocline. Our simulated SSW radiocarbon ages are too young compared to reconstructions (Fig. 1), consistent with inadequate isolation of these waters in the simulation. Further aging of these SSW could increase their regenerated signature (consuming O2, and decreasing δ 13 C of DIC) and steepen the chemocline while additionally improving agreement with the Bern3D model ∆DIC estimate and radiocarbon age. Indeed, both increased stratification and enhanced biological pumping in the glacial Southern 35 Ocean likely played a role in deoxygenating SSW (Yamamoto et al., 2019).

and iii) Southern
Ocean sea ice cover. The simulation of a larger increase in LGM Southern Ocean sea ice extent (Fig. S12) and/or formation rate may lead to a more stratified and isolated SSW (Ferrari et al., 2014;Nadeau et al., 2019;Jansen, 2017;Jansen, 2017 and2019;Stephens and Keeling, 2000;Stein et al., 2020). This could create older waters with higher regenerated signatures and increase the air-sea disequilibrium of biogeochemical tracers 5 such as DIC and O2 (Gottschalk et al., 2016a;Khatiwala et al., 2019;Yamamoto et al., 2019), thereby further lowering atmospheric CO2 (Khatiwala et al., 2019) as well increasing SSW radiocarbon aging.
Regarding biological mechanisms that could underlie our model bias, it is instructive to first look at our simulation of a simultaneous decrease in ̅̅̅̅̅̅̅̅ and pCO2 atm (Sect. 3.2 and 4.1). This apparent contradictive result can be attributed to the fact that our model mostly captures the physical carbon pump, as driven by temperature and ocean 10 volume decrease. Indeed, decomposition of DIC into its components (DIC soft , DIC pref , DIC sat , DIC bio , DIC carb , and DIC diss : Fig. S16; Bernardello et al., 2014) shows that the increased ocean inventory of DIC (which is driving down pCO2 atm ) can be attributed to the increase in the physical DIC pref component. Specifically, DIC pref increases due to increased DIC disequilibrium (DIC diss ) of SSW with the atmosphere as well as a strengthened solubility pump as evidenced by increased DIC sat . The more isolated abyssal cell (and accompanying increase in DIC diss ) is central 15 for the lower glacial pCO2 atm as described throughout this manuscript. Nevertheless, our study shows that the increased DIC diss is not sufficient for the simulation of a full glacial pCO2 atm drawdown, and additional changes in the biological carbon pump must have played a role. In our simulation, DIC soft lowers the ocean inventory of DIC in line with the simulated decrease in ̅̅̅̅̅̅̅̅ but in disagreement with proxy records.
The lack of an increased marine C inventory through the increase in DIC soft is also evident from simulated O2 and 20 δ 13 C. In our simulation, O2 concentrations increase due to the lack of biological O2 consumption at depth in SSW and increased solubility in Northern Source Waters (Fig. 3). Such an O2 increase is in disagreement with proxy reconstructions (Jaccard and Galbraith, 2011) and argues for biological processes to explain the missing regenerated nutrients (i.e., remineralization) and low DIC soft in our simulation. As for O2, the widspread increase in simulated δ 13 C of DIC contradicts δ 13 C records from sediment cores, in which the strengthening of the vertical 25 gradient is a main feature (Fig. 2). Deep δ 13 C and its vertical gradient is for an important part governed by biological processes (Morée et al., 2018), further arguing for biological processes to contribute to our simulated biases. Nevertheless, δ 13 C of DIC is governed by both ocean circulation (ventilation rate) and the efficiency of the biological pump (respiration rate), and their relative importance depends on location in complex and highly spatially heterogeneous ways (Gruber et al., 1999;Schmittner et al., 2013;Eide et al., 2017;Morée et al., 2018). 30 A possible candidate to increase deep DIC soft is found in Southern Ocean export production, which has a large potential to affect interior and lower latitude nutrient and DIC concentrations (Sarmiento et al., 2004;Primeau et al., 2013). Iron fertilization likely played a central role for the increase in production and export fluxes as well (e.g., Winckler et al., 2016;Costa et al., 2016;Oka et al., 2011;Lambert et al., 2015;Yamamoto et al., 2019), and may be insufficiently captured by our model. Comparison between our simulation and proxy data point to an 35 underestimation of LGM export production in the ice-free Southern Ocean (Fig. 4), which we conclude may have contributed to the simulated low ̅̅̅̅̅̅̅̅ and low-regenerated signature in the biogeochemical tracer distributions.
The biological pump efficiency and hence the Δδ 13 C further depend on the details of the particle flux mode through the water column. This flux depends on the size of the sinking particles, their composition, aggregation, and weight. The plankton community structure plays an important role here (Bach et al., 2019). The vertical particle flux can be accelerated, e.g., through enhanced dust deposition (clay), biogenic opal, and calcium carbonate (Haake and Ittekkot, 1990;Armstrong et al., 2002;Klaas and Archer, 2002;Honda and Watanabe, 2010;Mendes and Thomsen, 2012). The production of transparent exopolymer particles (TEP) also modifies the vertical flux mode of organic matter substantially (De la Rocha and Passow, 2007;Mari etal., 2017). These details on marine particle fluxes still need to be accounted for in marine biogeochemical models and can have differential effects on the 5 regional biological pump efficiency.
Besides alleviating model biases in physical or biological processes, our simulation of the LGM ocean may be improved through model development: For example, the results of our simulations on pCO2 atm may be amplified through carbonate compensation in interplay with Southern Ocean stratification (Kobayashi and Oka, 2018), if we would include the ocean sediments in our simulations. The latter would also facilitate direct comparison with 10 sediment core records, especially if the carbon isotopes are included in the sediment model. Last, a decreased δ 13 C of DIC is possible without drastic increases in regenerated nutrients when considering changes in C:P ratios (Ödalen et al., 2020). Such C:P changes could contribute to resolve the model-proxy error and provide an interesting direction for future work.

The remaining model-data bias after ̅̅̅̅̅̅̅̅ adjustment
Changes between preindustrial and LGM ocean circulation fields as simulated by ocean models generally fail to 20 account for the full 100-120 ppm drawdown in atmospheric CO2 (taken the outgassing by the land biosphere into account) when used in global ocean carbon cycle models (Heinze et al., 1991;Brovkin et al., 2007). Indeed, our study gives an estimate of ~21 ppm drawdown with an unadjusted biological pump efficiency.
This remaining model-data δ 13 C mismatch of 0.07 ‰ (Fig. 5) as well as the possibly too low Pacific O2 (Fig. S10) indicate that projecting changes in ̅̅̅̅̅̅̅̅ onto the simulated glacial circulation field still does not fully align with 25 the actual LGM state -despite exploring different approaches for the redistribution of the regenerated nutrients (Sect. 2.4). Specifically, eOur offline adjustment of the ̅̅̅̅̅̅̅̅ ven though a to 75 % and % ̅̅̅̅̅̅̅̅ and the reorganised circulation captures most of the magnitude of the LGM-PI δ 13 C change (Fig. 5, 6 and S17). However,, the strength of the glacial chemocline (the vertical δ 13 C gradient) remains too weak (Fig. 6). Other modelling studies of the glacial ocean show similar biases (e.g., Heinze et al., 2016;Schmittner and Somes, 2016). This 30 suggests that additional processes, which would allow stronger chemical stratification between intermediate and deep waters, are missing in the model(s) and are not explained in a simple way by intensification of the (softtissue) biological pump or our simulated changes in circulation. We recognize that optimizing the model ̅̅̅̅̅̅̅̅ to additional (quantitative) proxies such as the nitrogen isotopes could provide more constraints (Schmittner and Somes, 2016). Besides that, our theoretical (offline) approach of increasing the biological pump efficiency may be 35 too simplified to capture the vertical redistribution of regenerated nutrients -also suggested by the too weak chemocline.
Using their idealized model setup, Schmittner and Somes (2016) estimated an LGM ̅̅̅̅̅̅̅̅ of 54-59 % (41 % in PI) and a ∆DIC of 590-790 Gt C by exploring the effects of a uniform change in the maximum growth rate of phytoplankton. A direct comparison between these studies is complicated by the large differences between the models, but the differences indicate remaining uncertainties in the magnitude of both ∆DIC and LGM ̅̅̅̅̅̅̅̅ .
Nevertheless, LGM ̅̅̅̅̅̅̅̅ seems to have been about 18 (i.e. 59-41) % and up to 37 (i.e. 75-38) % higher than in the PI. Last, we note that the total marine ∆DIC (~1850 Gt C) estimate points towards a central role for the Pacific basin to store glacial carbon, if we consider Atlantic LGM storage >3 km depth to be in the order of 50 Gt C (Yu 5 et al., 2016).

Summary and 4 Cconclusions
We present a model simulation of the pre-industrial and Last Glacial Maximum (LGM) oceans. We use the simulations to explore a realization ofthe relative roles of physical and biological marine changes needed to 10 simulate an LGM ocean in satisfactory agreement with proxy data. Despite the good agreement between (qualitative and quantitative) proxy reconstructions and our simulation of different LGM and pre-industrial (PI) ocean circulation, our model is unable to simulate the complete set of biogeochemical changes implied by proxy data. Therefore, our results (mainly the lack of increased regenerated nutrients) confirm the idea that both biogeochemical (beyond those represented by the model) and physical changes must have been at play in the ocean 15 to create the LGM pCO2 atm drawdown (Heinze et al., 2016;Bouttes et al., 2011;Buchanan et al., 2016).
Comparison between a range of qualitative and quantitative proxy data and simulated biogeochemical tracers (specifically total dissolved inorganic carbon, regenerated phosphate, True Oxygen Utilization, O2 and δ 13 C) reveals that there is a too small signature of regenerated nutrients in our simulated LGM ocean. We conduct a theoretical exploration (offline) of the effects of changes in the global mean biological pump efficiency and 20 quantify its effect on the global mean absolute error between simulated δ 13 C and proxy δ 13 C data. The smallest error is found for an approximate doubling in the global mean efficiency of the biological pump ̅̅̅̅̅̅̅̅ (from 38 % in the PI ocean to ~75 % in the LGM). Such a change minimizes the simulated global mean absolute error for δ 13 C from 0.67 ‰ (for the originally simulated 33 % ̅̅̅̅̅̅̅̅ of the LGM) to 0.26 ‰ -only distinguishable by 0.07 ‰ from the δ 13 C data uncertainty. It additionally improves the agreement with both qualitative and quantitative 25 O2 reconstructions. Much of the remaining model-proxy δ 13 C data mismatch is due to a too weak vertical chemocline in glacial simulations. Therefore, scaling of the biological pump efficiency does not fully explain the glacial ocean proxy data using the modelled circulation field thought to most closely represent the glacial state.
We see different explanations of the bias that could strengthen the glacial chemocline, such as the further aging (by isolation and stratification) of interior waters through improved simulation of deep-water formation, which 30 would increase the regenerated signature of the interior LGM ocean. Especially Southern Source Waters would have a large potential (due to their large volume contribution) to increase global interior radiocarbon ages and regenerated signatures of the interior ocean.
The estimated ̅̅̅̅̅̅̅̅ increase to 75 % corresponds to an increase in oceanic DIC (∆DIC) of ~1850 Gt C. This lies in the range of estimates as derived by the Bern3D Earth System Model of Intermediate Complexity constrained 35 by single proxy targets (pCO2 atm , δ 13 C atm , marine δ 13 C of DIC, or deep equatorial Pacific CO3 2-). If all four targets are used as constraints, the range of ∆DIC estimates based on the Bern3D model is higher than 1850 Gt C (see Fig. 7 and S3 and SM3Appendix A). In order to disentangle and understand the processes contributing to ∆DIC, especially the large contribution from sedimentation-weathering imbalances, further work seems necessary.
Our results underline that only those coupled climate models that contain the processes and/or components that realistically change both ocean circulation and biogeochemistry will be able to simulate an LGM ocean in satisfactory agreement with proxy data. Such a simulation is also a test for Earth system models for their ability to 5 reproduce natural climate variations adequately as a basis for reliable future projections, including human-induced forcing. I.e., a satisfactory fidelity of Earth System Models in reproducing orbitally forced climate variations will increase our confidence in these models as tools for projecting future anthropogenic climate change. Therefore, future research should aim to identify the exact physical and biogeochemical processes that could have doubled the global mean biological pump efficiency (i.e., the interior regenerated signature) between the PI and LGM, with 10 a likely central role for Southern Source Waters.
The general agreement between our model results and proxy data for ocean circulation (within the uncertainty of reconstructions), after adjustments to the sea surface salinity field, demonstrates an advantage of our forced ocean model setup, and its flexibility, over fully coupled Earth System Models in exploring different circulation and biological pump scenarios for explaining glacial CO2 changes. We conclude that a large PI-to-LGM increase in 15 the efficiency of the biological pump (from 38 to ~75 %) as well as a reorganization of ocean circulation/stratification are essential for simulating an LGM ocean in optimal agreement with proxy data. Based on this, we expect that only fully coupled models that contain the processes and/or components that realistically change these aspects will be able to simulate an LGM ocean in satisfactory agreement with proxy data.

Experimental set-up of the Bern3D model
We estimate the LGM-PI change in marine DIC content (∆DIC) by applying a reduced form emulator based on forcing-response relationships obtained from factorial simulation over the deglacial with the Bern3D model (see Jeltsch-Thömmes et al., 2019). The mean over 21 ka BP to 19 ka BP and over 500 to 200 yr BP are used as LGM and PI time intervals for the calculation of ∆DIC. To estimate ∆DIC, the forcing-response relationships of seven 25 generic deglacial carbon cycle mechanisms in regard to LGM-PI changes in four observational targets (pCO2 atm , δ 13 C atm , marine δ 13 C of DIC, and deep equatorial Pacific CO3 2-) and DIC are investigated with the Bern3D. The seven generic forcings comprise changes in wind stress, air-sea gas exchange, the export rain ratio between organic and inorganic carbon, the remineralization depth of organic particles, coral growth, weathering of organic material, and carbon stocks in the land biosphere. These seven generic deglacial carbon cycle mechanisms were varied 30 individually by systematic parameter variations in addition to well-established forcings such as orbital parameters, greenhouse gas radiative forcing, land ice albedo, coral reef regrowth, and North Atlantic freshwater forcing. In a next step, Latin hypercube parameter sampling was used to vary the processes in combination and probe for nonlinear interactions and use the results to adjust the above forcing-response relationships. An emulator of the form 35 is derived from the forcing response relationships. (∆ ) is the sensitivity for each target T to the change in each mechanism ∆ . is the offset and the slope of the respective linear fit from the multi-parameter adjustment in order to correct for non-linearities in the addition of the seven generic mechanisms. We use the same half a million parameter combinations as Jeltsch-Thömmes et al. (2019) and investigate ∆DIC. The four proxy targets (pCO2 atm , δ 13 C atm , marine δ 13 C of DIC, and deep equatorial Pacific CO3 2-) are used as constraints.

Results from the Bern3D model
Applying single constraints only, yields similar ranges for ∆DIC for each constraint, covering the range of ∆DIC 5 values from about 0 Gt C to 6000 Gt C (Fig. 7). Considering all four proxy targets simultaneously shifts the estimate of ∆DIC to higher values such that the covered range amounts to about 2400 to 5500 Gt C and a median estimate of about 3900 Gt C (±1σ range of [3350,4480] Gt C; Fig. 7). This might seem large at first. Considering changes in the different carbon reservoirs that changed over the deglacial and eventually drew carbon from the ocean (atmosphere, land biosphere, coral reef regrowth, sedimentation-weathering imbalances) helps to put the 10 estimated ∆DIC into context (see also Fig. 7). The atmospheric increase over the deglacial is well constrained from ice-cores and amounts to about 200 Gt C. The change in the land biosphere carbon storage is more uncertain. A recent estimate, which integrates multiple proxy-data, derives an increase of the land biosphere over the deglacial of about 850 Gt C (median estimate) (Jeltsch-Thömmes et al., 2019), in good accordance with another recent estimate based on a dynamic global vegetation model (Müller and Joos, 2020). The growth of coral reefs during 15 sea level rise removed an additional 380 Gt C (Vecsei and Berger, 2004) to 1200 Gt C and more (Milliman, 1993;Kleypas, 1997;Ridgwell et al., 2003) from the ocean. Finally, imbalances in the sedimentation-weathering cycle of carbon need to be considered on the multi-millennial timescales. The contribution from sedimentationweathering imbalances to ∆DIC is uncertain and covers a range from a removal of almost 0 Gt C to more than 3000 Gt C from the ocean as estimated with the Bern3D model. The results thus point to an important role and 20 large contribution from sedimentation-weathering imbalances to ∆DIC estimates over glacial/interglacial timescales, however, with a large uncertainty. Finally, it has to be noted that in order to use the cost-efficient emulator and explore a large parameter space only the change between LGM and PI was considered. Including the CH and JS supervised the project, and provided the resources and funding acquisition needed for this work. 10 Competing interests. The authors declare that they have no conflict of interest.
Disclaimer. This article reflects only the authors' view -the funding agencies as well as their executive agencies are not responsible for any use that may be made of the information that the article contains.