Burial-nutrient feedbacks amplify the sensitivity of atmospheric carbon dioxide to changes in organic matter remineralisation

. Changes in the marine remineralisation of particulate organic matter (POM) and calcium carbonate potentially provide a positive feedback with atmospheric CO 2 and climate change. The responses to changes in remineralisation length scales are systematically mapped with the Bern3D ocean–sediment model for atmospheric CO 2 and tracer ﬁelds for which observations and palaeoproxies exist. Results show that the “sediment burial-nutrient feedback” ampliﬁes the response in atmospheric CO 2 by a factor of four to seven. A transient imbalance between the weathering ﬂux and the burial of organic matter and calcium carbonate lead to sustained changes in the ocean’s phosphate and alkalinity inventory and in turn in surface nutrient availability, marine productivity, and atmospheric CO 2 . It takes decades to centuries to reorganise tracers and ﬂuxes within the ocean, many millennia to approach equilibrium for burial ﬂuxes, while δ 13 C signatures are still changing 200 000 years after the perturbation. At 1.7 ppm m − 1 , atmospheric CO 2 sensitivity is about ﬁfty times larger for a unit change in the remineralisation depth of POM than of calcium carbonate. The results highlight the role of organic matter burial in atmospheric CO 2 and the substantial impacts of seemingly small changes in POM remineralisation.

Changes in the remineralisation of particulate organic matter (POM) have recently been suggested as a further feedback which may contribute considerably to past and future CO 2 changes. How POM is converted to inorganic nutrients in the water column depends on the speed with which the particles sink towards the ocean floor and on the local rate of remineralisation. An increase in temperature is expected to promote bacterial activities and to increase local remineralisation rate (Bendtsen et al., 2002), and to decrease viscosity and thus increase the speed of sinking particles (Taucher et al., 2014). In addition, changes in ecosystem structure in the surface ocean potentially affect the size and density distribution of the settling particles and thus sinking speed (Armstrong et al., 2002(Armstrong et al., , 2009 organic matter cycle influence the surface-to-deep transport of nutrients and carbon and thus respired carbon storage in the deep ocean. A deepening of the remineralisation depth of POM -corresponding to a lower rate of degradation as e.g. induced by lower bacterial activity -results in a drawdown of atmospheric CO 2 . The mechanism of a deepening of the POM remineralisation has been discussed and quantified in several studies.  estimated a CO 2 drawdown of ∼ 35 ppm for a 5 • C ocean cooling experiment and using a threedimensional ocean model. On the other hand, Chikamoto et al. (2012) finds only a small CO 2 sensitivity to a 5 • C ocean cooling (−9 ppm). Kwon et al. (2009) simulate a CO 2 drawdown by 10 and 27 ppm for an increase in the e-folding depth of POM remineralisation of 24 m in their nutrientrestoring and constant export-production model setups. Recently, Segschneider and Bendtsen (2013) estimated a reduction of anthropogenic CO 2 uptake rates of 0.2 gigatons of carbon (Gt C) per year by AD 2100 in future emission scenarios. Taucher et al. (2014) isolate the influence of viscosity changes in a global warming simulation and found that ocean carbon uptake is 17 % higher when considering temperaturedriven viscosity changes on particle sinking speed compared to a control. As a caveat, these studies either focus on the decadal-to century-scale response (Segschneider and Bendtsen, 2013) or neglect ocean-sediment interactions and the weathering-burial cycle.
Ocean-sediment interactions (including weathering and burial) can amplify or mitigate a perturbation in the oceanatmosphere system. For example, excess anthropogenic carbon is ultimately removed from the ocean-atmosphere system by calcium carbonate (CaCO 3 ) sediment buffering and burial. On the other hand, the change in atmospheric CO 2 in response to a change in calcite export is about four times larger in an "open" system considering sediment interactions as compared to a "closed" ocean-atmosphere system (e.g. Tschumi et al., 2011). Concenrning the role of the organic matter cycle, Tschumi et al. (2011) identified the importance of the nutrient-burial feedback for a range of mechanisms with the potential to explain low glacial CO 2 . They conclude that the long-term balance between burial of organic material and tracer input into the ocean through weathering, typically neglected in earlier studies, must be considered when investigating the glacial-interglacial evolution of atmospheric CO 2 and related tracers. However, these authors did not investigate changes in the remineralisation depth of POM.
As the mean remineralisation depth of POM is altered, the deposition of POM on the seafloor and burial rates are also expected to change. This leads to a transient imbalance between burial and continental weathering fluxes on multimillennial timescales. The consequences of these imbalances have not been discussed so far in a systematic manner and in the context of a three-dimensional, dynamic model. Matsumoto et al. (2007) discuss the effect of altered temperature on export and remineralisation rate. They applied an ecosys-tem model embedded in a dynamic ocean model and simulate a change in the export rain ratio between particulate organic carbon and CaCO 3 in response to altered temperatures. They argue that the effect of temperature on POM degradation has a more dominant impact on the export rain ratio than on POM remineralisation depth scale. Accordingly, only the change in rain ratio is prescribed in a box-type oceansediment model to estimate changes in atmospheric CO 2 , but the impact of a deepening of the remineralisation depth as well as changes in the burial of organic matter are neglected by Matsumoto and colleagues. Menviel et al. (2012) performed transient model simulations over the last glacial cycle with the Bern3D model applying a temperature-dependent remineralisation rate. They simulated a 31 ppm decrease in atmospheric CO 2 for a progressive increase of the POM and DOM remineralisation depth over the glaciation, while CO 2 increased by 21 ppm due to altered remineralisation over the deglaciation. This suggests that ocean-sediment interactions may contribute to the reconstructed atmospheric CO 2 variations, motivating the present study.
The goals of this study are to map systematically the spatio-temporal response to alterations in the remineralisation depth in a three-dimensional, dynamic model setting. The response in carbon isotopes, nutrients as well as in atmospheric CO 2 and its carbon isotopic signatures are analysed. A series of well-defined sensitivity experiments is performed, where the POM remineralisation depth is changed in a step-like manner. The Bern3D dynamic ocean model of intermediate complexity is applied with prognostic formulations of export production coupled to a three-dimensional sediment component, allowing for the quantification of the long-term sediment feedback associated with both POM and CaCO 3 burial on marine chemistry and atmospheric CO 2 and 13 CO 2 . In addition, changes in the calcite dissolution rate are considered. We show that previous estimates of the sensitivity of POM remineralisation rate changes are not applicable on glacial-interglacial timescales, as they do not include the long-term feedback, but point to a potential importance of this mechanism to explain low-frequency CO 2 and 13 CO 2 variations. The nutrient-burial feedback (Tschumi et al., 2011) amplifies the initial response to a deepening of the POM remineralisation by a factor of between 3 and 7 in our model.

Physical model
We invoke an updated (see Appendix) version of the Bern3D Earth system model of intermediate complexity (EMIC). Its frictional geostrophic balance three-dimensional ocean component is based on the model of Edwards and Marsh (2005) and is further improved as described in Müller et al. (2006). It includes an isopycnal diffusion scheme and Gent-McWilliams parameterisation for eddy-induced c) d) Figure 1. Preindustrial, annual mean distribution of (a) POM export, (b) calcite export, (c) opal export. The ratio of inorganic to organic carbon export (rain ratio) is shown in (d).
transport (Griffies, 1998). The horizontal resolution is now set to 41 × 40 grid boxes in the horizontal, while 32 logarithmically spaced layers in the vertical are used. Wind stress is prescribed according to the monthly climatology from NCEP/NCAR (Kalnay et al., 1996). The atmosphere is represented by a two-dimensional energy and moisture balance model (EBM) with the same horizontal resolution as the ocean described in Ritz et al. (2011a). Following Weaver et al. (2001), outgoing longwave radiative fluxes are parameterised after Thompson and Warren (1982) with additional radiative forcings due to e.g. CO 2 , other greenhouse gases and volcanic aerosols. The sea-ice model component is based on work by Semtner (1976) and Hibler (1979), and is similar to the sea-ice model of Edwards and Marsh (2005). Air-sea gas exchange for CO 2 and 14 CO 2 is implemented according to the OCMIP-2 protocol Najjar et al., 1999), but with a reduced scaling factor . 14 CO 2 is not fractionated during air-sea gas exchange, while 13 CO 2 is fractionated as a function of temperature and carbonate chemistry, as detailed in Müller et al. (2008).

