Ocean carbon cycling during the past 130 000 years – a pilot study on inverse palaeoclimate record modelling

What role did changes in marine carbon cycle processes and calcareous organisms play in glacial–interglacial variation in atmospheric pCO2? In order to answer this question, we explore results from an ocean biogeochemical general circulation model. We attempt to systematically reconcile model results with time-dependent sediment core data from the observations. For this purpose, we fit simulated sensitivities of oceanic tracer concentrations to changes in governing carbon cycle parameters to measured sediment core data. We assume that the time variation in the governing carbon cycle parameters follows the general pattern of the glacial–interglacial deuterium anomaly. Our analysis provides an independent estimate of a maximum mean sea surface temperature drawdown of about 5 C and a maximum outgassing of the land biosphere by about 430 Pg C at the Last Glacial Maximum as compared to pre-industrial times. The overall fit of modelled palaeoclimate tracers to observations, however, remains quite weak, indicating the potential of more detailed modelling studies to fully exploit the information stored in the palaeoclimatic archive. This study confirms the hypothesis that a decline in ocean temperature and a more efficient biological carbon pump in combination with changes in ocean circulation are the key factors for explaining the glacial CO2 drawdown. The analysis suggests that potential changes in the export rain ratio POC : CaCO3 may not have a substantial imprint on the palaeoclimatic archive. The use of the last glacial as an inverted analogue to potential ocean acidification impacts thus may be quite limited. A strong decrease in CaCO3 export production could potentially contribute to the glacial CO2 decline in the atmosphere, but this remains hypothetical.


