GEOCLIM reloaded ( v 1 . 0 ) : a new coupled earth system model for past climate change

Introduction Conclusions References


Introduction
In the sedimentary record, organic matter-rich strata (e.g. black shales, sapropels), significant isotopic excursions (e.g. rapid carbon isotope excursions) or unusual mineral deposits (e.g. capstone carbonates, banded-iron formations) bear the testimony of past extreme climate events that punctuated the evolution of the Earth climate system (e.g. the Correspondence to: S. Arndt (arndt@geo.uu.nl) Neoproterozoic glaciations, Jurassic and Cretaceous oceanic anoxic events, or the Palaeocene/Eocene Thermal Maximum). Over the past decades, concern about ongoing climate change has increased the interest in these events. Their causes are manifold and include among others plate tectonics, cyclic variations in orbital parameters of the Earth, varying intensities of volcanic activity, or meteorite impacts (e.g. Herget, 2006). Despite the different nature of these triggers, all extreme events in Earth history can be related to important variations in greenhouse gas concentrations, including atmospheric carbon dioxide (CO 2 ). This correlation finds its origin in the prominent role of atmospheric CO 2 concentrations in regulating the global climate. Over the past decades, considerable progress has been made in understanding and quantifying the feedbacks between the global carbon cycle and the climate system over geological timescales (e.g. Berner and Kothavala, 2001;Mackenzie et al., 2004;Donnadieu et al., 2004;Ridgwell and Zeebe, 2005;Donnadieu et al., 2006;Goddéris et al., 2008b) and references therein. It has been shown that the long-term evolution of the climate system is driven by a dynamic balance between CO 2 inputs from volcanoes, mid-oceanic ridges and metamorphic alteration of rock, and CO 2 sinks due to weathering and burial in marine sediments. On long timescales, the carbon cycle generally exerts a stabilizing mechanism on the climate system through a feed back mechanisms first described by Högbom (1894) (e.g. Berner, 1995 and now often referred to as the Walker thermostat (Walker et al., 1981). This planetary thermostat operates by adjusting the global CO 2 removal through chemical weathering and subsequent burial to the supply of CO 2 degassing from the deep, solid Earth. Such a self-regulating effect is a characteristic feature of complex systems and results from a high number of non-linear 452 S. Arndt et al.: GEOCLIM reloaded (v 1.0): a new coupled earth system interactions and feedbacks which, together, dampen the perturbations imposed on the system. Yet, important imbalances between these regulatory fluxes can arise on short timescales (10 5 − 10 6 yr) and complicate the response of this feedback loop to perturbations. The increasing temporal resolution of observational data, for instance, has revealed the intermittent decoupling between weathering rates and carbonate burial when tectonic processes lead to subsidence of platforms or disappearance of reefs (Compton and Mallinson, 1996;Opdyke and Walker, 1992;Ridgwell et al., 2003). The coupling between chemical weathering and organic carbon burial is even more prone to flux imbalances, since it is less direct (e.g. Berner, 2002). As a consequence, the short-term response of the coupled carbon-climate system to rapid perturbations is still poorly understood. Furthermore, the very different timescales that characterize elemental cycling through the reservoirs of the Earth system (atmosphere 10 1 yr, continents 10 3 yr, ocean 10 3 − 10 4 yr, sediments 10 6 − 10 8 yr) and the inability of observations to fully resolve all of the characteristic spatial and temporal scales of change challenge simple, straightforward answers.
Earth system models provide a powerful tool to resolve the coupling of (geo)physical, (geo)chemical and biological processes within the different compartments of the complex Earth system, at the appropriate temporal scales. When combined with the rapidly growing set of proxy and biomarker data this approach can also provide important insights into the functioning of the carbon-climate system. In particular Earth system models have improved our understanding and quantification of the effect of varying CO 2 levels on the longterm evolution of global biogeochemical cycles and climate (e.g. Goddéris and François, 1995;Berner and Kothavala, 2001;Wallmann, 2001;Mackenzie et al., 2004;Ridgwell and Zeebe, 2005). Yet, owing to the high computational needs required to perform numerical simulations over geological timescales, these models are generally characterized by a high degree of simplification. They often focus on specific aspects of the coupled carbon-climate system and rest on many simplified global relationships to constrain exchange and transformation processes. While these models have proven extremely useful to investigate the longterm climatic evolution, they are often not suitable to resolve the short-term response to perturbations. A major caveat of existing models is their inability to fully resolve diagenetic processes in marine sediment. Diagenetic processes link the rapid oceanic, atmospheric and terrestrial cycles to the slow geological rock cycle and are thus key components in the global carbon cycle (Arthur et al., 1988;Berner, 1990;Archer and Maier-Reimer, 1994;Mackenzie et al., 2004;Ridgwell and Zeebe, 2005;Ridgwell et al., 2007). They trigger a delayed response to changes in ocean geochemistry and control the temporal and permanent removal of carbon from the ocean reservoir. Yet, existing Earth system models do not resolve appropriately the benthic-pelagic coupling and carbon burial is thus generally described by simplified para-metric laws based on empirical relationships from modern day observations. However, these parametric laws do not account for the complex interplay between the different processes and forcing mechanisms that determine the diagenetic response and may not be applicable under the environmental conditions that prevailed during past extreme events. Here, we present a new version of the coupled Earth system model GEOCLIM (Donnadieu et al., 2004(Donnadieu et al., , 2006Goddéris et al., 2008b), the first numerical model that used a general circulation model to simulate the global carbon and alkalinity cycles in the geological past. Over the past years, GEO-CLIM has been successfully applied to explore the feedbacks between continental weathering and the evolution of the Mesozoic, Paleozoic as well as Neoproterozoic climate (Donnadieu et al., 2004(Donnadieu et al., , 2006Goddéris et al., 2008a;Le Hir et al., 2008a,b), late Devonian climate change (Goddéris and Joachimski, 2004), geochemical cycling during a 'snowball' glaciation (Le Hir et al., 2008c), the feedback between continental vegetation and climate (Donnadieu et al., 2009), as well as the influence of climatic conditions and the evolution of calcifying organisms . Although GEOCLIM is a powerful tool for reconstructing the long term evolution of the Earth climate, it only includes a low resolution box model for the global ocean and lacks an explicit description of diagenetic processes. Therefore, it is of limited applicability to explore rapid perturbations of the global carbon cycle and the short-term response of the climate system in the distant past. In addition, the description of biogeochemical cycles in the original version of GEOCLIM does not include the nitrogen cycle and important processes of the sulfur and methane cycles. The model can thus not resolve the biogeochemical cycling and the oceanic redox dynamics under anoxic conditions. The new version, GEO-CLIM reloaded is designed to examine the feedbacks between past climatic conditions, marine geochemical dynamics and diagenesis at short geological timescales (<1 Myrs). It combines the climate model FOAM 3-D GCM with a vertically resolved diffusion-advection model of the global ocean, a pelagic biogeochemical model and a fully formulated diagenetic model. The new model describes the global cycles of carbon, nitrogen, phosphorus and oxygen. In addition, it includes important aspects of the coupled sulphur and methane cycles. The model design accounts for the limitations set by present-day computing power, the poorly constrained paleoboundary conditions, the lack of comprehensive data-sets for validation and, thus, the inherent need for extensive sensitivity studies and scenario runs. It therefore aims for a generic, simple, yet realistic representation of transport and transformation processes within the different reservoirs of the Earth system and across their interfaces. In addition, to the existing capabilities of GEOCLIM to resolve past climate change in response to changes in sea level, continental configuration, continental weathering fluxes, vegetation and hydrothermal or volcanic activity, the new model GEOCLIM reloaded now provides a detailed description of the biogeochemical cycling in the coupled ocean-sediment system to elucidate redox processes under extreme environmental conditions such as ocean anoxia, euxinia (e.g. Arndt et al., 2010) or methane release events. Therefore, it allows exploring the mechanisms and feedbacks that drive the response of the fully coupled climate system to past extreme carbon cycle perturbations. This paper provides first a detailed description of the coupled atmospheric, continental, oceanic and sediment modules. Next, the model performance is tested for the modernday ocean, since comprehensive, global benthic and pelagic data sets that allow for a careful and comprehensive modeldata comparison are only available for the present-day. Simulation results are compared with the WOCE Hydrographic Program (WHP) bottle data that contains hydrographic, nutrient and tracer data from the one-time and repeat WOCE Hydrographic Program (*hy1.csv files from WHPO web-site as of 27 May 2002). Also included are more than 2000 stations from pre-WOCE cruises. In addition, simulation results are also compared to previously published integrated global flux estimates based on simulation results or observations. Due to the more laborious and costly work flow associated with sediment core retrieval and interstitial data analysis, benthic observations are much less comprehensive and a database comparable to the WOCE WHP bottle data is not available. In addition, benthic data and modelling studies generally reveal a bias towards highly specific environments, such as sediments underlying upwelling areas and continental margins, or hydrate-rich sediments, while deep ocean sediments, which represent the major part of the ocean surface are not well monitored (Thullner et al., 2009;Regnier et al., 2011). Therefore, independent benthic and pelagic estimates of fluxes at the sediment-water interface are often inconsistent. Furthermore, model estimates of benthic fluxes and processes are strongly dependent on the weakly constrained parameterization of organic matter degradation, and are thus inherently associated with a significant degree of uncertainty. Therefore, diagenetic simulation results are only compared to well-established trends in observed depth profiles and to existing global estimates derived from data compilation or other model results.

Model description
The original version of the GEOCLIM model (Donnadieu et al., 2004(Donnadieu et al., , 2006Goddéris et al., 2008b) couples the biogeochemical ocean-atmosphere model COMBINE (COupled Model of BIogeochemical cycles aNd climatE, Goddéris and Joachimski (2004)) to the climate model FOAM 3-D GCM. GEOCLIM uses climatic reconstruction to estimate the spatially-resolved continental weathering fluxes delivered to the ocean and, thus, accounts for the impact of the continental setting on global biogeochemical cycles. In the new version GEOCLIM reloaded (Fig. 1), the box model COMBINE is replaced by a vertically resolved diffusion- advection model of the global ocean, which is coupled to a pelagic biogeochemical model and a fully formulated diagenetic model. The new model describes the global biogeochemical cycles of carbon, nitrogen, phosphorus and oxygen, as well as sulfate reduction, pyrite formation and sulfide reoxidation. It allows investigating the dynamic interplay between simulated oceanic conditions and diagenetic processes, as well as the evolution of the redox structure of the global ocean. The following sections describe the main features of the atmosphere Sect. 2.1 and the continent Sect. 2.2 modules, which are adapted from the former version of GEOCLIM. Next, the new ocean Sect. 2.3 and sediment Sect. 2.4 modules are described in detail. Figure 1 shows a schematic overview of the transport and transformation processes within the different modules of GEOCLIM reloaded, as well as the exchange fluxes through their respective interfaces. Table 1 provides a general overview of the differential equations for all species incorporated in GEOCLIM reloaded.

Atmosphere
Climatic reconstructions are provided by the general circulation model FOAM1.5 (for a detailed description see: http: //www.mcs.anl.gov/research/projects/foam/index.html). The atmospheric module of FOAM (Jacob, 1997) is a parallelized Table 1. Differential equations for the evolution of concentrations in the different reservoirs of the Earth system. The transport term dC i dt | trans summarizes the transport processes as described in Sects. 2.3.1 and 2.4. The biogeochemical reaction network is described in Sects. 2.3.2 and 2.4 and synthesised in Table 5. Exchange fluxes are detailed in Sects. 2.2.1 and 2.3.3. V r denotes the volume of the respective reservoirs (r = ES, ED, IS, ID, PS, PD). x, y and z indicate the carbon, nitrogen and phosphorus contents of organic matter, respectively. * Note that shallow carbonate precipitation R 2,1 and weathering fluxes, F weathering , only affect the epicontinental surface ocean ES and equal zero in the differential equations for the polar and inner ocean surface layers. + Hydrothermal fluxes, F mor , only affect the evolution in the deep inner ocean ID and equal zero in the differential equations for the deep polar and epicontinental oceans. The oceanic reservoirs of sulfate, calcium and boron are imposed and concentrations are constrained based on available observations, or simulated estimates.  with the upgraded radiative and hydrologic physics incorporated in CCM3 v. 3.2. The atmosphere module runs at R15 spectral resolution (4.5 • × 7.5 • ) with 18 vertical levels. The ocean module (OM3) is a z-coordinate ocean model that has been optimized for performance and scalability on parallel computers. OM3 contains 24 vertical layers, a 128×128 horizontal grid (1.4 • × 2.8 • ) and solves the conservation equations using second order differencing and a fully explicit time step scheme. A coupler calculates and interpolates the fluxes of heat and momentum between the atmosphere and ocean models, which are linked to the continent and sea-ice models (Jacob, 1997). Within GEOCLIM reloaded, FOAM can be used either with the coupled dynamic ocean OM3 or in a mixed-layer mode. In the latter, the atmospheric model is linked to a 50-m mixed-layer ocean, which parameterizes heat transport through diffusion. It is assumed that the Earth's orbit around the Sun is circular (eccentricity = 0) and the Earths obliquity is 23.5 • . This setting leads to an equal annual insolation for both hemispheres. Solar luminosity evolves through time according to the stellar evolution models (Gough, 1981). Computational cost renders a direct coupling of the atmospheric module FOAM within GEOCLIM reloaded infeasible. Therefore, an offline coupling approach is adopted, which operates as follows. For a given paleogeography, FOAM is run to a steady state with a number of different atmospheric carbon dioxide concentrations that cover the range of plausible values. A look-up table is then constructed on the basis of the resulting climatic reconstructions. Air temperature and continental runoff, required to quantify continental weathering fluxes, can then be extracted by linear interpolation for any given atmospheric CO 2 concentration on a 48 × 40 grid. In its present form, GEOCLIM reloaded accounts for the atmospheric concentrations of oxygen, O 2 , and carbon dioxide, CO 2 . It assumes that the atmospheric reservoirs of nitrogen gas is never limiting. The atmospheric concentration of O 2 is determined by the consumption of oxygen through the continental weathering of kerogen, F kerw (see Eq. 7) and its exchange with the surface layer of the ocean, F air|sea,O 2 (see Eq. 51): where V A denotes the volume of the atmosphere. Atmospheric CO 2 is determined by the input of subaerial volcanic CO 2 , F vol,CO 2 (see Sect. 2.2.2), the consumption of CO 2 in the continental weathering process of granite rocks, F silw , basalt, F basw and carbonates, F carbw (see Sect. 2.2.1), as well as by its exchange with the surface layer of the ocean, F air|sea,CO 2 (see Sect. 2.3.3):

Continents
The continental module of GEOCLIM reloaded uses the climatic reconstruction of FOAM to calculate spatially resolved weathering fluxes on a 7.5 • long × 4.5 • lat grid within the given paleogeographic setting. The terrestrial vegetation is simulated with a dynamic vegetation model LPJ (Donnadieu et al., 2009). Vegetation exerts only an indirect influence on continental weathering fluxes through its effect on global average temperature and runoff. The geochemical impact of land plants on weathering is thus not directly accounted for.

Continental weathering fluxes
The description of continental weathering fluxes, F weathering follows the approach adopted in the COMBINE model (Goddéris and Joachimski, 2004;Donnadieu et al., 2006) is extended here by an explicit description of nitrogen delivery to the ocean. The parameter values used in the continental weathering module of GEOCLIM reloaded are detailed in Table 2. The continental weathering of granitic, F silw , and basaltic lithologies, F basw , depends linearly on continental runoff as well as on land area. In addition, the kinetics of mineral-rock interactions is reflected in a temperaturedependency which follows an Arrhenius type law whose apparent activation energy is directly derived from field measurements (Oliva et al., 2003;Dessert et al., 2003): where F n denotes the total atmospheric CO 2 consumption flux through granitic and basaltic weathering, k n is a time constant, n pixel is the number of continental grid cells, R is the perfect gas constant, E n a is the apparent activation energy, area(i), T (i) and run(i) are the surface area, the mean annual ground air temperature and the mean annual runoff for the grid cell i, respectively. The constants k basw and k silw were calibrated using a present-day control FOAM simulation, assuming that basalt weathering accounts for 30 % of the total CO 2 consumption through silicate weathering (Dessert et al., 2003). Furthermore, total CO 2 consumption through silicate weathering was fixed to match the present-day estimate of 13.6 × 10 15 mmol yr −1 proposed by Gaillardet et al. (1999) and updated by Dessert et al. (2003). Carbonate weathering depends on the equilibrium concentration of calcium, Ca 2+ , which is determined by the soil pCO soil 2 and the solubility product of calcite.
The time constant k cw is calibrated so that the presentday weathering flux results in a total supply of alkalinity (carbonate + silicate weathering) to the oceans of 60.0 × 10 15 meq yr −1 . The equilibrium concentration of calcium, Ca 2+ eq , is calculated at each time step and each continental grid cell i and accounts for the dependence of the calcite solubility product on air temperature, T (i) and soil pCO 2 soil . The pCO 2 soil of the soil depends on the net primary production of the vegetation, the organic carbon content of the soil, its moisture, as well as on temperature and additional soil characteristics. None of these parameters can be estimated precisely for the geological past. Therefore, GEO-CLIM reloaded adopts a simple parameterization which provides first a maximum value of pCO 2 soil as a function of runoff, based on the relationship between net primary productivity and rainfall for modern ecosystems (Lieth, 1984): The maximum partial pressure of CO 2 in the soil pCO 2 soil (max) is then weighted by an air temperature factor, assuming higher soil pCO 2 under warmer climate due to enhanced soil organic carbon degradation (Gwiazda and Broecker, 1994): The dependency of soil CO 2 on temperature and runoff results in low productivity and low soil CO 2 in arid as well as cold areas.
The carbon release and the oxygen consumption by continental kerogen oxidation, F kerw is proportional to the continental runoff. It is thus assumed that the feedback of atmospheric oxygen concentration on kerogen weathering is negligible (Chang and Berner, 1999): The constant k kerw is fixed so that the present-day atmospheric oxygen concentration simulated by GEOCLIM equals 1 present atmospheric level (PAL). The behavior of phosphorus and nitrogen during weathering is very different from that of the major cations. On a geological time scale, most of the phosphorus is released through the dissolution of apatite present as a trace mineral in most silicate and carbonate lithologies Mackenzie, 2000, 2003). In GEOCLIM reloaded, the continental weathering flux of phosphate released by silicate rock dissolution is described according to (Guidry and Mackenzie, 2000): where H soil,i is the proton concentration in soil solution at each continental grid cell i. It is calculated as the proton concentration (mol m −3 ) of a solution equilibrated with the partial pressure of CO 2 in the soil. Furthermore, the model accounts for an additional phosphorus flux associated with the dissolution of carbonate (C:P ratio of 1000) and the oxidation of kerogen carbon (C:P ratio of 250). Nitrogen delivery to the ocean through continental weathering is proportional to the total carbon dioxide production through kerogen oxidation (Eq. 7) with a nitrogen/carbon ratio of 0.03 and through silicate rock weathering (Eq. 3) with a ratio of 0.001 (Berner, 2006). The total dissolved nitrogen weathering flux F nw is partitioned between ammonium, F nh4w , and nitrate fluxes, F no3w , based on the best fit of riverine data compiled by Meybeck (1982):

Volcanic degassing
Solid Earth degassing has been investigated using various tracers, such as the inventory of basaltic plateau production, seafloor accretion rate or the inversion of the sea-level (Kominz, 1984;Gaffin, 1987;Engebretson et al., 1992). Although different studies came to contrasting conclusions and no clear consensus has been reached so far, it has been recently shown that the evolution of degassing rates is more constant and lower throughout Earth history than originally thought (Cogne and Humler, 2004;Rowley, 2002). Therefore, GEOCLIM reloaded assumes hence that the volcanic degassing flux of carbon dioxide, F vol,CO 2 , is close to its present-day value (6.8 × 10 15 mmol yr −1 ).

458
S. Arndt et al.: GEOCLIM reloaded (v 1.0): a new coupled earth system Volume of an inner deep ocean layer m 3 A I · z z grid size of the inner ocean m 10

Ocean
The ocean module of GEOCLIM reloaded ( Fig. 1) combines a description of the global ocean circulation Sect. 2.3.1 with a formulation of biogeochemical transformation processes Sect. 2.3.2, which are detailed in Table 1. The exchange fluxes through the air-sea and the sediment-water interfaces are described further in Sect. 2.3.3. The model design is numerically efficient, yet it captures the large-scale features of the global ocean circulation with a limited number of free parameters, as well as the main characteristics of the redox and acid-base dynamics.

Transport
The global ocean circulation is described by a onedimensional, advection-diffusion model of the inner ocean, I , (latitude < 60 • ), as well as box models of high-latitude, P , (latitude > 60 • ) and epicontinental, E, (waterdepth <200 m) oceans (Fig. 2). The model structure is based on the well-established HILDA model (high-latitude exchange/inner diffusion-advection; Joos et al., 1991;Siegenthaler and Joos, 1992;Shaffer and Sarmiento, 1995). The three oceanic regions are represented by well-mixed surface layers (ES, IS, PS) overlaying a stratified inner deep ocean (ID) and well-mixed epicontinental and polar deep ocean reservoirs (ED, PD). The inner and polar oceans are connected by the deep thermohaline circulation, characterized by a global recirculation cell. The formation of deep, dense water in the surface layer of the polar ocean and its subsequent upwelling in mid-and low latitudes is described by a constant upwelling velocity w. This recirculation loop is completed by a bilateral horizontal exchange rate q between the deep inner ocean and the deep polar ocean, which captures the effect of entrainment processes. In addition, a turbulent and convective exchange between to two polar layers contributes to the deep recirculation cell and is described by a bidirectional exchange, proportional to a constant exchange velocity u. Vertical exchange depends on the of the global ocean circulation is discussed in the text and the ocean's geometry is provided in Table 3.  Table 3. turbulent diffusion. A constant vertical diffusion coefficient k that includes the effects of small-scale turbulence and mid-latitude ventilation processes is applied here. We thus neglect the effect of sloping isopycnal surfaces on vertical fluxes. This effect is difficult to constrain since it depends on the tracer gradients across the isopycnal surfaces, as well as on the vertical density structure of the ocean (e.g. Harvey, 1995;Harvey and Huang, 2001). The epicontinental upwelling is scaled to a bilateral horizontal exchange rate cq between the epicontinental deep ocean and the intermediate waters (50-200 m) of the inner ocean. Finally, the surface and deep epicontinental boxes are connected through an intense vertical exchange at an exchange rate cu. Assuming continuity (e.g., A I · w i = A P · w P ), the conservation equations for the concentration of species i in the oceanic reagion n, C n i , are given by:

Polar deep ocean (PD):
anaerobic oxidation of CH 4 37 - where the continental weathering fluxes, F weathering , is described in Sect. 2.2.1 and the formulations of the sediment/water exchange flux, F sed , the air/sea flux, F air|sea and the hydrothermal flux, F mor are detailed in Sect. 2.3.3. The geometric parameters (depth, D, surface area, A and volume, V ) and their present-day values are defined in Table 3. The reaction terms are equal to the sum of all transformation processes j affecting species i, each specified by the stoichiometric coefficient of species i, s j i and the kinetically controlled reaction rate R j . The solution of Eq. 16 requires specification of external boundary conditions at the upper and lower boundaries of the model domain. A Dirichlet (concentration) condition is applied at the thermocline, D TC , while, at the bottom of the inner ocean, D, the boundary condition guarantees a balanced between out-fluxes and in-fluxes: Transport parameters of the original HILDA model are summarized in Table 4. They have been carefully calibrated and validated on the basis of tracer distributions in the ocean (Siegenthaler and Joos, 1992;Shaffer, 1996). In the polar ocean, the thermocline is very weak or non existent and thus facilitates a vigorous exchange (129 Sv) by convection, mixing and Ekman pumping (Sarmiento, 1983;Luyten,1983). The ventilation of the inner ocean is driven by deep recirculation (89 Sv) and deep upwelling (7 Sv). The amount of vertical mixing in the inner ocean is generally very small because the energy input away from the ocean boundaries is low. In situ estimates from dye releases, microstructure measurements, as well as inverse approaches indicate that vertical eddy diffusion coefficients fall typically within the range 1×10 −4 −1×10 −5 m 2 s −1 (e.g. Ganachaud and Wunsch, 2000;Munk and Wunsch, 1998;Ledwell et al., 1998). The shallow epicontinental ocean is subject to a complex interplay of different forcing mechanisms such as tides, wind, or freshwater discharge, that act on timescales ranging from hours to month. Therefore, shallow-water dynamics cannot be fully captured by a simplified model of global ocean circulation. Here, we assume that the vertical mixing between the surface and deep layers of the epicontinental ocean is facilitated by wind and tides. The wind-driven epicontinental upwelling flux is set to 5 Sv and the bidirectional vertical exchange to 27 Sv, a value which prevents the development of anoxia in the epicontinental ocean under present-day conditions.

Biogeochemistry
The biogeochemical oceanic module provides a mechanistic description of the redox and acid-base dynamics in the ocean. Table 1 lists the mass conservation equations for all state variables. Table 5 summarizes the influence of each biogeochemical reaction, R j , on the transformation of each chemical species. Table 6. Hydrothermal fluxes as given by Elderfield and Schultz (1996).

Carbon export production
Most biogeochemical models of the present-day ocean quantify the export fluxes of particulate carbon from the surface to the deep ocean on the basis of complex recycling and transformation pathways within the euphotic zone (e.g. Sarmiento and Gruber, 2006). However, high data requirements, excessive computational demands, as well as the limited transferability of such approaches to the geological timescale compromise the application of similar models in a paleo-context. Therefore, although the ultimate limiter of pelagic primary production at geological timescales remains to be unequivocally identified (Codispoti, 1989;Smith, 1984;Tyrrell, 1999;Falkowski et al., 1998), paleo-models generally relate the whole-community export production directly to the availability of nutrients within the euphotic zone. Geochemists generally assume that phosphorus is the limiting nutrient (e.g. Holland, 1984;Ingall, 1994, 1996;Slomp and Van Cappellen, 2007) while biologists tend to favor nitrogen (e.g. McElroy, 1983;Codispoti, 1989;Falkowski, 1997). Here, we assume that the export production depends on both dissolved nitrogen, DIN ≡ NO − 3 + NH + 4 , and phosphate, PO 3− 4 , concentrations in the surface layers of the ocean. In addition, as production can be supported by N 2 fixation when dissolved nitrogen concentrations and the N/P ratio are low (Ohki et al., 1986(Ohki et al., , 1991, the total export production includes the contribution of both phytoplankton, R 1,1 , and nitrogen fixers, R 1,2 , according to Fennel et al. (2005): where K PO 3− 4 and K DIN denote the half-saturation constants for phosphate and DIN. µ max,1 and µ max,2 are the maximum production rates for phytoplankton and nitrogen-fixers, respectively. Production in the surface layers reduces the concentrations of dissolved inorganic carbon, T C , DIN and PO 3− 4 according to the stoichiometric Redfield ratio. In addition, DIN consumption is partitioned between nitrate and ammonium in proportion to the availability of ammonium (Table 1). The export flux of organic carbon, F POC , and inorganic carbon, F PIC , requires depth-integration over the euphotic layer. In addition, it is assumed that F PIC is directly related to the export production of organic carbon through the export rain ratio, r PIC:POC : There is no simple means of deriving the value of r PIC:POC as it is value is strongly dependent on the local ecosystem structure (Shaffer, 1993). Highly parameterized formulations on the basis of ambient nutrient and temperature conditions have been proposed (e.g. Maier-Reimer, 1993;Archer et al., 2000;Heinze et al., 1999;Sarmiento et al., 2002). Here, we use instead a formulation that relates the flux of inorganic carbon to the saturation state with respect to calcite, ca , (e.g. Ridgwell and Zeebe, 2005): where ρ PIC:POC is an adjustable parameter which reflects the observed global variation of r PIC:POC (e.g. Sarmiento et al., 2002). The power of 1.7 is chosen according to Opdyke and Wilkinson (1993). The saturation state of oceanic waters with respect to calcite, Ca is given by the ratio of the in-situ ion concentration product over the apparent solubility product of calcite, K * sp,Ca : where K * sp,Ca is calculated according to the empirical, temperature and salinity dependent function proposed by Mucci (1983). The pressure effect on K * sp is estimated from the molal volume and compressibility of calcite according to Millero (1995). GEOCLIM reloaded accounts also for the formation of carbonates in the epicontinental ocean. Shallow-water carbonates involve mostly biotic sources such as corals, carbonate banks or reefs and to a much lesser bioturbation coefficient m 2 yr −1 0.0028 * , 0.00013 * * φ porosity at the SWI -0.8 * , 0.8 * * * epicontinental sediments, * * deep sea/polar sediments extent, abiotic sources such as ooids (e.g. Schneider et al., 2006). The understanding of the physico-chemical and biological factors controlling the formation and composition of shallow carbonates is still incomplete. In GEOCLIM, shallow-water carbonate formation R 2,1 , depends on the saturation state of the epicontinental ocean with respect to aragonite, Ar , as well as on the total shelf area available for the formation of carbonate platforms, A sc , (Goddéris and Joachimski, 2004): where k sc is a scaled constant for the precipitation rate per unit surface area. The exponent determines the response of the carbonate precipitation rate to a change in epicontinental ocean CO 2− 3 concentrations. Possible values range from 1.0 to 2.8 (e.g. Ridgwell and Zeebe, 2005). The saturation state, Ar depends on the apparent solubility product of aragonite, K * sp,Ar , which is calculated according to Mucci (1983). The pressure effect on K * sp,Ar is estimated from molal volume and compressibility of aragonite according to Millero (1995). GEOCLIM reloaded does not account for temperature effects on shallow-water carbonate formation. In general, reaction kinetics of carbonate minerals in the natural environment are poorly known. Multiple experiments on corals and coraline algae revealed a dependence of biogenic calcification on temperature (e.g. Mackenzie and Agegian, 1989;Marshall and Clode, 2004). However, these relationships are species-dependent and often subject to uncertainties associated with experimental design. In GEOCLIM reloaded, shallow water carbonate precipitation integrates both abiotic and biotic carbonate precipitation. Because of the uncertainties associated to the temperature-dependency of shallow carbonate formation, the evolution of the coastal surface area, simulated shallow water temperatures and the generic nature of GEOCLIM reloaded, a simplifying description is adopted here. Although this formulation cannot resolve the carbonate phase mineralogy, it allows for a simple and scalable description of shallow water carbonate fluxes and their impact on global carbon cycling in the light of uncertainties associated with simulating past climate change.

Organic matter decomposition
Sediment trap observations across the global ocean show that a large fraction of the organic carbon is efficiently mineralized as it sinks through the water column. Field observations suggest that the recycling process is exponentially faster in the mesopelagic layers of the ocean than in the deeper ocean (e.g. Martin et al., 1987;Berelson, 2001;Schlitzer, 2002;Lutz et al., 2002). Thus, organic carbon mineralized in the upper ocean rapidly returns to the atmosphere, while it may take hundreds of years to recycle organic carbon that reaches the deep inner ocean layers. In most paleomodels, the mineralization of organic matter is parameterized on the basis of an exponential fit to the observed decrease in organic carbon fluxes (e.g. Shaffer, 1996;Fennel et al., 2005;Ridgwell et al., 2007). However, this approach does not provide insights into the evolution of the reactivity of the settling organic matter, a factor that has been advocated to play a crucial role in carbon cycling during past extreme events (e.g. Demaison and Moore, 1980;Sinninghé Damsté et al., 1998;Arndt et al., 2009). Therefore, we explicitly describe organic matter mineralization by means of the reactive continuum model (Boudreau and Ruddick, 1991). This approach integrates the effect of compound-specific reactivities on organic matter degradation rates. It represents the total amount of organic matter as the integral of the distribution function of reactive types, g(k,t), over the entire range of possible degradation rate constants, k. Each reactive type is degraded according to a first order rate law. Boudreau and Ruddick (1991) proposed the Gamma Distribution as an initial distribution of reactive types, g(k,0): where the rate constant k is a continuous variable and (ν) is the Gamma function (e.g. Abramowitz and Stegun, 1972). The free parameters a and ν determine the shape of the initial distribution of organic matter reactivities. The parameter a describes the average lifetime of the more reactive components of the spectrum, while ν defines the shape of the distribution near k = 0. In general, low ν values reflect a predominance of unreactive types, while high ν values characterize a more uniform distribution. Derivation and substitution leads to the rate of organic matter degradation in the water column, R 3 (Boudreau and Ruddick, 1991): where a ocean denotes the age of organic carbon as it settles through the water column: where z is the depth and w POC denotes the average sinking velocity of organic particles in the ocean. Observed POC depth profiles show a large variability. In general, a stronger attenuation of POC fluxes is observed in the polar ocean (e.g. Berelson, 2001;Francois et al., 2002). This efficient mineralization is often related to the weak processing of organic matter in the euphotic layer, as well as the labile and loose particle packing associated with the strong, seasonal diatom blooms in this area. On the contrary, POC exported from carbonate-dominated, warm oligotrophic waters is tightly packed and more extensively processed by complex food webs in the euphotic zone. Therefore, we assume that organic matter exported from the polar ocean is more reactive than organic carbon from the epicontinental and inner euphotic layers (Fig. 3). The free parameters a and ν for the polar and the epicontinental/inner oceans are extracted from the depth evolution of observed organic matter fluxes in the global ocean (Fig. 3) by converting flux observation to concentrations, assuming an average particle sinking speed w POC of 50 m d −1 (e.g. Berelson, 2001).
In poorly ventilated environments, the degradation of organic matter is coupled to the sequential utilization of terminal electron acceptors (TEAs) other than oxygen, used in the order of their decreasing Gibbs energy yield (O 2 . After all TEAs are consumed, the remaining organic matter may be degraded by methanogenesis, a process which is independent of an external TEA and is sustained by the fermentation products of organic matter. The rate of TEA consumption assumes a Michaelis-Menten type relationship with respect to the TEA concentration (e.g. Boudreau, 1997). In addition, a specific metabolic pathway can be inhibited by the presence of energetically more favorable TEAs. Here, we only include aerobic oxidation R 4 , denitrification, R 5 , sulfate reduction, R 6 and methanogenesis, R 7 , since the contribution of metaloxide reduction pathways is generally small (Thullner et al., 2007(Thullner et al., , 2009: where K i denotes the half-saturation constant for the TEAs (i), oxygen O 2 , nitrate NO − 3 and sulfate SO 2− 4 .

Secondary redox reactions
Ammonium, total sulfide, TS, and methane, CH 4 , produced during organic matter degradation can be oxidized to nitrate, sulfate and carbon dioxide via a set of secondary redox reactions which sustain a continuous recycling of these redoxsensitive elements. Here, we consider nitrification, R 8 , total sulfide, T S ≡ HS − + H 2 S, oxidation by O 2 , R 9 , methane oxidation by O 2 , R 10 and anaerobic oxidation of methane by SO 2− 4 , R 11 . The rate of secondary redox reactions are formulated using bimolecular rate laws with a Monod-dependency with respect to the oxidant: where K i,redox denotes the Monod constants for oxygen O 2 , and sulfate SO 2− 4 .

Iron sulfide formation
The formation of iron sulfides in the water column and its subsequent burial in marine sediments can represent an important sulfur sink under anoxic and euxinic conditions (e.g. during oceanic anoxic events). Since the global iron cycle, which involves various transformations between particulate, dissolved and complexed phases, is not explicitly resolved in the model, we adopt a simplified description for iron sulfide precipitation. Following Slomp and Van Cappellen (2007), we make the simplifying assumption that the availability iron in the water column never limits the rate of iron sulfidization. The rate of pyrite formation, R 12 , is thus merely controlled by the concentration of total sulfide, T S : where k Fe,S denotes the rate constant of iron sulfide formation. Although this simple formulation does not provide a mechanistic description of the complex process of iron sulfide precipitation, it is nevertheless suitable for our purposes.

The carbonate system
The carbonate system is constrained by the dissolution and hydration of CO 2 , the following acid-base equilibria: and their corresponding mass action laws are: where K * 0 , K * 1 and K * 2 are the apparent constants for CO 2 solubility and dissociation of carbonic acid, HCO − 3 , and bicarbonate, CO 2− 3 , respectively. Thus, the carbonate system is mathematically described by three equations and four unknowns and, as such, underdetermined. The total alkalinity, T A , can be used to close the system: where the composite variables total dissolved carbon, T C ≡ CO 2 * + CO 2− 3 + HCO − 3 , total dissolved sulfide, T S ≡ HS − + H 2 S, total dissolved borate, T B ≡ B(OH) − 4 +B(OH) 3 , and total alkalinity, T A , are conserved during acid-base reactions. Contributions from ammonium, fluorine, phosphate, silicate and other minor species to T A are neglected here. Each species concentration contributing to T A , T C , T S and T B can be written as a fraction of the conserved quantities: , and the dissociation constants of water, K * 3 , sulfide, K * 4 , and boric acid, K * 5 , are based on a total concentration scale and, therefore, vary with temperature, salinity and pressure according to empirically derived equations (Dickson and Millero, 1987;Millero, 1995). The resulting system of equations can be solved by finding the chemically meaningful root of the function: Following the computationally fast approach proposed by Follows et al. (2006), the H + concentration can be determined by combining the expression for T C with the equilibrium expressions Eqs. (43) and (44). The H + concentration of the previous time step is used as a first guess to estimate the carbonate alkalinity, C A : C A is then used to evaluate the quadratic function in H + : If necessary, the new guess is subsequently used to calculate a new value for C A and the procedure is iteratively repeated until the root of f (H + ) has been found with sufficient accuracy (f (H + ) < 10 −6 ). For buffered systems, iteration is generally not needed.

Calcite dissolution
The mechanisms that control the dynamics of carbonate dissolution in the ocean are still elusive. Global observations indicate that the depth of the lysocline coincides with the saturation horizon (e.g. Sarmiento and Gruber, 2006), suggesting that thermodynamic constraints are a dominant control on calcium carbonate preservation. Yet, in-situ experiments carried out in the North Pacific (e.g. Berger, 1967) and laboratory studies (Morse and Berner, 1972;Berner and Morse, 1974) reveal that dissolution rates increase below the saturation horizon, emphasising a kinetic control on calcium carbonate dissolution. In addition, Iobservational evidences point to a partial dissolution of the sinking calcite in the supersaturated upper 1000 meters of the ocean (Milliman, 1993). Here, we follow the description proposed by Morse and Berner (1972) and, therefore, the model does not account for the dissolution of calcite above the lysocline, due to the uncertainty associated with this hypothesis. The rate of CaCO 3 dissolution , R 13 , is described by a non-linear rate law depending on the degree of the undersaturation of seawater with respect to calcite (e.g. Morse and Berner, 1972;Keir, 1980): where k diss,CaCO 3 is the rate constant of carbonate dissolution. The exponent n is a measure of how strongly the rate responds to a change in CO −2 3 concentration and, thus, of the buffering efficiency of the ocean with respect to changes in atmospheric CO 2 . The saturation state, Ca is defined as the ratio of the in-situ ion concentration product over the apparent solubility product of calcite, K * sp,Ca : K * sp,Ca is calculated according to an empirical, temperature and salinity dependent function proposed by Mucci (1983). The pressure effect on K * sp,Ca is estimated from the molal volume and compressibility of calcite according to Millero (1995).

Exchange fluxes Hydrothermal fluxes, F mor
Submarine hydrothermal activity at mid-ocean ridges has a global impact on ocean chemistry, but its evolution over geological timescales remains poorly quantified. Here, we use the primary axial high-temperature hydrothermal chemical fluxes reported by Elderfield and Schultz (1996) and German and Damm (2003) to estimate the hydrothermal input, F mor , of carbon dioxide, total alkalinity, sulfide and methane at mid-ocean ridges (Table 6). The net air-sea flux, F air|sea , of a gas, i, is quantified according to the stagnant film model (e.g. Whitman, 1923): where A denotes the ice-free surface area of the oceanic region in contact with the atmosphere (A E , A I , A P f ree ), k w,i denotes the piston velocity, C 0 i is the saturation concentration and C w i is the concentrations of the gaseous species. At equilibrium and assuming ideal behavior, C i is, according to Henry's law, directly proportional to the partial pressure of the corresponding gas p a i in moist air: The solubility of CO 2 , S CO 2 , is calculated via a solubility function, F CO 2 , which includes a correction for non-ideal behavior (Weiss, 1974). For oxygen, the saturation concentration, C 0 O 2 is obtained from the Bunsen coefficient, assuming near-ideal behavior (Weiss, 1974). The piston velocity k w,i is determined according to the formulation proposed by Wanninkhof (1992): assuming an average global windspeed, U , of 7.5 m s −1 at a nominal height of 10 m a.s.l. The Schmidt numbers, Sc i , are calculated as a function of temperature and salinity using the coefficients given in Wanninkhof (1992) and Keeling et al. (1998) for CO 2 and O 2 , respectively.

Sediment/water flux, F sed
The bidirectional exchange flux of dissolved species through the sediment-water interface (z = 0) is calculated by: where φ rsed denotes the porosity of the sediments underlying the oceanic region rsed, D i is the molecular diffusion coefficient of species i and D bio is the bioturbation coefficient. The concentration gradient through the sediment-water interface, dC r i dz , is given by the length-weigthed concentration difference between the porewater (see Eq. 2.4) and the overlaying bottom water.

Sediment
A sediment column, rsed, underlies each of the three deep oceanic regions, ED, ID, PD (Fig. 1). A vertically resolved, fully-formulated diagenetic model calculates the transport and the biogeochemical transformation processes at each grid point within these sediment columns as well as the sedimentary burial and recycling fluxes at the model boundaries (see Eq. 2.3.3). Tables 1 and 5 give a general overview of the reactive-transport equations incorporated in the benthic module of GEOCLIM reloaded. The one-dimensional mass conservation equation for dissolved or solid species is given by (e.g. Boudreau, 1997): is the concentration of species i in the sediment column rsed, D bio is the bioturbation coefficient; σ rsed i is equal to the porosity, φ rsed and to (1 − φ rsed ) for dissolved and solid species, respectively; D i denotes the molecular diffusion coefficient of dissolved species i (D i = 0 for solid species); ω rsed is the sediment burial rate and s j i denotes the stoichiometric coefficient of species i in the kinetically controlled reaction j , with rate R j . Porosity, φ rsed , is assumed to be constant over the entire sediment column. The activity of infaunal organisms in the upper decimeters of the sediment causes random displacements of sediments and porewaters and is simulated using a diffusive term (e.g. Boudreau, 1986), with a depth-dependent bioturbation coefficient D bio : where D bio (0) denotes the bioturbation coefficient at the SWI and erfc is the complementary error function. GEOCLIM reloaded assumes that bioturbation shuts down once bottom waters become anoxic (O 2 = 0.0 mmol m −3 ). The pumping activity by burrow-dwelling animals, the so-called bioirrigation, is neglected here as it is difficult to constrain on geological timescales. At the sediment-water interface, a Dirichlet (concentration) boundary condition is applied for dissolved species and a Neumann (flux) boundary condition is used to constrain the deposition fluxes of particulate organic carbon (POC) and inorganic carbon (PIC). The concentration of dissolved species and the solid deposition fluxes at the upper boundary are provided by the ocean module of GEOCLIM reloaded. At the bottom of the sediment column, a no-flux condition is applied.
The biogeochemical reaction network in the sediment, j s j i R j , includes primary and secondary redox processes, adsorption processes, carbonate dissolution and iron sulfide precipitation reactions, as well as the acid-base equilibria of the carbonate, sulfide, and boric acid systems. In the sediments, process formulations for the relative contributions of the different metabolic pathways of organic matter degradation, R 4 − R 7 , secondary redox reactions, R 8 − R 11 , and carbonate dissolution, R 13 follow the rate formulations described in Sects. 2.3.2 and are not repeated here. Similarly, the acid-base equilibria are also included according to the approach used for the ocean module. The following paragraph describes the process formulations of sedimentary organic matter degradation, R 3 , and iron sulfide formation, R 12 , that slightly differ from the pelagic formulations. Additional processes specified in the benthic module include adsorption equilibria for ammonium and phosphate, as well as apatite formation, R 14 . In the sediment, the bulk rate of organic matter degradation, R 3 is again determined by the reactive continuum model (Boudreau and Ruddick, 1991). However, the rate description reads here: where a sed (z) denotes the age of organic matter in the sediment: where ω denotes the sediment burial rate. In the bioturbated zone, a sed is assumed constant and equals a bio . Precipitation of iron sulfides, R 12 , only occurs in the oxic and suboxic layer of the sediment, where the upward diffusing sulfide produced in the anoxic sediment layers reacts with reactive iron hydroxides. The limited availability of reactive iron in anoxic sediment layers is assumed to inhibit iron sulfide formation. In the oxic and suboxic sediment layers, iron sulfide formation depends linearly on sulfide concentration: Adsorption of phosphate and ammonium is a controlling factor of nutrient porewater concentrations and benthic recycling fluxes (e.g. Krom and Berner, 1980; where PO 4 3− and NH 4 + denote the concentrations of adsorbed phosphate and ammonium. The adsorption of phosphate occurs predominantly on iron hydroxides and, therefore, is strongly dependent on the redox conditions in the sediment (e.g. Krom and Berner, 1980). In anoxic sediments, phosphate is released during the reduction of ferric oxyhydroxides and the ratio between adsorbed and dissolved phosphate is considerably less (2:1) than in the oxic and suboxic layers of the sediment (250 to 400:1) (e.g. Krom and Berner, 1980). Therefore, the adsorption coefficient for phosphate, K ads,PO 4 3− is higher when NO − 3 > 0 and lower in the anoxic layers.The precipitation of carbonate fluorapatite, R 14 , also reduces the benthic flux of phosphate (e.g. Jahnke et al., 1983;Van Cappellen and Berner, 1988). This process is described according to the approach proposed by Van Cappellen and Berner (1988), which expresses the rate as a function of the departure from the thermodynamic equilibrium PO 3− 4,eq : (64) where k ap denotes the effective rate coefficient of apatite precipitation and the equilibrium concentration PO 3− 4,eq is chosen according to Van Cappellen and Berner (1988).

Numerical solution
The different modules of GEOCLIM reloaded are solved stepwise within one model time step. First, climatic reconstructions are calculated on the basis of the simulated atmospheric carbon dioxide concentration. These reconstructions are then used to quantify the weathering fluxes to the epicontinental ocean. In the ocean module, transport and reaction equations are treated separately according to the operator splitting approach (Steefel and MacQuarrie, 1996). The continuity equations for the well-mixed ocean boxes are solved numerically using an Euler forward scheme, while the advection-dispersion equation for the vertically resolved, inner deep ocean is solved by applying again an operator splitting technique. At each time step of the numerical integration, the dispersion term is first discretized using the   Ruardij and Raaphorst (1995), Krom and Berner (1980)   to be continued Table 9.   Dickson (1990), Millero (1995) f(T,p,S) Dickson (1990), Millero (1995) semi-implicit Crank-Nicholson scheme (e.g. Press, 1992), followed by the calculation of the advective transport, using a 3rd order accurate total variation diminishing algorithm with flux limiters (e.g. Regnier et al., 1998). The reaction network is subsequently solved using the Euler scheme. The reactiontransport equation is discretized on a regular grid with a resolution z ocean of 10 m. The calculated concentrations at the last node right above the seafloor in the inner ocean and in the deep epicontinental and polar ocean boxes are then used as upper boundary conditions for the sediment module. Here, transport and reaction equations are also solved sequentially in a manner identical to that used for the inner ocean. The continuity equation (56) is however discretized on an uneven grid (e.g. Boudreau, 1997): where z(n) is the depth of the n-th grid point, L denotes the length of the simulated sediment column, ξ n is a point on a hypothetical even grid, and ξ c is the depth relative to which z(n) is quadratically distributed for ξ n ξ n and linearly distributed for ξ n ξ n . L and ξ n are chosen so that that grid size, z sed increases downcore from the sediment-water interface ( z(0) = 3 mm) to the maximum simulated sediment depth, nmax ( z(nmax) = 1 cm). The ocean transport model uses a time step t = 0.1 yr, while the biogeochemical process rates in the ocean are solved with a timestep t ocean = 10 −3 yr and the sediment algorithm is run with an automatically adjusting timestep t sed ≤ 10 −3 yr on the unevenly spaced grid. GEOCLIM reloaded is a comparably efficient model, due to its low spatial resolution of the global ocean and the offline coupling between the oceanic, continental and sedimentary sub-modules with the threedimensional climate module FOAM. The model can be run on a standard personal computer on a single CPU. Numerical efficiency is limited by the large range of characteristic temporal and spatial scales involved. Especially the diagenetic model requires a small timestep and, therefore, reduces the overall model performance. The model efficiency is approximately 100 kyrs per CPU day for the fully coupled model and reaches 1 Myrs per CPU day without the sediments. The model is thus not as efficient as a simple box model and is rather comparable to the computational efficiency of lowdimensional EMICs (e.g. Bern 2.5D). Assuming constant atmospheric CO 2 levels, the continent-ocean-sediment reaches a steady state after approximately 10 kyrs. The simulation time required to reach a "quasi" equilibrium for the coupled model is much longer and takes a few 100 kyrs depending on the initial conditions.

Model set-up
Steady-state Earth system simulations for present-day conditions are run to evaluate the overall performance ofthe new ocean and sediment modules within GEOCLIM reloaded. Since the model is designed to investigate extreme climate periods in Earth history, model simulations are realized with a generic set of model parameters (Table 8) and calibration effort is restricted to a minimum. Atmospheric carbon dioxide and oxygen concentrations are kept constant at their present day values of 365 ppm and 209 460 ppm, respectively. The constant CO 2 concentrations are used to determine time-invariant temperatures, runoff and continental weathering fluxes, which conduct the coupled continentocean-sediment system into a steady-state that can be compared with present-day conditions. Table 10 compares the global and regionalized organic and inorganic carbon export fluxes simulated by GEOCLIM reloaded with published estimates based on field observations or model predictions. Organic carbon export fluxes exert a dominant control on pelagic and benthic biogeochemical dynamics, while inorganic carbon export fluxes play a crucial role for the carbonate system, the alkalinity budget of the ocean, and the calcium carbonate compensation depth. Thus, considerable effort has been devoted over the past decade to quantify the global carbon export fluxes. For POC export flux estimations, a variety of different approaches, including satellite color (e.g. Laws et al., 2000), benthic oxygen fluxes (e.g. Jahnke, 1996), theoretical relationships (e.g. Dunne et al., 2005;Lutz et al., 2002), nutrient removal (e.g. Sarmiento et al., 2002) or global biogeochemical models  have been applied. Although published estimates of the total global POC export flux vary between 5.0 and 12.9 Pg C yr −1 (Table 10), a broad consensus has progressively converged towards an estimate of 10 ± 3 Pg C yr −1 . The simulated global POC export flux of 8.54 Pg C yr −1 falls within this range. In addition, the regionalized estimates of GEO-CLIM reloaded agree well with the sattelite-derived export fluxes from Dunne et al. (2007). In particular, the simulated POC export flux per unit surface area is higher in the polar ocean (4.54 mol C m −2 yr −1 ) than in the inner ocean (1.13 mol C m −2 yr −1 , Fig. 4a), in agreement with prominent maxima in the subpolar to polar regions and minima in the inner ocean gyres reported by Sarmiento and Gruber (2006) and Dunne et al. (2007). Published PIC export fluxes are either based on sediment trap data, estimated from seasonal alkalinity drawdowns, or based on molar PIC/POC rain ratios and organic carbon export fluxes. The three methods are subject to considerable uncertainty (Berelson et al., 2007) and there is thus no well established agreement concerning the magnitude of global PIC export fluxes. Estimates of global PIC export fluxes show a large variability (Table 10) and fall within the broad range between 0.4 and 1.8 Pg C yr −1 . In GEOCLIM reloaded, the PIC export (0.48 Pg C yr −1 ) is calculated on the basis of the simulated PIC/POC rain ratios (Eq. 24) and POC export fluxes. It falls within the   lower range of published values. Nonetheless, it agrees well with the recent regional estimates reported by Dunne et al. (2007), characterized by maximal PIC export in the carbonate-dominated equatorial regions and minimal export fluxes in polar regions where surface ocean productivity is dominated by diatoms (Table 10, Fig. 4b). The quantification of shallow-carbonate production not only requires the determination of production rates, butalso of global reef and platform surface areas. Published values are sparse (Milliman, 1993;Milliman and Droxler, 1996;Vecsei, 2004;Schneider et al., 2006) and indicate that the rate of shallow-water carbonate accumulation could be as high as ∼ 60 % of the global PIC export production. The simulated shallow-water carbonate accumulation rate of 0.25 Pg C yr −1 agrees well with this estimate.

Mineralization
The preferential consumption of labile compounds during the settling and burial of organic matter leads to a continuous decrease in reactivity and, therefore, in mineralization rates from the surface ocean down to the deep sediment (e.g. Hedges and Keil, 1995;Boudreau and Ruddick, 1991). Simulation results show that 82% (7.05 Pg C yr −1 ) of the global POC export production is mineralized within the water column (Fig. 4). Benthic mineralization removes another 14% (1.22 Pg C yr −1 ), a value which is equivalent to 79% of the global POC sedimentation flux (1.54 Pg C yr −1 ). Therefore, only 4 % (0.32 Pg C yr −1 ) of the global organic carbon export production is ultimately buried in sediments. The water-column estimates compare well with the 7 Pg C yr −1 proposed by Dunne et al. (2007) for the mineralization in the world ocean and the 6.6 Pg C yr −1 reported by Hansell and Carlson (2002) for the deep ocean mineralization. In contrast, the simulated global benthic mineralization rate of 1.22 Pg C yr −1 is lower than previously published estimates. Thullner et al. (2009) calculated a value of 2.52 Pg C yr −1 on the basis of a fully-formulated diagenetic model that resolves the global hypsometry. Experimentally-determined, total oxygen utilization (TOU) rates, assumed to be equivalent to the mineralization flux, are also higher and fall within the range of 2.28-2.76 Pg C yr −1 (Jørgensen, 1983;  Smith and Hollibaugh, 1993;Canfield et al., 2005). This discrepancy is partly due to the low resolution at which the seafloor is resolved in GEOCLIM reloaded (3 sediment columns). In addition, the simulated global POC deposition of 1.54 Pg C yr −1 , falls at the lower end of the spectrum of sedimentation fluxes (2.3±0.9 Pg C yr −1 ) reported by Dunne et al. (2007) and could also explain why the calculated benthic mineralization rates are lower than previous estimates. Note that the simulated POC burial flux of 0.32 Pg C yr −1 falls within the the range (0.29 ± 0.15 Pg C yr −1 ) reported by Dunne et al. (2007). In addition, the reasonable fit between carbon, nitrogen and phosphorus Vvertical concentration gradients in the ocean and observations (see following sections) indicates that GEOCLIM reloaded provides a realistic estimate of pelagic and benthic mineralization rates. Figure 5 shows that simulated oceanic depth profiles compare well with global observations from the WOCE Hydrographic Program. The oxidation of organic matter leads to a net consumption of oxygen (Fig. 5a), and a release of phosphate (Fig. 5b) and ammonium (not shown) below the euphotic surface layer of the oceans. The ventilation of the deep ocean results in a rapid oxidation of ammonium to nitrate and, therefore, nitrate (Fig. 5c) is the dominant form of dissolved inorganic nitrogen. In the epicontinental and polar oceans, comparably high and constant oxygen concentrations are maintained by the intense vertical mixing (Fig. 5a), a result in agreement with field observations. In the deep inner ocean, oxygen concentration increase again below the oxygen minimum zone due to the entrainment of oxygen-rich,   deep polar waters. A second decrease close to the sediment water interface can be identified in the simulation as a result of the high benthic oxygen demand. In contrast to the well-ventilated ocean, the supply of oxygen in marine sediments is strongly limited by transport and, therefore, these systems are generally characterized by a sequential depletion of terminal electron acceptors (TEAs) used in the process of organic matter mineralization. The sequential depletion of TEAs used during primary redox reactions shapes, in concert with the re-oxidation of reduced components, the redox-zonation of the sediments. Figure 6 illustrates simulated depth profiles of POC (Fig. 6a), oxygen (Fig. 6b), nitrate (Fig. 6c) and sulfate (Fig. 6d), as well as of the characteristic by-product of organic matter mineralization ammonium (Fig. 6e), phosphate (Fig. 6f) and sulfide (Fig. 6g). Epicontinental sediments receive compara-bly high fluxes of labile POC fluxes (Fig. 6a), which control the drawdown of oxygen and nitrate concentrations within the upper first centimeters and the subsequent consumption of sulfate ( Fig. 6b-d), leading to a significant build up of ammonium and sulfide in the anoxic layers ( Fig. 6e and g). The reduced species diffuse towards the sediment-water interface and are efficiently re-oxidized within the very narrow oxic zone, thus contributing significantly to the total benthic oxygen consumption. In addition, the limited availability of iron oxides in the anoxic sediment reduces the adsorption of phosphate, resulting in comparably high diffusive phosphate fluxes to the water column (Fig. 6f). Because of the longer transit time of organic matter through the deep inner ocean, underlying sediments receive lower amounts and much less reactive POC (Fig. 6a). Therefore, oxygen penetrates deeper into the sediment (Fig. 6b), favoring nitrification and, thus, higher nitrate concentrations. The latter are also sustained by comparatively higher concentrations in the overlying water column of the deep ocean (Fig. 6c). As a result, sulfate reduction is inhibited throughout most of the sediment column and sulfate concentrations only show a very small decrease below the depth of nitrate penetration (Fig. 6d). The build-up of ammonium and sulfide in the anoxic sediment is limited and these species are entirely re-oxidized in the expanded oxic layer ( Fig. 6e and g). Phosphate concentrations remain low due to the adsorption in the oxic and suboxic zones (Fig. 6f). In the polar ocean, the intense pelagic organic matter mineralization (Fig. 3) reduces even further the magnitude and reactivity of the POC deposition fluxes. Oxygen and nitrate concentrations remain therefore almost constant throughout the entire sediment column ( Fig. 6b and c) and sulfate reduction is completely inhibited ( Fig. 6d and g). The ammonium produced during organic matter degradation is efficiently converted to nitrate by nitrification, a process that occurs everywhere within the oxygenated sediment (Fig. 6c and e). Because of the low recycling rates and the enhanced adsorption onto iron oxides in the sediment column, benthic phosphate concentrations are also lower than in the other oceanic regions. These simulated global patterns are in overall qualitative agreement  with observations and results from global diagenetic models (e.g. Hensen et al., 2006;Thullner et al., 2009). The simulated oxygen penetration depths for the inner ocean are nevertheless deeper than those reported by Meile and Cappellen (2003) and Thullner et al. (2009) for average ocean conditions, but analysis of porewater profiles shows a high variability in oxygen penetration depths. For instance in sediments underlying inner ocean gyres, oxygen concentrations decrease in the topmost layers and then remain constant at non-zero values below 20-40 cm (e.g. Murray and Grundmanis, 1980;Sarmiento and Gruber, 2006;Fischer et al., 2009). Sachs et al. (2009) also report strongly variable oxygen penetration depth, between <15 cm and several decimeters in sediments from the Southern Ocean. Table 11 summarizes the regionally-resolved diffusive flux of oxygen, nitrate and phosphate through the sediment-water interface. The values of the benthic exchange fluxes mostly reflect the magnitude and the quality of the POC deposition fluxes. Simulated oxygen uptake is highest in the epicontinental sediments and decreases in the inner ocean and polar ocean sediments. The magnitude of the simulated diffusive oxygen fluxes com-pares well with global observations. According to benthic flux estimates from a number of different sites in the global ocean, oxygen fluxes vary between 600-900 mmol m −2 yr −1 in shallow oceanic regions and 100 mmol m −2 yr −1 underneath inner ocean gyres (Jahnke, 1996). Measured benthic oxygen fluxes across the continental shelf in the Northeast Pacific increase from 0-200 mmol m −2 yr −1 in deep ocean sediments to 200-1000 mmol m −2 yr −1 on the continental slope, as a function of POC deposition flux and bottom water oxygenation . In addition, simulated oxygen fluxes from a global diagenetic model decrease from 1000 mmol m −2 yr −1 in shallow sediments to 100-200 mmol m −2 yr −1 in the deep ocean (Thullner et al., 2009). Nitrate fluxes into the sediment are also higher in the shallow epicontinental ocean than in the inner ocean. In contrast, polar ocean sediments act as a small nitrate source to the overlaying water column. Observed nitrate fuxes generally reveal a strong variability, which reflects the complex nature of benthic nitrogen cycling. In the shallow ocean, nitrate depth-profiles are generally characterized by a shallow local maximum, whose location and magnitude depends bservations using the in-situ temperature, pressure, T C and T A (Fig. 7) The colour code indicates the number of times a concentration s reported for a given water depth.  (Fig. 7) The colour code indicates the number of times a concentration is reported for a given water depth.

Nutrient cycling and redox dynamics
Simulated steady state depth profiles of saturation state (Ω calcite ) in the epicontinental, inner and polar ocean (b results are compared to calcite saturation calculated based on WOCE Hydrographic Program (*hy1.csv files from s of May/27/2002) observations using the in-situ temperature, pressure, T C and T A (Fig. 7). The colour code ind times a concentration is reported for a given water depth.   (Fig. 7). The colour code indicates the number of times a concentration is reported for a given water depth.
on the intensity of nitrification in the narrow oxic zone of the sediment. As a consequence, observed nitrate fluxes show a large variability between 1000 mmol m −2 yr −1 and −250 mmol m −2 yr −1 (Middelburg et al., 1996). In the deep ocean, the low ammonium release combined with the presence of an expanded oxic zone leads to much lower (−5-igure 10: Simulated steady state depth profile of a) particulate inorganic carbon, b) total dissolved arbon, c) total alkalinity, d) pH and e) pore water saturation state with respect to calcite in picontinental, inner and polar ocean sediments. 5 mmol m −2 yr −1 ) and less variable nitrate fluxes (Middelburg et al., 1996). In deep-sea sediments, the increased availability of oxygen and nitrate inhibits sulfate reduction and diffusive sulfate fluxes are thus close to zero. Noticeable sulfate diffusion is thus only simulated for the shallow epicontinental ocean sediments, in agreement with simulation results of Thullner et al. (2009). The benthic phosphate fluxes simulated by GEOCLIM reloaded are equal to 0.7 mmol m −2 yr −1 , 0.1 and 0.01 mmol m −2 yr −1 for the epicontinental, the inner and the polar ocean, respectively. These estimates fall within observed fluxes in the eastern South Atlantic (0-2 mmol m −2 yr −1 ), although measured fluxes may locally reach values as high as 8 mmol m −2 yr −1 . Figure 7 compares the oceanic depth profiles of dissolved inorganic carbon (T C ), total alkalinity (T A ) and pH with global observations from the WOCE Hydrographic Program. In general, simulation results compare well to global observations. The vertical distribution of T C is mainly controlled by the organic carbon pump and to a lesser extent by air-water exchange and the inorganic carbon pump (e.g. Sarmiento and Gruber, 2006). In the vertically resolved ocean, T C concentrations rapidly increase below the thermocline where net mineralization rates are maximal. The same trend is observed in global observations which substantiate further the representation of the organic carbon pump in GEOCLIM reloaded. Observations and model results for the epicontinental waters reveal high T C concentrations. In the polar ocean, gas exchange tends to remove T C from the carbonrich waters and, therefore, reduces the vertical T C gradient produced by biological processes. The vertical distribution of T A is affected to a small extend by the regeneration of organic nitrogen, but is dominantly controlled by the efficiency of the carbonate pump. Primary production and nitrification reduce T A in the surface ocean layers, while carbonate dissolution increases T A in the deep ocean layers. The shape of the inner ocean T A profile shows a gradual increase throughout the water column and is consistent with the depth variations of the saturation horizons and the carbonate dissolution which are discussed below. The vertical distribution of oceanic pH reveals a sharp drop below the surface layer immediately followed by a local minimum. In the modern, well-ventilated ocean, this local pH minimum can be attributed to maxima in aerobic mineralization and secondary redox-reactions rates (Jourabchi et al., 2005;Soetaert et al., 2007). In the deep ocean, calcite dissolution buffers the drop in pH. The epicontinental and polar oceans reveal a similar pH profile. Despite considerable scatter, the same general features can also be identified in the global ocean pH data. The model performance can be assessed further by comparing the profiles of the inorganic carbon species HCO − 3 , CO 2− 3 and CO 2 * (Fig 8). A realistic description of the carbonate system is decisive for the performance of the earth system model, as the surface ocean CO 2 * concentration controls the air-sea gas exchange and the CO 2− 3 concentration determines the saturation state of oceanic waters with respect to calcium carbonate. In particular, the CO 2− 3 ion distribution results from the competing effects of mineralization and carbonate dissolution. In the upper ocean, aerobic organic matter mineralization increases T C concentrations but exerts, due to its coupling with nitrification, a negative net-effect on T A concentrations. Thus CO 2− 3 concentrations rapidly decrease in the upper ocean layers (Fig. 8). In the deep ocean, the dissolution of calcium carbonates and the associated release of CO 2− 3 counteracts the effect of mineralization and results in low and almost constant CO 2− 3 concentrations. Figure 9 compares the resulting variations in the saturation state of oceanic waters with saturation states derived from global observations. The solubility product of calcite and thus the CO 2− 3 saturation concentrations increase with increasing pressure. As a consequence, surface waters and the shallow epicontinental ocean are generally supersaturated with respect to calcium carbonate, while deep waters are generally undersaturated. In the inner ocean, the saturation horizon that separates supersaturated and undersaturated waters is located at a waterdepth of ca. 2500 m. Global field observations reveal that the saturation horizon in the global ocean is far from uniform and shoals from 4000 m depth in the South Atlantic to 1000 m in the North Pacific. The deep polar ocean box, on the other hand, is slightly undersaturated with respect to calcite. The simulated undersaturation is partly an artifact of the box structure of the polar ocean that does not allow resolving the pressure-induced depth-variations of the solubility product.

Inorganic carbon chemistry
The inorganic carbon chemistry in the sediments ultimately controls the burial flux of carbon and thus the longterm global climate. Fig. 10 shows the benthic depth profiles of PIC (Fig. 10a), T C (Fig. 10b), T A (Fig. 10c), pH (Fig. 10d) and the saturation state of porewaters with respect to calcite (Fig. 10e). The T C and T A depth profiles show a build up with depth, which results from the combined effects of benthic organic matter mineralization and carbonate dissolution. Diffusive alkalinity and T C fluxes are thus highest in the shallow epicontinental ocean, where intense degradation rates efficiently recycle T C (Table 11). In the deep ocean, the decrease in the magnitude and quality of POC deposition fluxes and the associated low degradation rates result in lower diffusive alkalinity and T C fluxes out of the sediment. Simulated pH profiles reveal characteristic features that are often observed in sedimentary settings (e.g. Luff et al., 2000;Jourabchi et al., 2005). In response to to aerobic mineralization and secondary redox-reactions, pH decreases within the topmost sediment layers (Jourabchi et al., 2005). This decrease is followed by a relatively constant pH deeper in the sediment, where the effect of calcium carbonate dissolution exceeds the effect of decreasing organic matter mineralization and buffers pH. In the epicontinental sediments, the interstitial pH at depth is significantly lower than in the inner ocean sediments due to comparably more intense organic matter mineralization and associated secondary redox reactions in these environments. As a result, the carbonate ion concentrations rapidly decrease in the shallow sediment and the interstitial waters become undersaturated with respect to calcium carbonate (Fig. 10e). Further down, calcite dissolution together with sulfate reduction rapidly buffers the pore waters and leads again to conditions that favor calcium carbonate preservation. Due to pressure and temperature effects on CO 2− 3 concentrations and on the solubility product, inner ocean sediments are undersaturated with respect to calcite. The dissolution of calcium carbonates maintains the saturation state of pore waters with respect to calcite close to one and limits the pH drop. In contrast, the inorganic carbon system in the polar ocean is entirely controlled by organic matter decomposition since the PIC export flux is completely dissolved in the water column. Thus, pore water pH and saturation state show a regular decline with depth that is promoted by aerobic respiration.

Calcium carbonate dissolution
A fraction of the calcium carbonate that is exported from the euphotic surface layers of the ocean is dissolved in the deep ocean. The magnitude of the dissolved fraction is controlled by the saturation state of the ocean with respect to calcium carbonate and the kinetics of the dissolution reaction. Simulation results have shown that surface waters are oversaturated, while deep waters are undersaturated with re-spect to calcium carbonate (Fig. 9). The simulated PIC flux thus behaves conservatively above the saturation horizon and then gradually decrease until the average seafloor depth is reached (Fig. 4). The deep calcium carbonate dissolution controls together with organic matter mineralization the oceanic profiles of T A and T C (Fig. 7). Results indicate that 0.16 Pg C yr −1 dissolve in the deep polar and inner oceans, which is equivalent to 37 % of the pelagic export production (0.48 Pg C yr −1 ) and 22 % of the total carbonate production (0.73 Pg C yr −1 ). On a global scale, 0.32 Pg C yr −1 settle on the seafloor and 0.25 Pg C yr −1 are precipitated in form of shallow water carbonates. Because most of the PIC export reaches the seafloor, sedimentary processes play an important role in the global PIC budget. In the sediment, calcium carbonate dissolution is mainly controlled by the intensity of organic matter mineralization. Simulation results predict that another 0.11 Pg C yr −1 dissolve in sediments, a value which corresponds to 23 % of global export production and 19 % of the global carbonate deposition. Most of the dissolution flux occurs in deep ocean sediments, leaving, on a global scale, 0.27 Pg C yr −1 and 0.19 Pg C yr −1 for burial in the epicontinental and inner ocean sediments, respectively.
Published PIC export production rates reveal large discrepancies. Nevertheless, simulated fluxes and rates fall within the range of reported values. For instance, Milliman and Droxler (1996) estimated that the oceanic dissolution flux (0.5 Pg C yr −1 ) consumes half of the carbonate export production (1 Pg C yr −1 ), the remaining PIC flux sustains a benthic dissolution flux of 0.37 Pg C yr −1 and a burial flux of 0.13 Pg C yr −1 . The carbonate budget compiled by Berelson et al. (2007) proposes a global export production of 0.4-1.8 Pg C yr −1 of which 0.6 ± 0.4 Pg C yr −1 and 0.4 ± 0.3 Pg C yr −1 are dissolved in the water column and in marine sediments, respectively. Dunne et al. (2007) report a global export flux of 0.52 ± 0.15 Pg C yr −1 , which reduces to 0.29 ± 0.15 Pg C yr −1 at 2000 m water depth. Thus, in general, global estimates indicate that approximately 50 % of the calcium carbonate export production settles onto the sediment and that about 30 % of this depositional flux is ultimately buried in the sediment.

Conclusions
The new Earth system model GEOCLIM reloaded provides an improved description of the marine biogeochemical cycles of carbon, nitrogen, phosphorus and oxygen. It comprises a number of novel features that address weaknesses in previous coupled Earth system models for the distant past. GEOCLIM reloaded includes a fully coupled, verticallyresolved description of the diagenetic reaction network for redox and acid-base chemistry in marine sediments, thus affording a robust quantification of sedimentary burial and recycling fluxes. The ocean model provides an adequate representation of the general circulation with a minimum number of free parameters. Furthermore these parameters can easily be constrained on the basis of detailed GCM simulations, allowing the integration of information from cost-intensive coupled atmosphere-ocean GCMs into the computationally efficient framework of GEOCLIM reloaded. The biogeochemical model also incorporates an explicit representation of the decreasing reactivity of organic matter during settling and burial. The consistent and fully coupled description of the global biogeochemical cycles provides, without major calibration efforts, reasonable fits to global observations of concentration profiles, integrated rates and fluxes.
GEOCLIM reloaded is well suited to examine the feedbacks between climatic conditions, ocean biogeochemical dynamics and diagenesis. Its design allows exploring the sensitivity of the Earth system to a wide range of perturbations and forcings, including extra-terrestrial variations in cosmic ray intensity, changes in volcanic and mid-oceanic ridge degassing and aerosol emissions from major magmatic events, warming-induced variations of the meridional overturning circulation, fluctuations in continental weathering fluxes, or abrupt shifts in biota. The model addresses steadystate or transient regimes, thus affording the study of shortand long-term responses of the climate system to these perturbations. The proposed model will advance our understanding of the functioning of the global carbon cycle in the geological past. It will provide novel insights in the quantitative significance of the carbon burial flux in marine sediments, including its partitioning between the inorganic and organic carbon pools, the benthic recycling fluxes and their impact on the biogeochemistry of the ocean-atmosphere system, and the sedimentary carbon sink as a regulator of carbon dioxide levels and past climate. In addition, it has the potential to shed light on the geochemical conditions that led to unique depositional environments characterizing many extreme climate events. In a general sense, the newly developed Earth system model opens new perspectives for a large range of scientific questions in the broad field of palaeo-climate research.