The marine biogeochemical cycle and particle remineralisation
The marine biogeochemical module computes the cycling of carbon (C), alkalinity (Alk), phosphate (PO 4 ), iron (Fe), oxygen (O 2 ), silica (SiO 2 ), and the carbon isotopes 13 C and 14 C. Carbon is represented as dissolved inorganic carbon (DIC) and labile dissolved organic carbon (DOC).
Prognostic formulations link marine productivity and dissolved organic matter (DOM) to available nutrients (PO 4 , Fe, SiO 2 ), temperature, sea ice and light using Michaelis-Menten limiting terms (Doney et al., 2006) in the euphotic zone (uppermost 75 m), as described in  and Tschumi et al. (2011). POM is exported from the euphotic zone as a constant fraction (σ ) of new production ( new ) at the reference depth z 0 . The remainder (i.e. 1 − σ ) is transferred to the labile DOM pool. The export between CaCO 3 and particulate organic carbon (POC) -i.e. the export "rain ratio" -is constant (∼ 0.075) unless silicic acid is abundant; in this case diatom growth is favoured at the expense of calcifier growth. The rain ratio is therefore lowered in these SiO 2 -replete regions, which is primarily the case in the Southern Ocean upwelling regions (south of ∼ 60 • S; Fig. 1).
The carbon isotope fractionation factor, α [CO 2 ]−C org , between dissolved CO 2 and organic carbon, depends on the concentration of dissolved CO 2 ([CO 2 ]) in units of µmol kg −1 according to Freeman and Hayes (1992):  DOM remineralisation is implemented using a fixed decay time of 0.5 yr. The remineralisation of sinking POM results from the prescribed depth-dependent scaling of the flux of particulate organic phosphorus (POP) and POM (Fig. 2), known as the "Martin curve" (Martin et al., 1987): Here, F POP (z) is the downward flux of POP at depth z (z positive downwards) and F POP (z 0 ) the POP flux at the reference depth z 0 = 75 m. Other elements such as carbon, carbon isotopes, oxygen and alkalinity are coupled to POP and DOP by constant elemental (Redfield) ratios set to P : C : Alk : O 2 = 1 : 117 : −17 : −170. The exponent α defines the shape of the curve: a high value of α shoals while a low α deepens the depth at which organic matter is remineralised. α is set to a globally uniform value of 0.83 in the standard setup.
The empirical relationship of Eq.
(2) -a result of a fit through sediment-trap data -makes per se no statement on the rate of remineralisation, as both sinking speed and degradation contribute to the vertical profile of observed POM concentrations. The remineralisation of POP, R, can be written as k is the remineralisation rate coefficient in s −1 , [POP] the concentration of POP in the water, and v the mean sinking speed of POP. It follows, by rearranging Eq. (3), and with the help of Eq.
(2), for the apparent rate of remineralisation, k app (in units of m −1 ), that The apparent remineralisation rate (k/v) decreases inversely proportional to depth z (Fig. 2b). Although the POM degradation rate is not constant over depth, we define an efolding depth as a length scale (l POM ) for convenience; that is, the depth at which the downward POM flux has diminished to 1/e ≈ 37 % is l POM is 250 m for the standard value of α = 0.83. The flux of organic matter is sketched in Fig. 2. The downward flux of calcite (including other forms of calcium carbonate such as aragonite or high-magnesium calcite particles) (F calc ) decreases exponentially with depth with a length scale l calc : A 10-layer sediment diagenesis model (Heinze et al., 1999;Gehlen et al., 2006) is coupled at the ocean floor. It features the same horizontal resolution as the ocean model. It dynamically calculates the transport, remineralisation/redissolution and bioturbation of solid material within the top 10 cm of the seafloor as well as pore-water chemistry and diffusion as described in detail in Tschumi et al. (2011). Modelled tracers are the four solid components (CaCO 3 , opal, POM and clay) and the eight pore water substances (DIC, DIC-13, DIC-14, total alkalinity, phosphate, nitrate, oxygen and silicic acid). The pore water CO 2− 3 concentration determines whether, and at which rate, CaCO 3 dissolves. The inclusion of the dissolution and burial process of CaCO 3 is crucial for simulating the so-called carbonate compensation.
The oxidation rate of POM within the diagenetic zone depends linearly on the porewater concentration of O 2 and the weight fraction of POM within the solid phase. Denitrification is not taken into account in this version of the model. The corresponding reaction rate parameters are global constants and a decrease in the reactivity of organic material by aging within the diagenetic zone is not considered (Middleburg et al., 1993). Fluxes of carbon and related elements due to POC degradation are coupled by fixed Redfield ratios (P : N : C : O 2 = 1 : 16 : 117 : −170 for oxidation). The model assumes conservation of volume, i.e. the entire column of the sediments is pushed downwards if deposition exceeds redissolution into pore waters. In this manuscript, the term "burial" refers to the net tracer flux at the oceansediment interface, i.e. "burial = deposition − redissolution" of the particulate material. The burial efficiency, i.e. the burial/deposition ratio of a solid species, is controlled by (i) the rate of redissolution within the sediments and (ii) by the rain rate of solid species, which controls how fast the sediment column is pushed downwards.
Any solid material that is pushed out of the diagenetic zone (top 10 cm) disappears into the subjacent diagenetically consolidated zone. The fate of this material is of no further interest for this study (it is known that preferential degradation of POC versus that of POP and the conversion of POP to oxideassociated P and authigenic P within the consolidated zone cause C : P ratios of organic material to deviate substantially from the classical Redfield ratio, Anderson et al., 2001). Input of terrestrial organic matter into the ocean and burial of terrestrial organic matter is not explicitly considered. Similarly, the cycling of P associated with iron and other oxides is neglected, as estimates suggest that 97 % of the P delivered to the sediment-water interface is in the form of organic matter (Delaney, 1998). The specific chemical composition of the organic matter, particle grain size of the sedimentary material and available area for sorption of organic matter (Hedges and Keil, 1995), as well as spatio-temporal variations in mineral deposition rates or sediment porosity, which likely influence organic matter preservation and burial (Burdige, 2007), are neglected.
Here, no land biosphere module is coupled to the Bern3D model, as our interest is in the ocean-atmosphere response to changes in the remineralisation profile of POM and calcite.