Introduction
Since the mid-1980's the drawdown of the atmospheric CO 2 mixing ration from 280-290 ppm to 180 ppm is known from analysis of Antarctic ice core measurements (Neftel et al. (1982), Barnola et al. (1987), Siegenthaler et al. (2005)). Despite a series of attempts to explain this major fluctuation in the key atmospheric greenhouse gas CO 2 , an overall closed explanation is still missing (Broecker and Peng (1986); Heinze et al. (1991); Heinze and Hasselmann (1993); Archer and Maier-Reimer 20 (1994); Sigman and Boyle (2000); Archer et al. (2000); Brovkin et al. (2007)). An understanding of the positive feedback between climate and atmospheric CO 2 content (cf. Shakun et al. (2012)) is important also for quantifying the changing carbon cycle today. The glacial-interglacial CO 2 variations are mainly caused by a redistribution of carbon between the Earth system reservoirs ocean, atmosphere, and land. In contrast, the modern carbon cycle is perturbed by the addition of CO 2 to the Earth's surface reservoirs, which had been excluded from interaction with biogeochemical processes in these reservoirs for a long time. Nevertheless, a reliable reconstruction of glacial-interglacial carbon dynamics with Earth system models would enhance their predictability of future carbon cycle changes. The human-induced CO 2 emissions from fossil-fuel burning, 5 cement manufacturing, and land-use change have risen the atmospheric CO 2 mixing ratio to 400 ppm relative to a pre-industrial level of 278 ppm (Le Quéré et al. (2014)). The atmospheric CO 2 mixing ration would be considerably higher, if the ocean would not take up a part of the human-made excess CO 2 from the atmosphere, primarily through physico-chemical buffering (Bolin and Eriksson (1959); Maier-Reimer and Hasselmann (1987); Joos et al. (2013)). This marine CO 2 buffering has the side effect of reducing the seawater pH (Haugan and Drange (1996); Caldeira and Wickett (2003)). The impact of this progressing 10 ocean acidification on marine biota, ecosystems, and biogeochemical processes is subject to broad interdisciplinary research and harmful effects could be identified (Raven (2005); Gattuso and Hansson (2011)). Especially a reduction in the formation of calcium carbonate (CaCO 3 ) shell material is at the focus of ocean acidification impact research (Riebesell et al. (2007); Iglesias-Rodriguez et al. (2008)). Several potential geologic analogues for oceaan acidification have been suggested (during the Paleocene Eocene Thermal Maximum at 56 Ma, Zachos et al. (2005); CO 2 from the Siberian Traps Large Igneous Province, 15 Svensen et al. (2009); methane release from magmatic intrusion into coal beds, Retallack (2008)). E.g., lack of sedimentary information prior to 180 million years before present limit the information about ecosystems back in deep time. Another prominent example of significant carbon fluctuations is the last glacial interglacial cycle. For example, the lower CO 2 partial pressure at the last glacial maximum (LGM) around 21 kyr BP (21,000 years before present) -could be used as a reverse paleoanalogue to a high CO 2 world. Could we say anything about the average behaviour of marine biota with respect to CaCO 3 20 formation under glacial conditions in order to learn about possible future developments? Evidence for a better preservation of foraminifera shells is provided by Barker and Elderfield (2002) which would be in line with higher glacial carbonate ion concentrations in surface waters Yu et al. (2014).
In this study, we combine a coarse resolution biogeochemical ocean general circulation model (BOGCM) together with long sediment core record time series in order to determine the most likely change in ocean carbon cycle parameters during the 25 past climatic cycle. We include also the "rain ratio" CaCO 3 :Corg (carbon atoms bound to CaCO 3 shell material versus carbon atoms bound to organic carbon in biogenic particulate matter during biological production) in the analysis. Thus, information is provided on whether the CaCO 3 formation in the glacial low CO 2 world increased or not. A CaCO 3 production increase at the LGM would be consistent with a decrease of CaCO 3 formation in a high CO 2 world. The model simulates a time dependent sediment record at each model grid point for variations of single parameters over the past 130 kyr (see section 2). We determine 30 the resulting model sensitivities for the sediment core tracers in relation to a number of key governing parameters of the ocean carbon cycle (see section 3). The sensitivities are projected onto a linear response model. This simplified model is then used for a simultaneous fit of all carbon cycle parameters to available sediment core data from the observed paleo-climatic archive (section 4). We pursue a multi-tracer approach by considering data on marine CaCO 3 wt-%, BSi wt-% (BSi is biogenic silica or opal, rf. Ragueneau et al. (2000)), δ 13 C planktonic , δ 13 C benthic , and the atmospheric CO 2 record of the Vostok ice core in feasible limits ; ; Heinze et al. (2003); Heinze et al. (2009)). This biogeochemical ocean model requires as input three-dimensional velocity and thermohaline fields as well as ice cover data. For the spin-up of the model, a velocity field representing the pre-industrial conditions ("interglacial first guess", Winguth et al. (1999)) with an annually averaged circulation from the Large Scale Geostrophic dynamical ocean general circulation model (Maier-Reimer et al. (1993)) is employed. The effect of deep convective mixing at high latitudes is considered in this forcing, 10 which is used for transporting the passive tracer substances within the model water column reservoir. The model has a horizontal resolution of 3.5 o ×3.5 o while the vertical resolution decreases downward with 11 layers centred at 25,75,150,250,450,700,1000,2000,3000,4000, and 5000 m. The bioturbated top sediment zone, the "sediment mixed layer" is structured into 10 layers which are separated by "downcore" interfaces at 0, 0.3, 0.6, 1.1, 1.6, 2.1, 3.1, 4.1, 5.1, 7.55, and 10 cm. The simplifying assumption is made that no pore water reactions take place below 10 cm depth in the sediment (following Smith and Rabouille 15 (2002); Boudreau (1997)). HAMOCC2s parameterises the processes air-sea gas exchange, biogenic particle export production out of the ocean surface layer, particle flux through the water column and particle degradation by remineralisation and re-dissolution, advection as well as mixing dissolved substances with the ocean velocity field, deposition of particulate matter on the sea floor, pore water chemistry and diffusion, advection of solid sediment components, bioturbation, and sediment accumulation (export out of the sediment mixed layer or "burial"). The general concept of the sediment model closely follows 20 Archer et al. (1993). For the present model configuration, the following tracer concentrations are used as prognostic variables: the atmosphere includes 12 CO 2 (carbon dioxide), δ 13 CO 2 , and O 2 ; the water column contains DIC (dissolved inorganic carbon), POC (particulate organic carbon), DOC (dissolved organic carbon), CaCO 3 (calcium carbonate or particulate inorganic carbon) of 12 C and before, DIC and TAlk are used as "master tracers" from which the other derived inorganic carbon species such as the CO 2− 3 concentration, the carbonate saturation, and the pH value are computed through a Newton-Raphson algorithm. A further extension of the model as presented in Heinze et al. (2009) is the introduction of separate fields for POP (particulate organic phosphorus), DOP (dissolved organic phosphorus), and solid organic phosphorous in the sediment mixed layer to allow for variable stoichiometry ("Redfield ratios" C:P) in this study. The other ocean water column and pore water processes are represented 5 through parameterisations as described for in earlier model versions 2.1 Transport of age information and passive tracers in the sediment Next to the various solid sediment components, the model predicts also the age of each sediment component. The age of the amount of each component rained onto the model seafloor is set to the actual model time step. The age information is then handed down into the model sediment layer as a passive tracer and this information transport is also depending on the pore water 10 dissolution reactions, the deposition flux, the vertical sediment advection, and the bioturbation. Each solid sediment species has an own age tracer. The higher the sediment accumulation rate is the closer the age values of the four different sediment component ages are. For a schematic of the sediment model, please, see Heinze et al. (2009) (Figures 1 and 2 therein). The model sediment model is initialised at age 0 with only clay sediment layers and then spun-up. During the progressing integration, the model builds up its own biogenic sediment and clay sediment according to the dust deposition and the kinetic 15 control of the continental inputs of nutrients, carbon, and alkalinity as prescribed through bulk numbers. Simultaneously with the sediment distribution, the age structure of each solid weight fraction is established. In equilibrium, the global sediment accumulation rates are a function of the continental input rates. The local sediment accumulation rates and sediment compositions are fully prognostic variables. In order to also "flush" the sediment mixed layer in low accumulation regions of the world ocean, the model is spun-up for 100,000 years to achieve overall quasi-equilibrium (as documented by the stability of 20 the atmospheric pCO 2 level and the asymptotic approach of the simulated globally integrated sediment accumulation rates to the globally prescribed continental matter input fluxes).
The model includes the cycling of 13 C. The observed foraminiferal CaCO 3 sediment record reflects δ 13 C variations from both the surface and the deep waters (depending on the shell material analysis of the different surface and deep dwelling species). As an extension to Heinze et al. (2009), the present model version also includes sedimentary benthic δ 13 C. Normally, 25 only the surface ocean δ 13 C, i.e. the δ 13 C planktonic is included in the simulated CaCO 3 deposited onto the model top sediment layer and then treated subsequently during pore water reactions and sediment advection. The reason for this is that now explicit benthic production of foraminifera is realised in the model. Applying the basic concept of the passive age tracer treatment in the early diagenesis module, we also transport the δ 13 C benthic signal of the CaCO 3 fraction within the sediment mixed layer (where we set the benthic solid δ 13 C value for the depositing CaCO 3 material to the respective bottom water δ 13 C of the model 30 layer directly over the respective top sediment box). Thus, the following sedimentary variables included in the present model for direct comparison with values from sediment core analysis: the "foraminiferal" values of δ 13 C planktonic , δ 13 C benthic (as the latter would be recorded in reality, e.g., by the foraminiferal species Cibicidoides wuellerstorfi, cf. Zahn et al. (1986)), and also their vertical difference ∆δ 13 C=δ 13 C planktonic -δ 13 C benthic as an indicator of the strength of the biological pump (or the vertical nutrient gradient in the water column, see Shackleton and Pisias (1985). The model master tracers DIC and TAlk for the inorganic carbon cycle include the sum of 12 C and 13 C. Explicit 13 C concentrations are computed from the spin-up through the requirement that the pre-industrial atmospheric δ 13 C value of -6.5 ‰. The continental input rate of 13 C and the corresponding output rate through sediment accumulation where iteratively determined, so that in the final model spin-up the 12 C and 13 C are fully compatible with the pre-industrial δ 13 C value.

Control run
Important global bulk numbers for the quasi-equilibrium state of the 100,000 yr long control run are given in Table 1. For illustration of the control performance, examples for the export production rates and the respective sediment components are provided in Figure 1. Meridional cross sections for dissolved inorganic phosphate, the carbonate ion concentration, and δ 13 C of total dissolved inorganic carbon are shown in Figure 2. The model is able to reproduce the major characteristics of global 10 tracer distributions, with decreasing deep water carbonate ion concentrations and δ 13 values as well as increasing phosphate concentrations from the Atlantic to the Pacific Ocean. Simulated tracer patterns appear to be more smoothed as compared to the observations because of the coarse spatial resolution of the model. The CaCO 3 sediment distribution in the northern Pacific Ocean may show somewhat too high CaCO 3 wt-% as compared with observations. Based on experience with previous simulations, the CaCO 3 content of the sediment mixed layer appears to react quite sensitively to the degree of undersaturation 15 and the redissolution rate constants chosen. Therefore, a slight overestimate of CaCO 3 sediment in the control run avoids for many grid points a potential complete dissolution of CaCO 3 during the sensitivity experiments.

Sensitivity experiments
In order to determine, how the governing process parameters of the carbon cycle changed from the cold low CO 2 world to the warmer and higher CO 2 pre-industrial world, we first carried out a suite of sensitivity experiments with the full 3-D 20 carbon cycle model. In each run, we change only one of the governing parameters in order to determine the respective change of the biogeochemical system. We are in particular interested in the respective change of the inorganic as well as organic carbon cycle and how these changes are imprinted onto the simulated paleo climatic record. In a previous study Heinze et al. (1991) had carried out sensitivity experiments with an earlier version of HAMOCC and also compared the results with data from the observational paleo-climatic archive. The present study, however, differs significantly from the previous approach. In 25 Heinze et al. (1991) only equilibrium responses of the biogeochemical state to instantaneous parameter changes (permanent switch from one of the control run configurations) were tested, the alternative velocity field used did not resemble the glacial ocean very well, and the model did not include an interactive early diagenesis sub-model (only interactive bulk sediment reservoirs which did not allow for varying total system inventories). The latter issue thus allowed only a comparison between the modelled tracer values and sediment core derived "proxy data". Also the fixed total system inventory of carbon, nutrients, 30 and alkalinity allowed a limited range of system changes. Already in the first study of HAMOCC including an interactive 5 Clim. Past Discuss., doi: 10.5194/cp-2016-35, 2016 Manuscript under review for journal Clim. Past Published: 18 April 2016 c Author(s) 2016. CC-BY 3.0 License. diagenesis module by Archer and Maier-Reimer (1994) revealed the higher sensitivity of such a model configuration, as the water column tracer inventory could vary considerably more dynamically.
We first carried out a control run over the past 130,000 years re-starting from the standard spin-up of the model. This control run revealed the excellent mass conservation and also the almost drift free 13 C distributions in all reservoirs. Thereafter, a series of simulations with perturbed model parameters were conducted by restarting from the same re-industrial control 5 run. The runs were started from "Eemian" conditions which we assumed to be similar to first order to the pre-industrial state as resulting from the HAMOCC spin-up. We then scaled the various parameter perturbations with the δD temperature record of the EPOCA Dome C ice core project (Jouzel (2004)) assuming that all carbon cycle variations during the past climatic cycle are correlated with temperature or the rather similar temporal pattern of the atmospheric CO 2 variation. There are uncertainties associated with this assumption, as the Antarctic air temperature record is a localised record and may not 10 reflect global temperature change, and also because ocean temperatures may lag behind this signal. On the other hand, the local Antarctic δD temperature record is highly correlated with the global atmospheric pCO 2 signal (Siegenthaler et al. (2005); Barnola et al. (1987); Jouzel et al. (1987)). The seawater temperature change would probably not lag too long behind this signal if reasonably long time intervals are considered. Due to the close correlation between air temperature and atmospheric pCO 2 over the past glacial-interglacial cycle, we cannot strictly discriminate between drivers for ocean and land carbon cycle 15 changes coming from physical (temperature) and biogeochemical forcings (pCO 2 , carbonate saturation, pH value), but have to justify this independently for the different parameters under investigation. For the sensitivity experiments, we selected the following governing model parameters of the carbon cycle: (1) The sea surface temperature for computation of the chemical and biological constants, (2) the release/uptake of carbon from/by the terrestrial biosphere, (3) the degradation rate constant of POC in the water column, (4) the dissolution rate constant of BSi in the water column, (5) the export production rain ratio 20 CaCO 3 :POC, (6) the 3-D oceanic velocity field, (6) the glacial dust deposition and associated stimulation of biological export production, and (8) the Redfield ratio C:P.
These parameter changes are summarised in Figure 3 and Table 2 and further described more in detail below.

Sea surface temperature for computation of the chemical and biological constants
In this experiment, we reduced the sea surface temperature up to a maximum change of 5 K (with the consideration that the 25 minimum temperature stays at or above the freezing point of sea water) only for the computation of the temperature dependent chemical and biological parameters in the model, which usually all have their fixed control run values. The effect on the CO 2 solubility is relatively strong accounting for ca. 12 ppm change in atmospheric pCO 2 for a 1 K change in sea water temperature (cf. Broecker and Peng (1986)). The 5 K change in temperature is approximately the maximum change that can be expected on the basis of stable oxygen isotopes or other temperature indicators (e.g., Stute et al. (1995)). 3.2 Release/uptake of carbon from/by the terrestrial biosphere Also variations in the land carbon cycle need to be properly taken into account for the overall pCO 2 signal in the atmosphere and the imprint of 13 C on the sedimentary record. We do not explicitly model the terrestrial carbon cycle in our experiment. Rather 6 Clim. Past Discuss., doi: 10.5194/cp-2016-35, 2016 Manuscript under review for journal Clim. Past Published: 18 April 2016 c Author(s) 2016. CC-BY 3.0 License.
we assume here, that the land biosphere would have a smaller net primary production (Hoogakker et al. (2016)) and standing biomass stock and hence provide a release of carbon to the atmosphere under colder and dryer conditions. The amplitude for biomass carbon loss from the land biosphere to the atmosphere was set to -500 PgC, which is on the more modest side of respective estimates from terrestrial paleo-climate records according to Crowley (1995). The reason for this choice is, that a potentially increased storage of organic carbon under the ice sheets and at cold conditions (lower bacterial activity for plant 5 biomass degradation) is also being discussed (e.g., Zeng (2003)). Such storage could lower the overall carbon loss from the terrestrial biosphere under glacial conditions.

Dissolution rate constants of POC and BSi
We explored the temperature effect on both sinking particulate organic carbon as well as biogenic silica. We assumed that the degradation of both substances slows down under lower temperatures due to thermodynamic effects as well as decreasing bac-10 terial decomposition (Turley and Mackie (1995); Van Cappellen et al. (2002), Bidle and Azam (1999)). Through the respective parameter changes we tested separately a strengthening of the biological carbon and of the silicon pump which would have led to a deeper penetration of carbon and nutrients into the water column under glacial conditions ("vertical fractionation" as described by Boyle (1988a) and Boyle (1988b)). We chose a maximum decrease of degradation strength of 10%.
3.4 Export production rain ratio CaCO 3 :POC 15 A decrease in rain ratio has been suggested as one of the potential mechanisms to extract carbon from the atmosphere during glacial times (Berger and Keir (1984)). A reduction of the CaCO 3 production to zero could potentially explain half of the glacial pCO 2 drawdown from the atmosphere (Broecker and Peng (1986)). However, the experiments by Zondervan et al. (2001) and Riebesell et al. (2000) would suggest the opposite: In a low CO 2 world, the CaCO 3 :POC ratio would be expected to decrease.
This would make it more difficult to explain the glacial CO 2 uptake by the oceans. Heinze and Hasselmann (1993) showed 20 that this in principle is possible. We follow in our experiment here, however, the idea of an diminished CaCO 3 :POC export production rain ratio during glacial times. In order to get a pronounced signal, we decrease the amplitude to a -10% change at the LGM with respect to the control run.

3-D oceanic velocity field
Practically all geochemically relevant paleoclimate tracers in the ocean depend on the ocean circulation. Therefore, it is nec-25 essary to take a realistic glacial velocity field into account in order to provide the correct addition to tracer changes which may be induced by biological and chemical processes. We follow here the same approach as outlined in Heinze (2001), where we "mix" the velocity, thermohaline, convective mixing, as well as sea ice margin values from two different runs of the dynamical ocean general circulation model LSG (Large Scale Geostrophic Model, Maier-Reimer et al. (1993)). These runs are the respective interglacial and glacial first guess runs from Winguth et al. (1999) which show major features of the pre-industrial as 30 well as LGM distributions of the δ 13 C of dissolved inorganic carbon. As maximum amplitude for the velocity field change we assume a full switch to LGM conditions at the largest negative excursion of the δD EPICA Dome C ice core curve. It would be preferable to have respective atmospheric forcing fields available to have physically better constrained velocity fields at hand.
However, the considerable effort to realise such a run (which would need to be carried out with the full seasonally resolving LSG model) was deemed too large for our study here. The "kinematic" approach chosen here, should work for practical reasons as the respective δ 13 C distributions of DIC show realistic patterns.

Glacial dust deposition
We carried out a number of sensitivity tests concerning the effect of a change in continental dust deposition onto the sea surface. First we only tested the pure effect of clay addition through dust by applying the glacial dust deposition field by Mahowald et al. (1999) at the LGM and respective mixtures of interglacial and glacial dust deposition fields at other time intervals. As the dust records of marine sediment core records usually do not follow the pattern of the δD record (e.g., Figure 4 10 of Carter and Manighetti (2006)), we scaled the change in dust supply with the 6 th order of the glacial-interglacial difference of the δD curve, i.e. assumed a sudden increase of dustiness only at really low temperatures. Again we chose a maximum value of dust change to 100% glacial conditions at the LGM and smaller changes else (see Figure 3). The effect of this dust addition would be a change in the dilution of the sediment weight percentages of the reactive sedimentary material by inert clay and also a change in the local sediment accumulation rates. It had been suggested that a stimulation of biological export 15 production through increased dust and iron deposition to the ocean surface under glacial conditions (e.g., Martin et al. (1994); Berger and Wefer (1991)) could occur. In an additional experiment, the maximum uptake velocities for phosphorus, carbon, and silica were scaled with the glacial-interglacial difference of the dust deposition as provided by Mahowald et al. (1999). For the results of the fitting procedure as shown below, we considered here the run with the dilution effect only in order to separate sediment dilution effects from those caused by an increase in biological production. 20 3.7 Redfield ratio C:P Finally we allowed carbon over-or underconsumption in response to the interglacial to glacial change in environmental conditions. Riebesell et al. (2007) and Bellerby et al. (2008) postulated an increasing carbon overconsumption C:N under high CO 2 conditions resulting from a mesocosm experiment. In contrast, reversely, a change towards increasing carbon overconsumption during the glacial low CO 2 world had been suggested as a powerful mechanism to account for the glacial CO 2 25 drawdown and at the same time cause a plausible foraminiferal δ 13 C signal (e.g., Shackleton and Pisias (1985);Broecker (1982); Broecker and Peng (1986); Heinze and Hasselmann (1993)). In our sensitivity experiment we investigated the possibility for an increase in carbon overconsumption by imposing homogeneously over the ocean an increase in the Redfield ratio C:P. We specified a maximum change by 15% at the LGM (see Figure 3, lowermost panel).
We try to make a reasonable trade-off between the number of free parameters and the need for inclusion of the most essential 30 parameters which contribute to the tracer variations in the paleo-record. Thus we aim at isolating imprints of parameter variations in the observed sedimentary record which would otherwise be masked by independent further processes (e.g., it could be 8 Clim. Past Discuss., doi: 10.5194/cp-2016-35, 2016 Manuscript under review for journal Clim. Past Published: 18 April 2016 c Author(s) 2016. CC-BY 3.0 License. futile to derive changes in biological pump strength from the foraminiferal δ 13 C distribution in the ocean without taking into account respective changes in the ocean currents which also affect the marine δ 13 C distribution).
For determining the most likely choice of the governing parameters, we analysed modelled output for a series of paleoclimatic data sets. The model delivers first of all prognostic atmospheric pCO 2 values for each sensitivity experiment relative to the control run during each time step for comparison with the Antarctic ice core records. The model runs were started 5 from year 129536 BP and integrated until the present (the start year was chosen so that the first and last δD value are close to each other, no anthropogenic CO 2 emissions were imposed during the very last part of the integration). Each run requires around 4 days CPU time on an advanced UNIX work station. Using the methodology of Heinze et al. (2009) the model at each grid point accumulates sediment in a continuous way for the components CaCO 3 , opal, organic carbon, clay, and organic phosphorus. δ 13 C planktonic and δ 13 C benthic of the CaCO 3 fraction and the individual ages of each sediment component plus 10 the atmospheric CO 2 partial pressure were stored together with the actual accumulation rates at the base of the sediment mixed layer. In order to reduce the amount of data output, we only stored 100-yr averages. Over the integration time period, for each run a continuous temporally varying sediment stratigraphy is built up. After the model runs, we "drill" into this synthetic sediment and "recover" simulated sediment cores for comparison with the observations from real world sediment core data.
We can provide such modelled cores at every grid point of the model individually. 15 4 Linear response model, observational sediment core data base, and fitting procedure The inverse modelling procedure as applied and described here, draws in many aspects from the work by Heinze and Hasselmann (1993). We partly repeat some of the methodological issues here, to avoid that the reader has to refer too often to a separate study when reading this analysis here. The method is illustrated in Figure 4. We had carried out a total of eight sensitivity experiments with the full 3-D biogeochemical ocean general circulation model, providing 8 different data sets for respective 20 changes in 79 different tracer time series (see below) resulting formally in a total of 20,540 modelled tracer values and the same number of respective observed values (for each run). The eight sensitivity experiments represent the response of the full 3-D model for variations of n=8 governing input parameters x j (j=1, . . ., n) ( Table 2) resulting in m=20,540 paleoclimate tracer record changes g i (i=1, . . ., m) as induced in these experiments as a consequence to the changes x j of the parameters.
The resulting linear response model thus consists of the matrix describing the linearised relation between perturbations 25 of the input and output variables. The model output variablesĝ i in the sensitivity experiments correspond to the observed perturbations of the atmospheric CO 2 record and the sediment records of CaCO 3 , δ 13 C benthic , δ 13 C planktonic , and opal (BSi) (see following sub-section).
In the full non-linear biogeochemical ocean general circulation model, the paleo-climate tracers G i are functions of the governing model parameters X j , Clim. Past Discuss., doi: 10.5194/cp-2016-35, 2016 Manuscript under review for journal Clim. In the simplified linear response model, the relation (1) is Taylor expanded (and then truncated) about a reference state X j = X 0 j , G i = G 0 i as given by the defined by the standard run (with zero parameter changes): The response coefficients A ij form the elements of the model matrix A. These response coefficients resulted from the model sensitivity experiments under the assumption of a linear relation 5 between the parameter changes x j and the response vectorĝ i predicted in every sensitivity experiment j.

Data base of observations from the paleo climate record
For calibrating the free model parameters, we employ a data base of paleo-climatic records as summarised in Tables 3 and   4. Locations are indicated in Figure 5. Most of the δ 13 C benthic and δ 13 C planktonic data were taken from the compilation of Oliver et al. (2010). The majority of CaCO 3 wt-% data were taken from a compilation by Hoogakker (pers. comm., Table   10 4). Further marine sediment core data (δ 13 C benthic , δ 13 C planktonic , BSi, and CaCO 3 wt-% ) were taken from various literature sources (see references in Tables 3 and 4). For atmospheric pCO 2 we use the Vostok ice core signal from Antarctica (Barnola et al. (1987)). We did not explicitly synchronise the different paleo-climatic time series to a common age model here, but rather took the measurement/age combinations at face value. We estimate the error due to this to up to a few thousand years (compare the discussion in Oliver et al. (2010)). As we are not interested here in exact timing and phase shifts of signals in dif-15 ferent variables, and due to the overall errors in the modelling approach (coarse resolution model, only crude representation of circulation changes, assumption that pre-industrial model state corresponds also to the Eemian) we assume that this approach is reasonable. Some sediment core data even were extracted from published graphics by hand. Most of the δ 13 C benthic and CaCO 3 wt-% time series employed, however, were synchronised in respective data compilations. For the fitting procedure, we interpolated all observed data to regular time series with an increment of 500 years (260 equidistant data points in time). 20 Modelled annual global mean values of atmospheric pCO 2 were compared for comparison with the Vostok pCO 2 signal.
The corresponding model time series data of the different sediment paleoclimate tracer curves where extracted for the various model sensitivity experiment runs from the respective model grid point closest to the location of real sediment core extraction.
Modelled sediment data were interpolated with respect to age onto the same 500 years points as the observations. The age of each sediment component was simulated according to Heinze et al. (2009) during the sensitivity experiments and stored 25 together with the simulated depth in core of each modelled sediment variable.
For the fitting procedure, observed and modelled data where translated to changes relative to the preindustrial. Thus, we analysed the changes in the tracers relative to this reference while the respective parameter changes were taken as changes relative to the control run values. 30 Our general fitting procedure closely follows the method as described in Heinze and Hasselmann (1993)).

The general minimisation problem
We fitted the linearly modelled paleoclimatic tracer changesĝ i to the observed tracer changes g i through minimizing the mean square sum Q 2 of the individual tracer errors ǫ i =ĝ i − g i : The dimensions of the vectors and matrices are described by the following symbols (cf. Heinze and Hasselmann (1993)): The left marker of a matrix symbol denotes the number of rows, the right marker specifies the number of columns. In our specific case, the marker "·" corresponds to 1, "|" to the number m = 20, 540 of observed tracer data, and "|" to the number n = 8 of parameter changes x j to be estimated. The solution |x· which minimizes the mean square expression Q 2 is derived through the necessary condition dQ 2 dx = 0 : which also can be written as Through one realisation of the fitting procedure we thus, simultaneously, determined all eight optimal parameter changesx j (t) necessary for a best possible reproduction of the observed tracer changes. In addition the procedure yields a time series of 15 residual errors ǫ i (t) which quantify the expected remaining discrepancies between modelled and observed tracers after the optimisation (we write here "expected" as the non-linear model may result in somewhat different residuals than the linear response model). Our problem is formally overdetermined where m = 20540 tracer data points must be fitted with must be fitted with n = 8 parameters. Therefore, the least squares formulation of our solution procedure is appropriate.
The error expression of eq. (3) is defined with a simple unit matrix norm. The more general maximum likelihood norm 20 would be given by the inverse of the covariance matrix of the measurement errors (e.g., Martin (1971)). In eq. (3), we make the assumption that the errors are uncorrelated and have the same variance. At present, we have no conclusive means for estimating the error correlation. Therefore, the assumption of a simple unit matrix norm is reasonable here. However, we normalised the variables g i by the mean absolute values of the respective tracer over the past climatic cycle, so that differential weighting of different tracer records, e.g., already from the use of different physical units is smoothed out.