Model initialisation and preindustrial state
The model is spun up over 60 000 yr to a preindustrial equilibrium corresponding to AD 1765 boundary conditions. CO 2 is set to 278 ppm and δ 13 C of CO 2 to −6.3 ‰. This results in an oceanic DIC inventory of ∼ 37 400 Gt C, similar to the 37 510 Gt C as estimated from the GLODAP (Key et al., 2004) and World Ocean Atlas Antonov et al., 2010) data sets. The loss of tracers due to burial of particulate matter is compensated during the spin up by a variable weathering so as to conserve oceanic inventories of tracers. After the system equilibrated, the weathering fluxes are kept constant at the rate diagnosed during the end of the spin up. Globally integrated fluxes of POM export, deposition and burial are 11.7, 0.63 and 0.18 Gt C yr −1 for preindustrial conditions.  Tréguer et al. (1995); f Milliman and Droxler (1996).
The partitioning of POM burial between the deep ocean (> 1000 m water depth) and the continental margin (< 1000 m water depth) is 70 and 30 %, respectively, while observations indicate that 80-90 % of the POM burial is on continental margins. This model bias is likely linked to the coarse horizontal resolution and the simple continental runoff scheme, compromising the representation of nearcoast processes. Table 1 provides an overview of steady-state export, deposition, burial fluxes of POM, calcite and opal, and of the weathering-burial fluxes of different tracers. The model is able to represent the observation-based distribution of tracers in the ocean interior (Appendix, Fig. A1). Known shortcomings are a too sluggish formation of Antarctic Bottom Water in the Atlantic, and intermediate waters do not penetrate far enough towards the Equator.
Modelled POM and calcite export ( Fig. 1a and b) is low in the subtropical gyres as well as in the Arctic and around Antarctica. Generally high export fluxes are found in upwelling regions off Africa and South America, in the northern North Pacific and Atlantic, and in the Southern Ocean between ∼ 40 and 60 • S . The pattern is comparable to those 326 R. Roth et al.: Remineralisation simulated by state-of-the-art Earth system models (ESM) and to observation-based estimates (Schneider et al., 2008;Steinacher et al., 2010). Similarly as in ESMs, there is too much new production in the eastern part of the tropical Indian and too little in the Arabian Sea. Modelled opal export ( Fig. 1c) features the well-known belt along 60 • S and maxima in the tropical eastern Pacific and the northwestern Pacific.

Experiments
Sensitivity experiments are performed where l POM (or equivalently α) and/or l calc (Eq. 6) are changed in a step-wise manner (Fig. 3a). The model is then re-equilibrated for up to 200 kyr (some experiments are only integrated for 50 and 100 kyr, respectively) in a prognostic mode where atmospheric CO 2 and δ 13 C evolve freely. As weathering fluxes are kept constant in the perturbation experiments, global inventories of tracers in the atmosphere-ocean-sediment system are not conserved ( Fig. 3h and i). We refer to this setup as "open system". We also run the system in an atmosphere-oceanonly setup without sediment and no river input and burial, referred to as "closed system". The differences in results between the open and closed systems are the contribution by the sediment and weathering-burial feedbacks.
We performed a range of sensitivity experiments where l POM is varied between 225 and 375 m. Our discussion mainly focuses, as an illustrative example, on results from experiments where the exponent α is changed from 0.83 to 0.77, corresponding to a change in l POM from 250 to 275 m.
We restrict our analysis to idealised changes in in the remineralisation length scale (parameter α and l POM ) as the link between remineralisation and temperature changes is not well understood. Changes in ecosystem structure in the euphotic zone and corresponding changes in the quality and sinking speed of POM have the potential to influence remineralisation rates in addition to changes in bacterial activities in the twilight zone and the deep ocean. For illustrative purposes, we may describe a potential relation between temperature changes in the water column and changes in α. Let us assume a spatially uniform temperature change, a constant mean particle settling velocity and that the remineralisation rate of POM is temperature driven and can be de- . With the help of Eq. (2), a change in α from α 0 = 0.83 to α 1 = 0.77 translates, for a Q 10 value of 2, to a mean change in temperature T = 10 • C ln( α 1 α 0 )/ ln(Q 10 ) = −1.1 • C. Changing concentrations of atmospheric CO 2 as induced by biological changes are not fed back to the radiative code in the standard setup, as this so-called oceanic climate-carbon feedback has been systematically analysed elsewhere (Arora et al., 2013;Plattner et al., 2008;Joos et al., 1999). This implies that the physical part of the model (e.g. ocean circulation) remains unchanged during our experiments.