The SVD least squares solution
In many inverse problems it is not possible to determine a physically meaningful least squares solution directly from eq. (4). The solution can be contaminated by noise if the model matrix |A| is badly conditioned. This leads in practice to unrealistically high variations of the optimised parameter changes if these depend only on very small variations in the tracer changes. Beforehand, it cannot be seen whether A is badly conditioned or not. Therefore, we applied the singular value decomposition (SVD) technique which provides a quantitative treatment of the noise problem. For an introduction to this technique, please, see the very useful work by Wunsch (1989). The model matrix A is decomposed into a product of three matrices (SVD; see Lanczos (1961)): The column vectors of |U | and |V | represent orthonormal bases |u i · , i = 1, . . . , m and |v j · , j = 1, . . . , n for the tracer (here 5 tracer change) and parameter (here parameter change) spaces, respectively, associated with the model matrix |A| . |Λ| The diagonal rectangular matrix consists of the singular values λ k . The eigenvectors |u k ·, |v k · and the eigenvalues λ k form the solution of the coupled eigenvalue problem: with the following properties: 15 The number p ≤ min(m, n) of non-zero singular values defines the rank of |A| .
The unknown parameter changes |x· and the modelled tracer changes |ĝ· can now be written in terms of the linearly independent eigenvectors |v j · and |u i · , respectively: Clim. Past Discuss., doi: 10.5194/cp-2016-35, 2016 Manuscript under review for journal Clim. where the α j are unknown and β i = ·u T i | |g·. The observed tracer changes are correspondingly given by |g· ≡ m i=1 β i |u i ·.) Equations (4), (7), and (9) together yield the solution for the optimal parameter changes: with α i = β i /λ i . Because only p = n < m vectors |u i · are combined in the overdetermined case, residuals between modelled tracer changes |ĝ· and observed tracer changes |g· occur. Thus, the observed tracer variations cannot be completely 5 recombined by the least squares solution of eq. (12).
In the exactly determined case we would have where δ ij = 1 f or i=j 0 f or i =j . For the overdetermined case we arrive at with |R I | = |Û | |Û T | being the "resolution matrix" (Wiggins (1972); Wunsch (1989)) of the changes in the tracer concentrations and |Û | = ( |u 1 · , . . . , |u n · ) is the truncated "model" version of the tracer base matrix |U | .
According to eq. (12), the coefficient α i in the linear combination of the parameter space eigenvectors |v i · grows to very large values if λ i becomes very small. These very small λ i values lead to an unstable solution: Only unreasonable extremely 15 large changes in selected parameters would be able to create the tracer changes as observed. The contributions to the solution from these small eigenvalues should accordingly be filtered out or considered with a significantly reduced weight.
For the resulting parameter changes from our fitting procedure, we discarded those components in eq. (12) associated with the smallest singular values step by step starting with the smallest singular value, followed by the next smallest one , and so on) to filter out unrealistically large parameter changes. The associated approximate solution becomes then The number p of eigenvalues retained is called the effective rank of |A| . In our case here the full rank solution would be achieved for a rank=8 (including all parameter changes). For an effective rank smaller than min(m, n), the formally overdetermined problem is now changed into an underdetermined problem. Nevertheless solution eq. (15) is the unique least squares solution with minimum length of the solution vector for the underdetermined case (e.g., Matsu'ura and Hirata (1982)). Such combinations. The full rank solution (p = n ≤ m) would be: The solution for the rank deficient case (p < n ≤ m) is now: 5 |R P | = |V ≀ ≀V T | is the "resolution matrix" of the parameters with |V ≀ = (|v 1 · , . . . , |v p · ). This resolution matrix |R P | projects the full (in art insufficiently constrained) set of parameters onto the sub-space of stable solutions for the parameter changes.
If for fixed i the diagonal element R P ii of the matrix |R P | is close to 1 and the other elements R P ij , j = i are close to zero, parameter x i is still well resolved even in the underdetermined case. Less well resolved parameters, are instead replaced by a stable linear combination of the original parameters. Through analysis of the resolution matrix for the parameters, we can 10 therefore determine which parameter changes can be uniquely determined by the analysis and for which parameters no clear conclusion can be made.

Results and their discussion
We carried out numerous runs with the linear response model testing different combinations of parameters and observed time series of the paleoclimate as constraints. Overall we experienced, that the paleoclimatic information and the system sensitivity 15 as provided by the biogeochemical ocean model lead to quite consistent results, if several parameters and paleoclimatic tracers are taken into account.
We present here results where we included all parameter changes and all paleoclimatic tracer curves (see summaries in Table   5 and Figures 6-8). The 8 unknowns (namely the excursions of the parameter amplitudes relative to interglacial conditions) where determined by an optimal fit to the 20,540 data points form the observational records. Formally, such a system is heavily 20 overdetermined. No a priori knowledge and no artificial limits where imposed to the free parameters to be determined. Are the results for the free parameters reasonable? Figure 6 shows the parameter time series for the full rank solution (compare to the first guess values as shown in Figure 3), which also gives the best fit to the observations. Please, note that the model parameter changes by purpose show the same temporal pattern (special case dust as discussed above) and only the maximum amplitude of the parameters was determined. The a posteriori parameter changes remain relatively close to the initially chosen parameter 25 changes (first guesses), yet with one exception. The rain ratio change shows a dramatic decrease in pelagic CaCO 3 production.
Such a change may not be fully out of scope (see, e.g., the discussion in Broecker and Peng (1986)) but one would possibly expect rather an increase in CaCO 3 production at low ambient pCO 2 and high CaCO 3 saturation. The suggested strong change in the rain ratio may be an unstable part of the solution, where the fitting procedure needs huge changes in one (or possibly several parameters) to achieve minor changes in the paleoclimate tracers simulated. (Nevertheless, the corresponding SVD solution is the unique solution of the system for full rank with altogether minimum deviation from the control run.) We, therefore, started to reduce the rank from 8 to first 7 and then 6, i.e. we still aimed at determining all 8 parameter changes simultaneously, but partially in linear combination of each other. In the solution for rank 7 (Figure 7), the change in sea surface temperature seems to be unrealistically high. A global decrease by about 8 • C would imply an unrealistic widespread ocean freezing. The solution for rank 6 ( Figure 8) finally arrives at overall realistic parameter excursions during the past climatic  (2013)). Interestingly, our estimate is closest to the value originally determined by Shackleton (1977), which may be due to the marine δ 13 C constraint used in that study which also is included in our analysis. Overall, our study confirms earlier estimates on the LGM terrestrial carbon storage through an 25 independent by our approach.
Our fitting procedure confirms that the glacial ocean circulation should have somehow resembled the simulated glacial ocean conditions as deduced by Winguth et al. (1999). However, in an ideal case with a perfect simulation of the glacial ocean circulation, the switch to this circulation at glacial conditions should be close to 100 % and not only 50 % as in this study. This may be due to still some deviations of the simulated glacial ocean velocity especially for the Pacific Ocean, but also due 30 to deficiencies in the biogeochemical model including its sediment module. Nevertheless, the tendency of a reaction towards glacial physical boundary conditions resulting from the fitting procedure is encouraging and consistent.
The most difficult part in interpreting the results are the reactions of the ecosystem parameters. The small changes in POC and BSi degradation or dissolution rates may be realistic in view of the overall still modest change in seawater temperatures.
The almost vanishing rain ration change in the rank 6 solution and the suggested rain ratio decrease (corresponding to CaCO 3 export production decrease) cannot be easily explained. The lack of better information on CaCO 3 production during glacial times can here not give any clear answer concerning potential ocean acidification impacts. The reaction of the system during 5 glacial times remains diffuse in our analysis (see also below for the resolution of parameters analysis).
For some CaCO 3 tracer records a remarkable improvement in fit resulted when employing increased dust deposition and the related dilution effect (not shown). However, on the average the inverse approach does not suggest a strong effect of the dust deposition changes for a better reproduction of the paleo-climate tracers here. This can possibly be due to the simple modelling approach which we take here and regional variations in dust deposition which so far cannot be resolved by the input fields. 10 The increase in carbon overconsumption (change in stoichiometry of carbon to nutrient elements) as suggested in our study confirms earlier results (Heinze and Hasselmann (1993)). This does not necessarily mean that such a change in carbon to nutrient ratio is realistic, at least for two reasons: Ecosystem changes to achieve such a change would need to be dramatic. The glacial-interglacial atmospheric CO 2 change corresponds in amplitude to about the pre-industrial to present increase. So far, no respective carbon underconsumption due to the anthropogenic CO 2 release has been identified. Rather, the opposite has been 15 suggested for further dramatic increases in atmospheric and sea surface CO 2 (Riebesell et al. (2007)). The increase of carbon to nutrient ratios under cold and low-CO 2 conditions would contradict the evidence from mesocosm experiments in a Norwegian coastal setting (Riebesell et al. (2007)). The experimental results from the mesocosm studies may potentially be influenced by the specific experiment setting and thus may not be valid at other locations. Earlier, Heinze and Hasselmann (1993) could not well separate the effect of stimulation of POC export production by a change in stoichiometry or by an increase of the 20 oceanic inventory of dissolved phosphate. Thus, the requirement of a carbon overconsumption in this study may indicate that the biological organic carbon production may have been stimulated by other processes than carbon overconsumption , e.g., fertilisation of ocean productivity by dust-derived micro-nutrients.
Moreover, we look at the impact of the rank reduction on the parameter change quantification and the goodness of fit to the observations. The reduction of the rank of our linear system involves that not all unknown parameter changes can be determined 25 in an independent way. The resolution matrix of the parameters (cf. eq. 17 identifies which free parameters still can be uniquely determined after the rank reduction and which ones can only be given as linear combinations of each other. Instead of giving the matrices in terms of numerical values, we visualised the matrices through circles in Figure 9, for the solution with rank 8 (full rank), rank 7, and rank 6. The full rank resolution diagram shows a perfect diagonal with all parameters formally perfectly resolved. For the rank reduction, the resolution for the rain ratio change, and also the BSi dissolution rate 30 change deteriorates strongly. This shows, that the linear system is in fact poorly constrained, and that after all hardly anything can be firmly said about the glacial-interglacial rain ratio change. Unfortunately, implications for future ocean acidification impacts cannot be deduced from our analysis. This is important information, however. It may be that the quite comprehensive sediment record collection which we employed here, may still be inadequate to address the ocean acidification problem and possibly also the atmospheric glacial CO 2 draw-down problem. Parts of the real system may not be represented -being it in 35

16
Clim. Past Discuss., doi: 10.5194/cp-2016-35, 2016 Manuscript under review for journal Clim. Past Published: 18 April 2016 c Author(s) 2016. CC-BY 3.0 License. the observations, in the model, or both. Further sedimentary tracer types may be needed in order to make the data base more complete and to resolve all and possible more carbon cycle parameters properly. The trade-off between realistic parameter estimates and goodness of fit becomes obvious when one checks, e.g., the overall estimated reproduction of the atmospheric Vostok ice core curve for the full rank and rank deficient solutions ( Figure 10). While the full rank solution gives a quite striking fit to the atmospheric CO 2 curve, already the rank 7 solution shows a strong decrease in fit, while finally the rank 6 5 solution (and thus also all solutions with smaller ranks) reveals an almost complete deterioration of the fit. This provides a clear dilemma.
In general, the improvement of the fit for a simultaneous optimisation of all parameters is quite weak except for atmospheric pCO 2 (see summary in Table 6). However, our 3-D biogeochemical ocean model is not detailed enough to capture all aspects of the glacial ocean correctly. For example, our representation of the circulation variations is simplified. Also it may be possible 10 that we have missed one or more key processes, which may have caused or contributed to the simultaneous changes in the paleoclimate tracer distributions. Further, we did not carry out any regional differentiation concerning the perturbation of governing parameters but used global mean changes. More detailed data assimilation schemes and the use of higher resolution coupled biogeochemical-physical models will be an option on the future for better reproducing the sediment core data, however, for the price of many more degrees of freedom to constrain correctly. It remains a possibility, that the glacial biological CaCO 3 shell 15 material production indeed was considerably smaller than the interglacial production in the case that this parameter change did not cause any major imprint on the paleo-climate tracer combination which was used in this study. In addition, also some of the sedimentary records may reflect specific local conditions which spatially cannot be resolved by our coarse model. Such special local conditions can include sediment focusing due to specific bathymetric features, smaller scale dynamic flow conditions such as localised upwellings, and meandering frontal systems. 20 We also analysed the change in simulated tracer records as compared to the control run (Table 7). As one can see, in spite of the poor fit (except for the atmospheric pCO 2 record), indeed deviations from the tracer output records of the model control run occurred. These tend to get smaller with reduction of the rank.

Conclusions
In our study we combine a comprehensive sedimentary data base from the paleo-climatic sediment core record with a coarse 25 resolution BOGCM. We assume that the governing carbon cycle parameter variations over the past climatic cycle follow the same temporal pattern. As the paleo-record of atmospheric pCO 2 and surface temperature show a strong correlation (Siegenthaler et al. (2005)), we cannot decide in general on whether the parameter variations are caused by physical forcing (temperature), chemical forcing (CO 2 ), or both simultaneously. With such approximation, we can reduce the search for unknown parameter changes to the maximum amplitude of these changes at the LGM. By employing 79 observed paleo-climatic 30 tracer data curves (with altogether 20,540 formal data points) for determination of 8 unknown maximum parameter excursions with respect to the preindustrial situation, we arrive at a formally completely overdetermined linear response model for estimating the parameter changes. We can confirm the quantification for mean global maximum SST change (-5 K) and carbon loss from terrestrial systems  to the atmosphere as suggested by earlier studies. Further, our model-data combination clearly identifies a substantial increase in the biological organic carbon export during glacial times, though the underlying process may not be a carbon to nutrient overconsumption. The response of glacial-interglacial changes in the biological CaCO 3 production remain uncertain. CaCO 3 :C org rain ratio reductions may have contributed to the glacial CO 2 reduction in the atmosphere, but such rain ratio changes cannot be resolved with our linear response model. This can be due to model deficiencies in 5 our approach, or due to lack of information from the present data base for this purpose. With our combination of model results and observations no overall answer concerning the mechanisms behind the glacial pCO 2 drawdown can be provided. Simple synchronous temperature dependent or atmospheric pCO 2 dependent forcing fields and governing parameter time series of the same temporal pattern appear to not be able to reproduce the marine sediment core record satisfactorily.  Table A1.
Air-sea gas exchange 20 The atmospheric model reservoirs atmosphere and ocean surface layer exchange CO 2 and oxygen. The flux of CO 2 across the air-sea interface is simulated as follows: where F CO2 is the carbon dioxide flux across air/sea interface and k CO2 is the specific gas exchange rate. The CO 2 partial pressure is calculated from the free carbon dioxide concentration (sum of aqueous CO 2 and carbonic acid) in seawater with 25 Henry's law through use of the solubility α (Weiss (1974)): For the oxygen flux the following approach is pursued: where F O2 is the net gas flux between sea surface and atmosphere, k O2 is the mean gas transfer velocity, C O2,equilibrium and modelled) value, f (T, S) is the measured solubility as a function of sea water temperature and salinity, and C O2,atmosphere (t) is the tropospheric O 2 concentration, where f (T, S) is taken from Weiss (1970). O 2 gas exchange is resulting from the analytical solution of the differential equation where C is the gas concentration in [mole/cm 3 ], F the flux in [mole(cm 2 · s)] of gas into (or out of) a control volume of 1 cm 5 length, 1cm width and thickness ∆z in [cm]. The numerical formulation is where ∆t is the time step (here 1 year) and ∆z is the thickness of the layer affected by gas exchange (∆z is set equal to 50 m in cases of hydrostatic stability and to the maximum depth of the convective layer else). The model atmosphere is 10 represented through a 1 layer box over each grid point. At every time step, zonal averages are determined for the atmospheric concentrations. Gas transport is simulated through meridional diffusion only (as in reality the intrahemispheric tropospheric mixing time is much shorter than the time step of 1 year applied here).
Biogenic particle export production and particle flux through the water column 15 As a consequence of the annual averaging, the model simulates only export production of biogenic particles. Particle production is assumed to exclusively take place in the model surface layer representing the euphotic zone. Phosphate serves as ultimate biolimiting nutrient, the nitrogen cycle is not simulated explicitly. POC (particulate organic carbon) and opal export production rates are predicted using Michaelis Menten kinetics for nutrient uptake (e.g., Parsons and Takahashi (1973)): where P P OC and P opal are the POC and opal export production rates [molel −1 yr −1 ], Red(C : P ) is the Redfield ratio C:P, V P OC max and V opal max are the maximum uptake rate of phosphate and silicic acid from the water column [yr −1 ], and K P OC s as well as K opal s are the respective half saturation constants. The parameters V P OC max , V opal max , K P OC s , and K opal s are prescribed as a function of sea surface temperature following Heinze et al. (2003). 25 CaCO 3 export production depends on the local production ratio P opal /P P OC . CaCO 3 export starts to gradually increase gradually when P opal /P P OC sinks below a threshold value S opal , i.e., when the silicic acid concentration is becoming too

19
Clim. Past Discuss., doi: 10.5194/cp-2016-35, 2016 Manuscript under review for journal Clim. small to support exclusively diatom growth: f or P opal /P P OC < S opal ; f or P opal /P P OC ≥ S opal 5 where R is the maximum possible rain ratio C(CaCO 3 ):C(POC) and S opal is the threshold value of P opal /P P OC for gradual onset of CaCO 3 production.
Particle fluxes and particle degradation are simulated through mass balance equations for sinking particulate matter M settle , where M settle stands for the different particle species P OC settle , CaCO 3settle , opal settle , and clay settle respectively. The general formulation for the mass balance M settle becomes then: (f or other layers) 15 where ∆z 0 is the thickness of the euphotic zone [m], w M is the particle settling velocity [m yr −1 ], r M is the reaction rate constant [yr −1 ] for degradation of particulate matter, and P M is the export production rate in the top layer, where again "M " denotes one of the particle species POC, CaCO 3 , opal, and clay respectively. P clay is the dust input from the atmosphere which is prescribed for the control run according to the modern dust deposition field from (Mahowald et al. (1999)). Clay material is considered as chemically inert. The system of equations (A9) is solved through an implicit numerical scheme. 20 CaCO 3 dissolution in the water column is performed depending on the prevailing carbonate saturation using the following rate constant: where k CaCO3 is a fixed standard rate constant, and (1 − Ω) = ([CO 2− 3 ] sat − [CO 2− 3 ])/[CO 2− 3 ] sat is the degree of undersaturation. Similar to Heinze et al. (2009) we apply a minimum value of 0.062 for the undersaturation thus allowing for some 25 CaCO 3 dissolution in oversaturated waters (as a simplification our model simulates only calcite and not the meta-stable form aragonite or high Mg calcites). Opal dissolution is formulated to depend on the deviation from the saturation concentration for opal in seawater (following Ragueneau et al. (2000); where in the formulation below the effect of changes in free surface area of reactive opal is formally lumped into k opal here): 30 20 Clim. Past Discuss., doi: 10.5194/cp-2016-35, 2016 Manuscript under review for journal Clim. with k opal a standard upper limit opal degradation rate, Si saturation the water column saturation value for silicic acid, and Si actual the actual silicic acid value at a model grid point (the division by 1molel −1 is needed in order to match the units correctly).
The losses in particulate matter according to eq. (A9) for settling particulate material are mirrored by respective source terms for the dissolved species within the water column (for TAlk, DIC, phosphate, oxygen, and silicic acid). The remainder The mass balance equations for a solid sediment component expressed in moles per unit volume of sediment S * and the related concentration of the corresponding dissolved substance C within the pore water are: where D B is the diffusion coefficient for bioturbation, w the vertical advection velocity of solid compound, and G the reaction 20 rate (with S * , C, and G reported here in relation to full sediment volume of a given sediment layer for sake of simplicity; within the model, the varying volumes for dissolved and solid substances following the porosity profile are taken into account).
The prognostic equations for solid sediment concentrations S * (organic carbon, CaCO 3 , opal, and clay) and for dissolved components C (TAlk, DIC, phosphate, oxygen, and silicic acid) are coupled through the respective reaction rates. The diffusive transport of porewater tracer concentrations is simulated via the respective deviations from saturation in parallel with the 25 reduction of this deviation due to chemical porewater reactions: where U is the deviation of the saturation concentration C sat from the actual concentration C) [mole l −1 ], G is the reaction rate (sink for U from dissolution of solid material) [mole l −1 yr −1 ] and D W is the diffusion coefficient. As diffusion coefficient D W we employ a general value of8·10 −6 cm 2 /s for all porewater tracers (this value is in the range of the values given for 30 various pore water species in the work of Li and Gregory (1974). The solid sediment components change due to porewater reactions and particle deposition corresponding to eq. A14: Clim. Past Discuss., doi: 10.5194/cp-2016-35, 2016 Manuscript under review for journal Clim. where S * is the solid sediment component expressed in the same units as U [mole l −1 ], G is the reaction rate (sink due to dissolution) [mole l −1 yr −1 ], and Q is the deposition flux from particle rain [mole l −1 yr −1 ] (the latter occurs only in the top sediment layer). The amount G of solid matter which can be dissolved per unit of time is a function of parameter r * c , the deviation from saturation concentration U , and the amount of solid material available S * : where for r * c we have:

Bioturbation
Bioturbation is simulated as in Heinze et al. (2009) through vertical "diffusion" of the solid substances. This step is carried out after the vertical advection step described further below. The solid matter compounds are slowly mixed in proportion to the 20 prevailing weight fractions in the layers overlying each other: with D B being the bioturbation coefficient. Non-local mixing (e.g., Boudreau and Imboden (1987)) is neglected here. burial > rain − redissolution (accumulation of sediment) and (b) burial < rain − (erosion of sediment), please, see Figure   2 in Heinze et al. (2009). The vertical sediment advection velocity w and also its sign (see eq. A12) depend on the sediment deposition from the water column and the strength of the pore water chemical reactions.

23
Clim . Past Discuss., doi:10.5194/cp-2016-35, 2016 Manuscript under review for journal Clim.    3-D oceanic velocity field: 100% glacial field glacial dust deposition and associated stimulation 7 of biological export production: 100% glacial field 8 Redfield ratio C:P: +15% Table 3. Marine sediment core data and ice core data which were used for determining the parameter changes through a fit of the linear response model. References in parentheses denote the original source of the data in case the data has been taken from a compilation.  threshold value S opal of particle 0.6 export production ratio C(opal):C(POC) for the onset of CaCO3 production (see eq. 5) maximum possible production rain ratio C(