Ocean response
We first analyse changes in fluxes and tracer distributions for an increase in l POM by 25 m from 250 to 275 m. This results in a decrease in the fraction of the POM export that is remineralised above ∼ 250 m and in an increase below (Fig. 2c). These changes affect tracer distributions in the model and in turn surface nutrient availability, productivity, export, deposition, and burial fluxes of POM and calcite (Fig. 3). Both in the open and closed system, the remineralisation of POM at deeper depths leads to a larger surface-to-deep ocean gradient in nutrients, to less nutrients in in the upper thermocline and to less nutrient input into the euphotic zone (note that we only consider changes in remineralisation below the euphotic zone). As a result, surface water concentrations of PO 4 decreases leading to a decrease in global export production of POM and CaCO 3 by ∼ 12 % ( Fig. 3b and c). The reduced export triggers some minor changes in the ecosystem structure in the Southern Ocean, where phosphate depletion favours the growth of diatoms, leading to a local decrease in the export rain ratio. The net effect of the production changes is thus a surface-ocean increase in DIC and Alk, except for the Southern Ocean, where surface Alk remains constant.
The step-wise change in remineralisation depth leads to an initial spike-like increase in POM deposition (Fig. 3d). Then deposition is decreasing in parallel with export, but remains above the initial deposition rate throughout the simulation. The higher deposition leads to a temporary excess of POM burial compared to the closed system (or the weathering input) ( Fig. 3f and h; P and C fluxes scale with the Redfield ratio of 1 : 117) and consequently to a decrease in the ocean inventory of phosphate further reducing organic matter production in the euphotic zone. The decrease in POM and calcite exports occurs in two phases. First, the deepening of the remineralisation depth leads to a redistribution of nutrients within the ocean on a multi-decadal to century scale. Second, the whole ocean phosphate inventory and particle export slowly decrease on millennial timescales due to excess burial of POP to reach a new equilibrium when the balance between burial and weathering input of phosphorus is reestablished around 40 to 50 kyr after the step change (Fig. 3f, dashed line).
Next, we discuss the processes and the time evolution of POM and calcite burial, here taken as the flux leaving the ocean at the sediment interface. Increased POM deposition (rain) tends to increase POM burial and to alter POM oxidation. The change in the amount of POM oxsidised in the sediments varies in space and time: in the deep ocean, reduced oxygen availability in the porewater decreases local remineralisation. On the other hand, increased oxygen levels at coastal margins promote the remineralisation in the sediments in these regions. At, equilibrium POM burial rate has to balance input of phosphorus by weathering. The result is an initial spike in POM burial by almost 80 %, about four times larger than the relative initial increase in deposition, in response to an initial oxygen reduction in the pore water. Afterwards, POM deposition and POM burial decreases quickly within a few centuries and then more slowly to approach steady state over the next millennia. As POM deposition stabilises on a higher level while POM burial relaxes to the initial value, the burial efficiency slightly decreases in the long run by ∼ 5 %. Similarly, the temporal evolution of calcite burial is driven by changes in deposition and redissolution fluxes. Calcite deposition decreases according to the global decrease in primary production and export. As the lysocline shoals as a consequence of the disturbed seawater chemistry -CO 2− 3 initially decreases in the deep ocean -more sediments are exposed to undersaturated waters, leading to an increased redissolution flux for the first ∼ 5 kyr. Thereafter, the lysocline starts to deepen until burial matches weathering input; in the long run, the lysocline deepens by ∼ 800 m in the Pacific. The result is an initial drop in calcite burial by about 10 % and a very slow recovery to the initial burial rate over the next 100 kyr. The imbalance in the burial and weathering fluxes of POM and calcite causes a removal of only 4.2 × 10 14 mol of carbon (5 Gt C) from the ocean and the addition of about 55 × 10 15 equivalent of alkalinity (a mean Alk increase of ∼ 39 µmol L −1 ) ( Fig. 3h and i). The reduced loss of carbon by decreased calcite burial is largely offset by the enhanced burial of POM, while both a reduced calcite burial and an enhanced POM burial tend to increase ocean alkalinity. Overall, the ocean inventory in DIC is also increased as the ocean absorbs about 70 Gt C from the atmosphere (see next section). The whole ocean phosphorus inventory is affected by the excess POP burial, and decreases by 2.2 × 10 14 mol P.
The adjustment to a new equilibrium takes longer for the alkalinity inventory, co-governing calcite burial than for the phosphorus inventory, co-governing POM burial. 63 % of the final perturbation in the phosphorus and alkalinity inventory is reached after 8.5 kyr (Fig. 3f) and after 32 kyr (Fig. 3i), respectively. The difference in timescales is linked to the difference in mean residence time of 23 kyr for phosphorus (3030 T mol/0.13 T mol yr −1 ) and of 253 kyr for alkalinity (3290 P eq/0.013 P eq yr −1 ).
Regional changes in the surface concentration of PO 4 , POM export and POM deposition are distinct (Fig. 4). After 50 kyr PO 4 is reduced almost everywhere in the surface ocean (Fig. 4a). Smallest changes are found in the oligothrophic regions where PO 4 is already strongly depleted in the control simulation. Relatively large reductions in the Southern Ocean and in the upwelling areas off South America and Africa are simulated. These are linked to the upwelling of less PO 4 rich waters as caused by the reduction in the total phosphate inventory of the ocean by excess POP burial. Relatively large reductions are also simulated in the subtropical gyres in the Pacific and Atlantic as well as in the eastern Indian where the model features too high PO 4 concentrations in the standard setup.
POM export (Fig. 4b) is hardly reduced in the Southern Ocean, while substantial reductions are simulated in equatorial regions and in the northern North Pacific and North Atlantic. POM deposition (Fig. 4c) increases along the productivity belt in the Southern Ocean and in the upwelling regions off South America and Africa and in the northwestern North Pacific. There is also an increase in POM deposition in the eastern tropical Indian Ocean, which is linked to too excessive productivity in this area. In the North Atlantic, the influences of the reduction in POM export and the increase in remineralisation depth on POM deposition nearly cancel.
Next, we discuss the redistribution of tracers within the ocean for the new equilibrium in the closed and open system. The distribution changes in the closed system occur in the first two millennia, while it takes many millennia to approach the equilibrium in the open system. Figure 5 shows the resulting equilibrium anomalies w.r.t. the control simulation for PO 4 , oxygen, DIC, Alk and CO 2− 3 in a transect through the Atlantic, Southern Ocean and Pacific (the transect is highlighted in red in Fig. A2b).
In the closed system (i.e. without the sediment module), the ocean's PO 4 inventory is conserved, and a deepening of the remineralisation depth leads to a depletion in PO 4 in the upper ocean and an increase in the deep (Fig. 5a). The pattern reflects the major transport pathways in the ocean showing reduced PO 4 concentration in North Atlantic Deep Water and enhanced concentration in Antarctic Bottom Water in the Atlantic and Pacific. As expected, this pattern is similar in subsurface waters for DIC (Fig. 5e) and inverse for δ 13 C (Fig. 7d) and oxygen (Fig. 5c) as these tracers are linked to PO 4 by constant Redfield elemental ratios in biological fluxes. At the same time, CO 2− 3 and Alk ( Fig. 5g and i) are slightly increased in the upper part of the water column as less POM is remineralised in the upper ocean (POM remineralisation consumes Alk). CO 2− 3 (and Alk) is lowered in the deep as less calcite and more POM is remineralised.
In the setup with the sediment model, downward fluxes of POM (and CaCO 3 ) are not remineralised at the seafloor, but are deposited onto the model's sediment layer, and undergo bioturbation, remineralisation and vertical transport. The ocean inventory of PO 4 , and Alk change substantially in response to excess POM and reduced calcite burial as discussed above. As a result, PO 4 is reduced and Alk increased in the entire ocean compared to the control simulation. The changes in DIC are more subtle and regionally distinct. Reduced remineralisation in the thermocline leads to negative anomalies in the Atlantic and Pacific. The increase in Alk causes a change in the partitioning between dissolved CO 2 , bicarbonate, and carbonate ion. The result is a decrease in pCO 2 , an increase in carbonate ion, and an increase in DIC in the surface ocean. These positive DIC anomalies are communicated to the deep in the North Atlantic and Southern     whole-ocean increase in CO 2− 3 leads to a decrease in seawater pCO 2 at the surface.
How does this influence the atmospheric CO 2 concentrations? In Fig. 6a, time series of the initial 50 kyr response are shown for both the closed (black) and open system (red). The magnitude of changes in atmospheric CO 2 is ∼ 5 times higher in the setup with sediments compared to the setup without sediments at year 50 000. Atmospheric CO 2 decreased by 33.4 ppm in the open system and by 8.1 ppm in the closed system for a change in l POM from 250 to 275 m. The CO 2 decrease is more than 100 ppm for a change in l POM to 375 m. Atmospheric CO 2 is stabilised after ∼ 2 kyr in the closed system, while at the same time the response in the open system is only ∼ 40 % of the equilibrium response. 62 % of the equilibrium CO 2 change are reached after ∼ 10 kyr and 93 % after ∼ 50 kyr. The effective e-folding timescale, defined as the time when 63 % of the final anomaly is reached, is ∼ 10 kyr in the open system and only 0.5 kyr in the closed system.
Global export production decreases when increasing the remineralisation length scale in the standard model setup. We run a sensitivity experiment to explore the influence of the reduced export flux vs. the change in the remineralisation profile. The export production climatology from the control run is prescribed in these perturbation experiments. As to be expected, changes in CO 2 increase by a factor of ∼ 2-4 w.r.t. to the standard setup with prognostic production (see Fig. 6b). The constant export flux leads to a constant drain of nutrients and carbon from the ocean as burial fluxes continuously exceed weathering fluxes in the open system. In turn, the oceanic inventories slowly drift towards zero and negative values. Therefore, no equilibrium changes in CO 2 can be stated for the prescribed production model (Fig. 6b shows anomalies after 50 kyr). Next, we assess the sensitivity of atmospheric CO 2 to l POM changes in different oceanic regions ( Table 2). The hypothesised change in the remineralisation depth may be different in different regions. l POM is changed in 3 latitudinal bands separately, namely 90-30 • N, 30 • N-30 • S and 30-90 • S. The sum of the CO 2 changes for the three regional experiments is with 35.6 ppm only slightly larger than for the global experiment (33.4 ppm). The tropical band covers the largest sea surface area (52 %) and most of the global POM export flux (∼ 55 %). For a change of l POM from 250 to 275 m, the tropical band alone contributes ∼ 60 % of the CO 2 change. The effectiveness of l POM changes per unit area or per unit POM export to reduce CO 2 is different for different regions (Table 2). Normalised sensitivities are highest in the northern extratropical band, likely linked to an efficient propagation of anomalies by North Atlantic Deep Water.

Carbon isotopes
Next, we further discuss the response in δ 13 C of DIC and of atmospheric CO 2 . δ 13 C of atmospheric CO 2 is closely coupled through air-sea exchange to δ 13 C of DIC in the surface ocean. We start with the closed system response. The spatiotemporal response in δ 13 C of DIC mirrors that of PO 4 as the marine biogeochemical cycles of P and C and 13 C are coupled. The deepening of l POM leads to less remineralisation of isotopically light (∼ −20 ‰) POM in the upper ocean and more remineralisation in the deep. In turn δ 13 C of DIC in the deep ocean decreases and δ 13 C at depths shallower than ∼ 1000 m increases (Fig. 7d). In addition, the decrease in surface water concentration of CO 2 leads to slightly less negative values of the fractionation factor for photosynthesis and thus for the transformation of inorganic carbon to POC; this tends to slightly decrease δ 13 C in the surface ocean. The result is an increase in surface ocean and atmospheric δ 13 CO 2 by ∼ 0.03 ‰ within 2 kyr ( Fig. 7a and d). Additional processes are important for the evolution of δ 13 C of DIC and CO 2 in the open system. (i) The excess burial of isotopically light POC (δ 13 C ≈ −20 ‰) during the first ∼ 50 kyr tends to increase the average δ 13 C signature in the ocean-atmosphere system; the δ 13 C signature of the POC burial flux is, at −20 ‰, about 7 ‰ lower than the signature of the total carbon weathering-burial flux (−12.6 ‰). (ii) Similarly, burial of isotopically enriched calcite (δ 13 C ∼ 3 ‰) is reduced (relative to initial conditions and the weathering flux) during the first ∼ 100 kyr; this tends to increase δ 13 C during this period; the calcite burial flux is enriched by about 15 ‰ compared to the average signature of the weathering-burial flux. (iii) The reduction in surface CO 2 is, at 33 ppm, about 4 times larger in the open system than in the closed system for an increase in l POM to 275 m. This results in a less negative fractionation factor for POM formation (Freeman and Hayes, 1992) and tends to decrease δ 13 C in the surface. (iv) At the new equilibrium, the input of 13 C by weathering must balance the loss of 13 C by burial, as is the case for the flux of POC and calcite (Fig. 8). The consequence is that the positive shift in the POC fractionation factor leads to a loss of 13 C from the ocean until it is balanced by a corresponding negative anomaly in 13 C of DIC in surface waters communicated to the POC and calcite burial fluxes.
In the surface ocean and the atmosphere, δ 13 C increases during roughly the first 10-15 kyr by 0.08 ‰ in response to the closed system processes and in response to the anomalies in POC and calcite burial. Then, δ 13 C decreases and the anomalies in δ 13 C turn negative around 50 kyr followed by a further decrease until the end of the simulation in response to (iii) and (iv). The fractionation factor (expressed here as = 1000 · (α − 1)) for photosynthesis (POC formation) changes roughly in parallel with surface ocean and atmospheric CO 2 . It increases by 0.8 ‰ from −13.2 to −12.4 ‰ in the global average during the first ∼ 100 kyr in the simulation where l POM equals 275 m (Fig. 8a). In turn, the δ 13 C signature of the POM burial flux increases by about 0.55 ‰ and starts to decreases only after about 50 kyr (Fig. 8b). δ 13 C of the burial flux of calcite closely follows the evolution of δ 13 C in the surface-atmosphere system. The whole ocean's 13 C inventory is driven by changes in the δ 13 C burial flux, expressed as the burial flux of C times its isotopic signature in ‰. The 13 C burial ( Fig. 8c and d) is thus a convolution of the total C burial, the particle's δ 13 C and the POCto-calcite ratio in the burial flux. The new equilibrium with zero δ 13 C anomaly in the total burial flux is only achieved after the end of the simulation (200 kyr) (Fig. 8d). At the end of the simulation, the whole-ocean 13 C of DIC (and DOC) is still relaxing from its initial perturbation and a stabilisation is expected to take another 200 kyr or so (Fig. 8e). These long response timescales are a result of the convolution of multiple slow processes and weak negative (stabilising) feedbacks.
The effect on the radioactive carbon isotope 14 C is small compared to past natural variations. 14 C of atmospheric CO 2 increases by ∼ 12 ‰ (∼ 3 ‰ in the closed system). A constant rate of atmospheric 14 C production is applied. The increase can be attributed to the lower atmospheric CO 2 concentration, or more generally, to the lower carbon inventory in the (open) system. The biologically mediated fluxes of 14 C  within the ocean are small compared to the fluxes driven by advection, convection, and mixing. This explains the small sensitivity of 14 C to changes in the marine biological cycle.

l calc changes
It has been suggested that changes in the rate of CaCO 3 dissolution in the upper ocean will be a significant feedback affecting atmospheric CO 2 concentrations and future climatic changes (e.g. Barrett et al., 2014). To this end, we prescribe in a further set of sensitivity simulations a step change in the e-folding dissolution length scale l calc governing the dissolution profile of CaCO 3 particles within the water column (see Eq. 6), l calc is changed at the end of the spin up from its standard value of 2900 m to values ranging from 2100 to 3700 m; then the run is continued for another 50,000 years with the new value of l calc . Export fluxes both of POC and CaCO 3 remain constant as changes in the dissolution of CaCO 3 do not affect productivity in our model.
The experimental setup with the assumption of an efolding remineralisation profile for CaCO 3 particles (Eq. 6) is highly idealised. The mechanisms for the dissolution of calcite and other forms of CaCO 3 (e.g. aragonite or highmagnesium calcite) within the water column are quantitatively not well understood. Generally, dissolution of CaCO 3 particles within the water column is thought to be linked to low (undersaturated) concentration of carbonate ions in the surrounding water. However, considerable CaCO 3 dissolution may occur in the upper ocean (Berelson et al., 2007) in waters that are saturated with respect to CaCO 3 in the mineral form of calcite or even of aragonite, perhaps due to the dissolution of high-magnesium calcite or due to dissolution in specific microenvironments such as zooplankton guts, fecal pellets or organic aggregates (see e.g. Barrett et al., 2014, and references therein). Here, we do not explicitly take into account such effects, and the dissolution is for example not affected by changes in the calcite saturation horizon (the depth where the concentration of calcium times the concentration of carbonate equals the solubility product of calcite; Mucci, 1983).
The results for these additional sensitivity experiments are as follows. A deepening of l calc (i.e. a lower calcite redissolution rate) leads to an increase of dissolution products -DIC and Alk -in the deep ocean (below ∼ 3 km) and a decrease above. This tends to hinder redissolution of calcite from sediment in the deepest layers. Calcite deposition on sediments and burial both increase in response to the increased dissolution length scale. As a consequence, Alk and DIC burial and thus the removal from the ocean increase with a ratio of 2 : 1. In the long run, the balance between weathering input and burial is reestablished by the following process. The removal of Alk and DIC causes an upward shift in the lysocline. In turn, a larger fraction of the calcite deposition flux is redissolved, and less is buried.
The ocean inventories of DIC and Alk are both perturbed towards lower values when l calc is increased. This is different to a deepening of the remineralisation depth for POM, where the Alk and DIC budgets are perturbed in opposite directions. As a result, atmospheric CO 2 is less sensitive to changes in l calc than to changes in l POM . For example, a global deepening of l calc from 2900 to 3300 m results in an increase in CO 2 by 12.2 ppm in the open system (and only 1.8 ppm in the closed system). Although the absolute sensitivity is lower than for l POM , the amplification by weathering-burial dynamics is even more pronounced, with an ∼ seven-fold higher response in the open than in the closed system. As oceanic δ 13 C of DIC is primarily controlled by the organic matter cycle, changes in l calc only marginally influence 13 CO 2 , such that the signal is too low to separate it from model-internal variability and drift. Table 3. Equilibrium changes in atmospheric CO 2 and 13 CO 2 (expressed as δ 13 C) induced by a global change in the mean remineralisation depth of POM and the calcite redissolution length scale.

Sensitivities and sediment feedback quantification
In the following, sensitivities to changes in remineralisation depth are calculated from an experiment where l POM is reduced by 25 m and an experiment where it is increased by 25 m (Table 3). The difference in CO 2 or δ 13 C between the two experiments is divided by 50 m, the total change in the e-folding depth. This corresponds to a sensitivity calculated symmetrically around the standard (best-guess) value. In general, atmospheric CO 2 responds more sensitively to a shoaling of l POM and l calc than to a deepening. The equilibrium sensitivity of atmospheric CO 2 to changes in l POM is ∼ 1.68 ppm m −1 in the open system, and about 4 times lower, 0.38 ppm m −1 , in the closed system. Sensitivities in 13 CO 2 (Table 3) are less amplified by the weatheringburial loop due to pCO 2 -induced changes in the fractionation factor. Next, combined changes in l POM and l calc are examined. We sampled the two-dimensional parameter space and integrated each model for 20 kyr. To get equilibrium changes in atmospheric CO 2 , a linear combination of four exponential functions with different timescales is fitted to the response of the first 20 kyr. The equilibrium response is then found analytically for t → ∞. The result of this exercise is shown in Fig. 9a and b. As indicated by the green line in Fig. 9b, LGM CO 2 levels could e.g. be reached by changing l POM to ∼ 345 m. This number is reduced if l calc is decreased simultaneously. In the closed system on the other hand, LGM CO 2 levels are far from being reconcilable within the bounds of our parameter sampling. The amplification in the response of CO 2 is calculated as the ratio CO 2,open / CO 2,closed and is shown in Fig. 9c. The sediment amplification range is 3-7 for changes in l POM and 7-8 for changes in l calc .

Discussion and conclusion
We systematically mapped the spatio-temporal responses for a variety of tracers and proxies to changes in the remineralisation depth of organic matter and calcium carbonate. The Bern3D dynamic ocean model including sediment interactions with organic matter, opal and calcium carbonate is applied. It is shown that on long timescales, ocean-sediment interaction and the weathering-burial cycle strongly amplify the ocean-only response in atmospheric CO 2 , δ 13 C, ocean tracers and particle export fluxes. These processes also lead to sustained changes in whole ocean nutrient and alkalinity inventories. The adjustment to the new steady state occurs on a number of timescales. Ocean-only reorganisations take place within decades to centuries, and most of the changes in particle export occur within less than 2 kyr, while the adjustment of the phosphate and alkalinity inventories occurs on a typical timescale of about 10-30 kyr. Even longer timescales are involved in the evolution of δ 13 C. δ 13 C signatures are not yet in equilibrium after 200 kyr, the end of our simulations. This implies, as already highlighted by Tschumi et al. (2011), albeit for other mechanisms, that ocean-sediment interactions and the burial-weathering cycles must be considered when discussing the long-term evolution of atmospheric CO 2 and δ 13 C.
Our model study is associated with a number of limitations. The experiments are idealised and not intended to track the evolution of palaeoproxy signals directly. Changes in remineralisation rates of particulate organic matter (POM) and calcium carbonate (CaCO 3 ) are treated in a parameterised, globally uniform way using a power law and an exponential scaling for the attenuation of POM and CaCO 3 fluxes with depth. It is a task for future research to advance the mechanistic understanding of the processes governing particle fluxes and remineralisation rates. Furthermore, the Bern3D model is a coarse resolution model with simple parameterisations for the productivity of organic matter, opal, and CaCO 3 , and fixed Redfield ratios are applied. Productivity is linked to phosphate, as well as to iron and silicic acid, but not to nitrogen compounds. Thus, responses in the nitrogen cycle are not taken into account and these may alter in particular the response in export production. Implicitly, the setup corresponds to the assumption that nitrogen fixers make up for the loss of nitrogen by POM burial to the extent needed to support production. Our model does not include a comprehensive ecosystem model, and may underestimate the shift in ecosystem structures due to the nutrient drawdown as induced by the slowed POM remineralisation. For example, Segschneider and Bendtsen (2013) find an opposite effect of nutrient (i.e. nitrate) depletion on surface Alk using a more complex biogeochemistry model.
There are also limitations regarding the sediment model. For example, the spatio-temporal variability in the deposition of mineral particles or the influence of particle grain size on organic matter preservation is neglected. The coarse resolution hampers the representation of coastal and continental boundaries, where most POM deposition, remineralisation and burial occurs (e.g. Wallmann et al., 2012). The model does not resolve river deltas and estuaries and their carbon cycle (e.g. Regnier et al., 2013). Another caveat is that denitrification within the sediment is not represented by our model, eventually leading to a bias in the long-term response of POM degradation and thus burial efficiencies. Therefore, our findings are to be confirmed and refined by a higherresolved ocean model with a more complete representation of sediment processes.
The starting point for all sensitivity experiments is a steady state corresponding to preindustrial conditions, and sensitivity of results may vary somewhat between different background states. A further caveat in our model could be that weathering rates are kept constant. On timescales of several 10 kyr riverine input of tracers are thought to adapt to changes in climatic conditions as well as in pCO 2 (e.g. Colbourn et al., 2013, and references therein). As this negative feedback is not included in our model, the long-term response may be overestimated. Slowed bacterial degradation of POM may also imply a change in the mean lifetime of DOM, which we kept constant in our experiment. The impact of a corresponding increase in DOM lifetime is found to weaken the overall response in CO 2 to a change in remineralisation depth of POM by about 15 % (Menviel et al., 2012). Finally, we did not take into account the radiative feedback associated with altered CO 2 to ease interpretation of results. The climate-carbon feedback is positive in our model, and its inclusion would yield an even larger change in atmospheric CO 2 . Despite these numerous caveats, we expect that our experiments will provide the first-order response to changes in the remineralisation depth of POM and CaCO 3 .
The results show a high sensitivity of atmospheric CO 2 to changes in the POM remineralisation depth. A new element, compared to earlier work with dynamic three-dimensional ocean models, is that changes in the POM remineralisation rate not only change the atmosphere-ocean partitioning of carbon, but also impact the long-term carbon cycle by disturbing the weathering-burial balance of carbon, alkalinity and nutrients. This temporary imbalance amplifies the ocean-atmosphere-only response by a factor of ∼ 3-7 on multi-millennial timescales. A subtle shift of 1 m (or 4 ‰) in the length scale of POM remineralisation causes a change in atmospheric CO 2 of almost 2 ppm. Our results are in quantitative agreement with Kwon et al. (2009) and when sediment interactions are neglected.  find a sensitivity in atmospheric CO 2 corresponding to about 7-8 ppm per degree warming and assuming a Q 10 factor for remineralisation of 2. Kwon et al. (2009) find -for l POM deepening of 24 m -a CO 2 reduction of 10 ppm using a nutrient restoring formulation of export production (and 27 ppm with constant export) with a model of comparable complexity, but without taking into account ocean-sediment interactions. In our model featuring prognostic export, we find a reduction of 8 ppm in our closed system setup, i.e. without sediment interactions, for l POM deepening of 25 m. In the open system, weatheringburial dynamics increase the response four-fold to 33 ppm. It has been suggested that changes in the rate of CaCO 3 dissolution in the upper ocean may be a significant feedback affecting atmospheric CO 2 concentrations and future climatic changes (Barrett et al., 2014). In our sensitivity experiments, we varied the e-folding dissolution length scale between 2300 and 3500 m, and thereby the amount of CaCO 3 dissolution in near-surface water and the upper ocean. These variations have a relatively small influence on atmospheric CO 2 in our model on both century (Fig. 9a) and glacialinterglacial timescales (Fig. 9b). This is consistent with the results of Gangstø et al. (2011), who applied the PISCES ecosystem model with CaCO 3 production in the form of calcite and aragonite over the industrial period and in future scenarios. However, these authors restricted the dissolution of CaCO 3 to waters that are undersaturated with respect to calcite and aragonite and thus to deeper layers (in particular for calcite), and did not apply a sediment model. Taken together, changes in upper ocean and deep ocean remineralisation rates of CaCO 3 appear to exert a small influence on atmospheric CO 2 and thus climate.
The idea of the glacial drawdown of nutrients and carbon from the surface to the deep ocean -either by a more stratified ocean and more sluggish ventilation of the deep ocean or other alterations of the marine biological cycle -and as a consequence, the build-up of old respired carbon, has been widely discussed (Boyle, 1988). Evidence for increased glacial storage of respired carbon comes from oxygen-sensitive sediment proxies which show that the oxygenation was decreased during glacials. As deep glacial Cd/Ca ratio in the North Pacific, a proxy for PO 4 , was not increased, these findings lead to the "respired-carbon deepening hypothesis" which denote an increase in remineralised : preformed nutrient supply to the deep ocean (Jaccard et al., 2009;Bradtmiller et al., 2010;Jaccard and Galbraith, 2011). Figure 10 shows the ratio of remineralised : preformed nutrients (PO 4,rem /PO 4,pre ) in our model -based on apparent oxygen utilisation -as well as the change in this ratio induced by an increase in l POM . A 25 m increase in the e-folding depth of POM remineralisation leads indeed to a widespread increase in PO 4,rem /PO 4,pre by ∼ 3 % in the global average. The response in CO 2 to a deepening of the remineralisation depth evolves on timescales that are much longer than the typical duration of the transition from a glacial to an interglacial state, e.g. ∼ 7 kyr for termination I. As noted by Menviel et al. (2012), this mechanism can thus not explain, at least not in isolation, the relatively rapid increase in CO 2 over glacial terminations. In agreement, the whole ocean increase in Alk and CO 2− 3 resulting over long timescales from a deeper remineralisation depth is not observed in proxy records, which show little glacialinterglacial changes in CO 2− 3 (Yu et al., 2010(Yu et al., , 2013. The results confirm that the nutrient-burial feedback identified by Tschumi et al. (2011) should not be neglected (e.g. Sigman et al., 1998;Matsumoto et al., 2007;Chikamoto et al., 2008) when discussing long-term carbon cycle changes. This nutrient-burial feedback arising as marine biological productivity in the long run is limited by the   riverine or aeolian input of nutrients into the ocean. The burial flux, not only of CaCO 3 but also of POM, must relax to the input flux of tracers into the ocean. This nutrientburial feedback amplifies the initial response to a deepening of the POM remineralisation by a factor of between 3 and 7 in our model. There are small changes in 14 C in the atmosphere and the ocean in our POM-deepening experiments.
On the other hand, we find substantial changes in alkalinity, DIC, oxygen, and atmospheric CO 2 . This suggests that the relationship between 14 C and DIC observed in the modern ocean (Sarnthein et al., 2013) may not always be applicable to infer changes in DIC from changes in 14 C. We have shown that a change in the rate of organic matter degradation constitutes an important potential feedback for atmospheric CO 2 to temperature change. Small changes in the remineralisation depth lead to surprisingly large changes in atmospheric CO 2 and related tracers in particular over glacial-interglacial timescales. While the available proxy records do not support the conclusion that a deepening in the POM remineralisation depth played a dominant role in the CO 2 increase during the last termination, the high sensitivity implies that subtle changes in the remineralisation depth on global to regional scales cannot be easily disregarded for the explanation of the long-term evolution of atmospheric CO 2 and δ 13 C, the distribution of nutrients, alkalinity, oxygen, carbonate ions or δ 13 C within the ocean, as well as for marine productivity, particle fluxes, and sedimentation. Recently, the physical core of the Bern3D ocean component underwent some major changes, leading to notable improvements (Fig. A1). Before this update, only horizontal grids with a latitudinal spacing proportional to the sine of the latitude were possible, such that the surface area of all gridboxes is equal. Based on this constraint, many equations were simplified in the original Bern3D model (Müller et al., 2006) (which in turn is based on the GOLDSTEIN model; Edwards and Marsh, 2005), as gridcell areas were omitted in the equations.
All equations were generalised to be used with arbitrary rectangular grids. Several horizontal grid configurations were tested in an early stage. A grid with 41 × 40 gridboxes in the horizontal featuring a higher latitudinal resolution in the Southern Ocean as well as an increased longitudinal resolution in the Atlantic has been chosen. The old and new grids are shown in Fig. A2. This particular 41 × 40 grid offers a good balance between computational efficiency and increased resolution in important locations. Although the number of boxes only modestly increased, the timestep (δt) had to be decreased by a factor of two in order to maintain numerical stability.
In addition to the horizontal resolution, some minor changes were incorporated into the ocean, EBM and biogeochemical modules, as described in the following. As a result of a proper re-tuning of all model components, there is some notable improvement in the modelled ocean distribution of tracers, as illustrated in the Taylor diagram shown in Fig. A1. For a more detailed description and the present-day solution, see Roth (2013, Chap. 2).

A1 Physical component
The following changes in the physical part were made (parameters w.r.t. the 36 × 36 setup as described in Ritz et al., 2011a andupdated in Ritz et al., 2011b): -new grid with 41 × 40 gridboxes in the horizontal. -Ocean: -new parameter set (see Table A1) -max. isopycnal slope has been increased to prevent spatial oscillations of the tracer fields in the Southern Ocean -shuffling convection is now applied for all tracers (instead of only to T and S). -Atmosphere: -new parameter set (see Table A1) -no Atlantic-to-Pacific freshwater correction flux -zonal resolved winds are now used for moisture and heat advection -zonal eddy-diffusive transport of heat is now a function of latitude (see Table A1) -moisture advection is now reduced by a scaling factor β moist = 0.5 -land albedo changes through snow are now taken into account by default. Changes in snow albedo calculated by the snow albedo parametrisation of Ritz et al. (2011b) are accounted for as anomalies to the fixed surface albedo. Minimum snow albedo is the MODIS 1 land surface albedo from which the snow contribution has been removed -air-sea momentum transfer (i.e. wind stress) is reduced by fractional sea-ice cover (F ice ) by a factor of (1 − F 2 ice ).
-Sea ice: -sea ice is solved with a separate time-stepping: 10 steps per atmospheric time step -sea ice advection is reduced to 0.3 times the surface ocean current for better consistency with observations.
-Ice sheet: -new ice-sheet parametrisation: Greenland and Antarctica are now assumed to have temporally constant ice-sheet cover. Due to this, the climate is no longer affected when ice sheets are activated at a present day state (albedo of Greenland and Antarctica remain at the values of the MODIS climatology) -updated ice-sheet mask (Peltier, 2004).
In the 41 × 40 version, the key ocean parameters remained largely unchanged, with the exception of an ad hoc doubling of diapycnal diffusivity (k D ) from 10 −5 m 2 s −1 to 2 × 10 −5 m 2 s −1 . This change was motivated by the lowerthan-observed rates of anthropogenic CO 2 uptake in the old model version (Gerber and Joos, 2013). With the new model, the cumulative uptake of CO 2 from AD 1765 to 2011 is 156 Gt C, which compares well with observational estimates of 155 ± 31 Gt C (Le Quéré et al., 2013).

A2 Biogeochemical component
The following changes in the biogeochemical component were made (parameters w.r.t. the 36 × 36 setup as described in , and Tschumi et al. 2008: -Ocean biogeochemistry: -new parameter set (see Table A2) -new windspeed climatology (from NCEP instead of PO.DACC) for the gas exchange formulation -new present-day dust input fields from Mahowald et al. (2006) -virtual fluxes are now applied to all biogeochemical tracers. -Sediment: -new parameter set (see Table A2) -Alk : P ratio for organic matter formation and remineralisation is set to −17 according to Paulmier et al. (2009) -weathering input is now distributed along the coastline instead of the entire surface ocean -organic matter denitrification disabled per default to avoid model drift.  Figure A2. Comparison of the horizontal grid for the previous (36 × 36) (a) and updated (41 × 40) (b) versions of the Bern3D ocean model. Both versions feature 32 depth layers. Latitudinal resolution near the Equator remained largely unchanged, while resolution was increased at the high latitudes. The ocean, atmosphere and sediment models share the same horizontal grid. The red line in (b) shows the path of the transect plot used in the main text. The light colours depict the definition of ocean basins as used to calculate depth gradients and zonal means.