Presentation, calibration and validation of the low-order, DCESS Earth System Model

. A new, low-order Earth System Model is described, calibrated and tested against Earth system data. The model features modules for the atmosphere, ocean, ocean sediment, land biosphere and lithosphere and has been designed to simulate global change on time scales of years to millions of years. The atmosphere module considers radiation balance, meridional transport of heat and water vapor between low-mid latitude and high latitude zones, heat and gas exchange with the ocean and sea ice and snow cover. Gases considered are carbon dioxide and methane for all three carbon isotopes, nitrous oxide and oxygen. The ocean module has 100 m vertical resolution, carbonate chemistry and prescribed circulation and mixing. Ocean biogeochemical tracers are phosphate, dissolved oxygen, dissolved inorganic carbon for all three carbon isotopes and alkalinity. Biogenic production of particulate organic matter in the ocean surface layer depends on phosphate gas inputs and other anthropogenic and natural forcing. Long term, transient model behavior is studied with a set of 100 000 year simulations, forced by a slow, 5000 Gt C input of CO 2 to the atmosphere, and with a 1.5 million year simulation, forced by a doubling of lithosphere CO 2 outgassing.

sedimentation velocities at the base of the bioturbated layer. Bioturbation rates and oxic and anoxic remineralisation rates depend on organic carbon rain rates and dissolved oxygen concentrations. The land biosphere module considers leaves, wood, litter and soil. Net primary production depends on atmospheric carbon dioxide concentration and remineralization rates in the litter and soil are related to mean atmospheric temperatures. Methane production is a small fraction of the soil remineralization. The lithosphere module considers outgassing, weathering of carbonate and silicate rocks and weathering of rocks containing old organic carbon and phosphorus. Weathering rates are related to mean atmospheric temperatures.
A pre-industrial, steady state calibration to Earth system data is carried out. Ocean observations of temperature, carbon 14, phosphate, dissolved oxygen, dissolved inorganic carbon and alkalinity constrain air-sea exchange and ocean circulation, mixing and biogeochemical parameters. Observed calcite and organic carbon distributions and inventories in the ocean sediment help constrain sediment module parameters. Carbon isotopic data and carbonate vs. silicate weathering fractions are used to estimate initial lithosphere outgassing and rock weathering rates. Model performance is tested by simulating atmospheric greenhouse gas increases, global warming and model tracer evolution for the period 1765 to 2000, as forced by prescribed anthropogenic greenhouse gas inputs and other anthropogenic and natural forcing. Long term, transient model behavior is studied with a set of 100 000 year simulations, forced by a slow, 5000 Gt C input of CO 2 to the atmosphere, and with a 1.5 million year simulation, forced by a doubling of lithosphere CO 2 outgassing.
Published by Copernicus Publications on behalf of the European Geosciences Union.

Introduction
Earth System Models are needed as tools for understanding past global changes and for projecting future global change.
The atmospheric concentration of carbon dioxide, pCO 2 , is a key feature of the Earth system via greenhouse radiative forcing of climate and this concentration is determined for different time scales by different Earth system components. For example, ocean circulation and biogeochemical cycling are most important for pCO 2 on time scales of hundreds to thousands of years whereas lithosphere outgassing and rock weathering are most important on time scales of hundreds of thousands to millions of years. Thus carbon cycle components needed in Earth System Models will depend upon time scales to be addressed in the modelling work. On the other hand, time scales that can be addressed by Earth System Models are restricted by their complexity and spatial resolution, given available computational resources. For example, comprehensive Earth System Models of high complexity and spatial resolution built around Atmosphere-Ocean General Circulation Models can nowadays be integrated over thousands of years and have been extended to include land and ocean biogeochemical cycling, important over these time scales (Cox et al., 2000;Doney et al., 2006;Schmittner et al., 2008). Only very recently have comprehensive Earth System Models of intermediate complexity and spatial resolution been developed that can be integrated over tens of thousands of years and that have been extended to include ocean sediments, important over these longer time scales (Ridgwell and Hargreaves, 2007;Brovkin et al., 2007).
Here we describe a new, comprehensive Earth System Model of low complexity and spatial resolution called the DCESS (Danish Center for Earth System Science) model. This model includes a lithosphere module with outgassing and climate dependent-rock weathering in addition to atmosphere, ocean, land biosphere and ocean sediment modules. The sediment module features a new, semi-analytical treatment and accounts for depth-and composition-dependent porosity, oxic and anoxic remineralization, calcite dissolution and burial of organic carbon and calcite. As such, this module is sufficiently flexible to deal with shallow, high productive as well as deep, low productive regions. Our model is thus suited for investigating Earth system changes on scales of years to millions of years and it can easily be integrated over these time scales.
The DCESS model has been developed in the spirit of the high-latitude exchange/interior diffusion-advection (HILDA) model in the sense that model parameters are calibrated to the greatest extent possible by fitting model output to Earth system data (Shaffer, 1993(Shaffer, , 1996Shaffer and Sarmiento, 1995). As demonstrated by the applications of the HILDA model over recent years, results from a wellcalibrated, low order model are useful and trustworthy within the bounds of their low spatial resolution (Gerber et al., 2004;Friedlingstein et al., 2006). Fast, low-order models like the DCESS model are also well suited for sensitivity studies and for hypothesis testing that can provide guidance for the application of more complex models. This paper is organized as follows: in Sect. 2, we present descriptions of the atmosphere, ocean, ocean sediment, land biosphere and lithosphere modules. In Sect. 3, we carry out a pre-industrial, steady state calibration of the full model to appropriate Earth system data. Starting from this calibration, we document model performance by simulating atmospheric greenhouse gas increases and global warming from 1765 to 2000, as forced by prescribed anthropogenic greenhouse gas inputs and other anthropogenic and natural forcing. We also present results from 100 000 year and 1.5 million year forced simulations. In Sect. 4, we supply a short discussion with conclusions. Finally, we present a more detailed description of the new, semi-analytical ocean sediment model in Appendix A.

Model description
The DCESS model contains atmosphere, ocean, ocean sediment, land biosphere and lithosphere components. Sea ice and snow cover are diagnosed from estimated meridional profiles of atmospheric temperature. The model geometry consists of one hemisphere, divided into two 360 • wide zones by 52 • latitude (Fig. 1). Values for global reservoirs, transports and fluxes are obtained by doubling the hemispheric values. The model ocean is 270 • wide and extends from the equator to 70 • latitude. In this way the ocean covers 70.5% of the Earth surface and is divided into low-mid and high latitude sectors in the proportion 84:16 as in the one hemisphere, HILDA model (Shaffer and Sarmiento, 1995). Each ocean sector is divided into 55 layers with 100 m vertical resolution to maximum depths of 5500 m (Fig. 2). Each of the 110 ocean layers is assigned an ocean sediment section. Ocean layer and sediment sector widths are determined from observed ocean depth distributions.
2.1 Atmosphere exchange, heat balance and ice/snow extent We use a simple, zone mean, energy balance model for the near surface atmospheric temperature, T a ( • C), forced by yearly-mean insolation, meridional sensible and latent heat transports and air-sea exchange. In combination with the simple sea ice and snow parameterizations, the model in-cludes the ice/snow-albedo feedback and the insulating effect of sea-ice. Prognostic equations for mean T a in the 0-52 • and 52 • -90 • zones, T l,h a , are obtained by integrating the surface energy balance over the zones. Thus, where A l,h are zone surface areas and ρ 0 C p b l,h are zone heat capacities, expressed as water equivalent capacities, whereby C p is the specific heat capacity [4×10 3 J (kg • C) −1 ], ρ 0 is the reference density of water (1×10 3 kg m −3 ), b l,h are thicknesses (b l =5 m, b h =20 m, from Olsen et al., 2005), a is the Earth's radius, θ is latitude and ξ is longitude. Furthermore, F merid is the loss (low-mid latitude) or gain (high latitude) of heat due to meridional transports across 52 • and F toa and F T are the vertical fluxes of heat through the top of the atmosphere and the ocean surface. A no flux boundary condition has been applied at the equator and at the pole.
Latitudinal variations of T a in the model are represented by a second order Legendre polynomial in sine of latitude (Wang et al., 1999), T a (θ ) = T 0 + (T 1 2) 3 sin 2 θ − 1 With T 0 and T 1 determined by matching the area-weighted, zone mean values of T a (θ ) to the prognostic mean sector values, T l,h a (θ ), in each hemisphere. The temperatures and temperature gradients entering Eqs. (3), (4), (5) and (8) below are obtained via Eq. (2).
Observations show that eddy heat fluxes in the midlatitude atmosphere are much greater than advective heat fluxes there (Oort and Peixoto, 1983). By neglecting the advective heat fluxes, Wang et al. (1999) developed suitable expressions for F merid and the associated moisture flux, E merid , in terms of T a and ∂T a /∂θ, ∂T a ∂θ m−1 ∂T a ∂θ (4) where K t is a sensible heat exchange coefficient, K q is a latent heat exchange coefficient and L v is the latent heat of condensation (2.25×10 9 J m −3 ). From observations, m is found to vary with latitude (Stone and Miller, 1980) and on this basis we take m to be 2.5 at 52 • . E merid enters (leaves) the high (low-mid) latitude, ocean surface layer. For the heat flux at the top of the atmosphere we take, where the outgoing long wave radiation is A+BT a (Budyko, 1969), whereby A and BT a are the flux at T a =0 and the deviation from this flux, respectively. This simple formulation includes implicitly the radiative effects of changes in cloud cover and in atmospheric water vapor content. Greenhouse gas forcing is modeled by taking A to depend on deviations of (prognostic) atmospheric partial pressures of the carbon dioxide, methane and nitrous oxide (pCO 2 , pCH 4 and pN 2 O; see Sect. 2.2) from their pre-industrial (PI) values such that where expressions for A CO 2 , A CH 4 and A N 2 O are taken from Myhre et al. (1998). For example, A CO 2 = 5.35 ln pCO 2 /pCO 2,PI .
We take the year 1765 as our pre-industrial baseline and pCO 2,PI , pCH 4,PI and pN 2 O PI then to be 278, 0.72 and 0.27 µatm, respectively, from ice core observations (Meure et al., 2006). Furthermore, α in Eq. (5) is the planetary albedo, equal to 0.62 for ice and snow-covered areas and to 0.3+0.0875 3 sin 2 θ−1 otherwise, including the effects of mean cloud cover and lower solar inclination at higher latitudes (Hartman, 1994). Finally, Q is the yearlymean, latitude-dependent, short-wave radiation, taken to be (Q 0 /4) 1+(Q 2 /2) 3 sin 2 θ−1 where Q 0 is the solar constant, at present 1365 W m −2 , and Q 2 is −0.482 (Hartmann, 1994). For air-sea heat exchange we take, from Haney (1971), where L o is the direct (solar) heating of the ocean surface layer, taken to be 30 and 0 W m −2 for the low-mid latitude and high latitude sectors, respectively, as a good approximation (Haney, 1971), K AS is a constant bulk transfer coefficient, taken to be 30 W m −2 • C −1 but set to zero for areas covered by sea-ice, T l,hni a are mean atmospheric temperatures for the low-mid latitude sector and for the ice-free part of the high latitude sector, and T l,h o are the zone mean, ocean surface temperatures (see below).
Finally, we take the sea ice and snow line latitudes to be located where (prescribed) atmospheric temperatures T ice and T snow are found in the atmospheric temperature profile (Eq. 2). T ice and T snow are taken to be −5 and 0 • C, respectively.

Atmosphere chemistry and air-sea gas exchange
We take the atmosphere to be well mixed for gases and consider the partial pressures of 12,13,14 CO 2 , 12,13,14 CH 4 , N 2 O and O 2 . The prognostic equation for the partial pressure of a gas χ is taken to be where υ a is the atmospheric mole volume (1.733×10 20 moles atm −1 ), A l o is the low-mid latitude, ocean surface area, A hni o is the ice free part of the high latitude, ocean surface area, l,h S are the low-mid and high latitude, air-sea gas exchange fluxes and I are sources or sinks within the atmosphere or net transports to or from the atmosphere via weathering, volcanism, interaction with the land biosphere and, for recent times, anthropogenic activities.
Air-sea exchange for 12 CO 2 is written where the gas transfer velocities k l,h w are 0.39u 2 Sc l,h /660 0.5 whereby u l,h are the mean wind speeds at 10 m above the ocean surface and Sc l,h are CO 2 Schmidt numbers that depend on prognostic temperatures of the ocean surface layers (Wanninkhof, 1992), and η l,h CO 2 are CO 2 solubilities that depend on prognostic temperatures, and to a lesser degree, on prognostic salinities of the ocean surface layers (Weiss, 1974). Also, pCO l,h 2,w are the prognostic CO 2 partial pressures of the ocean surface layers, whereby pCO l,h 2,w = [CO 2 ] l,h /η l,h CO 2 where [CO 2 ] l,h are the prognostic dissolved (or aqueous) CO 2 concentrations of the ocean surface layers as calculated from ocean carbonate carbon chemistry (see Sect. 2.4). Our atmosphere module does not calculate wind speeds so we adopt here, in both model sectors for simplicity, a global average u of 6.6 m s −1 (Archer and Jacobson, 2005). The estimate for the gas transfer velocity at 20 • C is then 17 cm h −1 . For simplicity, we often drop the superscript "12" when referring to carbon in 12 C form since about 99% of all carbon is in this form.
Air-sea exchange for i CO 2 with i=13 and 14 may be written where i α k are the kinetic fractionation factors ( 13 α k =0.99912, Zhang et al., 1995), DI i C l,h and DIC l,h are the prognostic total inorganic carbon concentrations in the ocean surface layers (see Sect. 2.4), i α l,h aw are fractionation factors due to different i CO 2 solubilities. The coefficients i α l,h wa are the overall fractionation factors due to fractionation in the dissociation reactions of ocean carbonate chemistry such that where i α l,h HCO 3 and i α l,h CO 3 are individual fractionation factors for the species HCO − 3 and CO 2− 3 and values for the ocean surface layer concentrations of these species follow from the ocean carbonate chemistry calculations (see Sect. 2.4). Fractionation factors i α l,h aw , i α l,h HCO 3 and i α l,h CO 3 all depend upon ocean surface layer temperatures, taken from Zhang et al. (1995) for 13 C. For these fractionation factors, i α k and all other fractionation factors considered below, we assume that 14 α=1−2(1− 13 α).
Air-sea exchange for O 2 may be written with k l,h w as above but with substitution for the O 2 Schmidt numbers that depend on prognostic temperatures of the ocean surface layers (Keeling et al., 1998). The O 2 solubility, η l,h O 2 , was converted from the Bunsen solubility coefficients that depend on prognostic temperatures, and to a lesser degree, on prognostic salinities of the ocean surface layers (Weiss, 1970) to model units using the ideal gas mole volume. The quantities [O 2 ] l,h are the prognostic dissolved oxygen concentrations in the ocean surface layers (Sect. 2.4). For simplicity we assume no air-sea exchange of methane species or of nitrous oxide.
The model includes the following sources/sinks of atmospheric CO 2 : net exchange with the land biosphere (Sect. 2.6), oxidation of atmospheric methane (see below), volcanic input, weathering of "old" organic carbon in rocks and weathering of carbonate and silicate rocks (see Sect. 2.7). In recent times, there have been additional, anthropogenic CO 2 sources due to fossil fuel burning and sources/sinks due to land use change (mainly deforestation). All the above sources and sinks are also included for atmospheric 13 CO 2 . For atmospheric 14 CO 2 , the above sinks and sources are included except that there are no 14 CO 2 sources from volcanoes, organic carbon weathering and fossil fuel burning, all sources of old carbon. Radiocarbon is produced in the atmosphere by cosmic ray flux and, in recent times, by atomic bomb testing. The cosmic ray source of atmospheric 14 CO 2 may be expressed as A l +A h P14 C /A vg where P14 C is the 14 C production rate (atoms m −2 s −1 ) and A vg is the Avogadro number. A small amount of 14 C enters the land biosphere and decays there radioactively. The still smaller atmospheric sink is λ C14 υ a p 14 CO 2 where λ C14 is the radioactive decay rate for 14 C (3.84×10 −12 s −1 ). By far most of the 14 C produced enters the ocean via air-sea exchange and by far most of this decays within the ocean (Sect. 2.4). A small amount of 14 C from the ocean surface layer enters the ocean sediment via the rain of biogenic particles and part of this returns to the ocean after remineralization/dissolution in the sediment (see Appendix A).
Methane is produced by the land biosphere (Sect. 2.6) and, in recent centuries, by human activities. "Melting" of methane hydrate in the arctic tundra and in ocean sediments may represent yet another methane source. The main atmospheric sink of methane is associated with reaction with OH radicals in the troposphere. Since this reaction depletes the concentration of these radicals, atmospheric lifetimes for methane grow as methane concentrations increase. We model this effect, together with the effects of associated chemical reactions in the troposphere and stratosphere, by fitting a simple model to results from a complex atmospheric chemistry model (Schmidt and Shindell, 2003). Thus, we take the atmospheric methane sink to be λ CH 4 pCH 4 where where M≡ pCH 4 −pCH 4,PI /pCH 4,PI , s y is the number of seconds in a year and L CH 4 , L CH 4 ,PI are atmospheric lifetimes of methane in years. We found a good fit to the results of Schmidt and Shindel (loss rates in their Table 1) for a=0.96 and b=6.6. The atmospheric methane concentration and the total natural + anthropogenic methane sinks have been about 1.77 µatm and 0.581 Gt (CH 4 ) yr −1 in recent years ( Atmospheric sinks of O 2 are associated with weathering of organic carbon in rocks and oxidation of reduced carbon emitted in lithosphere outgassing (Sect. 2.7). In the model, a long term, steady state of pO 2 is achieved when these sinks balance the O 2 source associated with burial of organic matter in ocean sediments (Sect. 2.5). This O 2 source leads to net, long term transport of O 2 from the ocean to the atmosphere via air-sea exchange. Additional atmospheric sinks (sources) of O 2 are associated with decreasing (increasing) biomass on land and, in recent times, with burning of fossil fuels.
Finally, we consider atmospheric cycling of the oxygen isotopes of water. Fractionation during evaporation enriches the low-mid latitude surface ocean and depletes low-mid latitude atmospheric moisture in 18 O. Atmospheric moisture is further depleted via condensation upon poleward transport and associated cooling. Upon precipitation, this moisture leaves the high-latitude ocean depleted in 18 O. Here we use the approach described in detail in Olsen et al. (2005) to model these processes, making use of the atmospheric temperature at the latitude dividing the two model zones (52 • ), as calculated from the meridional temperature profile (Eq. 2).
A key result of these calculations is an estimate for the atmospheric content of 18 O at the dividing latitude, δ 18 O a (θ=52 • ) where δ 18 O is defined in the conventional way relative to Standard Mean Ocean Water (SMOW).

Ocean circulation and mixing
The HILDA model serves as the point of departure for the DCESS model formulation of ocean physics and biogeochemical cycling. As in HILDA, four physical parameters characterize ocean circulation and mixing. For the DCESS model with continuous vertical stratification in both sectors these four parameters are 1) a transport, V , associated with high latitude sinking and the deepest, low-mid latitude upwelling, 2) constant horizontal diffusion between the zones, K h , associated with the wind-driven circulation and deep recirculation (Shaffer and Sarmiento, 1995), 3) strong, constant vertical diffusion in the high latitude zone, K h v , associated with high latitude convection and, 4) a weak, depthdependent, vertical diffusion in the low-mid latitude zone, K l v (z). A discussion of real ocean physical analogues of related parameters in the HILDA model is given in Shaffer and Sarmiento (1995).
The deep overturning circulation, V , equals the poleward flow in the model surface layer. The transport down out of the high latitude mixed layer, equatorward between the deepest model layers at 5500m and upward into the low-mid latitude surface layer is V +E merid , the sum of ocean surface layer and atmosphere water transports from the low-mid latitude to the high latitude model zone (note that V E merid ).
The water transported in the atmosphere contains no dissolved substances. This approach leads to realistic forcing of surface layer concentration/dilution of dissolved substances (like salt and alkalinity) and avoids the use of artificial salt fluxes and the need for dealing with salinity-normalized, dissolved substances. The depth-dependent vertical velocity for each zone, w l,h (z), is calculated from continuity using model bathymetry such that w l, are observed low-mid and high latitude zone ocean areas as functions of ocean depth (Fig. 2).
In the low-mid latitude zone, vertical diffusion is calculated as where K l v,0 is a vertical diffusion scale, N l obs (z) is the observed mean Brunt-Väisäla frequency profile for the low-mid latitude zone and K l v,0 , N l obs,0 are corresponding values at 100 m depth. N l obs (z) is equal to g ρ 0 ∂ρ l obs (z) ∂z 0.5 where g is gravity, ρ 0 is mean water density and the observed water density profile, ρ l obs (z), has been calculated from observed mean profiles of temperature and salinity from this zone (Fig. 3a, b). The K l v (z) parameterization is consistent with diapycnal mixing via breaking of internal waves (Gargett and Holloway, 1984). We found a good fit to the results of this calculation with which describes a five fold increase in the vertical diffusion from the surface to the bottom layer, similar to the Bryan and Lewis (1979) vertical diffusion profile often assumed in Ocean General Circulation Models. With the above physics and for each of the two ocean zones, conservation of an ocean tracer ϕ may be written where S denote air-sea exchange of heat, fresh water and gases, B denote exchange of dissolved substances with the ocean sediment (see also Sect. 2.5 and Appendix A) and I denote internal source/sink terms. For the surface layers, the internal source/sink terms include biological production, river inflow of dissolved substances, and the direct solar heating (low-mid latitude only). For the ocean interior, these terms include remineralization and dissolution of biogenic matter produced in the surface layers. From the assumptions above, the meridional velocity is only defined in the surface and the deepest model layers and is set by V and V +E merid , respectively.
The tracers temperature, salinity and 18 O content of water are forced at or near the ocean surface only. Ocean temperature, T l,h w (z), is forced by air-sea heat exchange and by direct solar forcing of the low-mid latitude surface layer (Eq. 8). Concentration/dilution of the surface layer as described above provides the forcing for ocean salinity, S l,h (z). Ocean mean salinity is 34.72 from observations. Ocean distributions of 18 O content in water, δ 18 O l,h w (z), are forced by sources/sinks in the mid-low/high latitude surface layers equal to ±E merid δ 18 O a (θ=52 • ). Ocean mean δ 18 O w is taken to be zero. We also calculate the distribution of the 18 O content in biogenic carbonate, δ 18 O l,h c (z), whereby for benthic carbonate deposits, (Bemis et al., 1998). The same relation is used for pelagic carbonate but with the use of surface layer δ 18 O w and T .

Ocean biogeochemical cycling
The biogeochemical ocean tracers considered here, phosphate (PO 4 ), dissolved oxygen (O 2 ), dissolved inorganic carbon (DIC) in 12,13,14 C species, and alkalinity (ALK), are all forced by net (new) production of organic matter and CaCO 3 shells in the lighted surface layers. In addition, O 2 and DI 12,13,14 C are forced by air-sea exchange and PO 4 , DI 12,13,14 C and ALK are forced by river inputs to the surface layer and concentration/dilution of this layer by evaporation/precipitation. In the subsurface layers, the biogeochemical ocean tracers are forced by remineralization of organic matter and dissolution of CaCO 3 shells in the water column as well as exchange with the ocean sediment. DI 14 C is affected by radioactive decay in all ocean layers. For simplicity, we have neglected explicit nitrogen cycling, i.e. phosphate is assumed to be the basic limiting nutrient, and have assumed that all biogenic export from the surface layer is in the form of particles and that all CaCO 3 is in the form of calcite. New production of organic matter in the surface layer is parameterized in terms of phosphorus (mol P m −2 s −1 ) as  -Reimer, 1993;Yamanaka and Tajika, 1996) where A l,hni o are the ice-free ocean surface areas, z eu is the surface layer depth (100 m), sy is the number of seconds per year, PO l,h 4 are the phosphate concentrations in the surface layer, and P 1/2 is a half saturation constant (1 µmol/m 3 ). L l,h f are efficiency factors, taken to be 1 for the low-mid zone and to some lower value for the high latitude zone, as determined by model fit to ocean data (Sect. 3). This is the way that the model accounts for light and iron limitation in this zone. For simplicity, we neglect dissolved organic matter such that the rate of export of particulate organic matter (POM) down out of the surface layer is equal to new production.
Sources/sinks in the surface layer due to new production are −NP l,h for PO 4 , −r CP NP l,h for DI 12 C, (r OCP +r ONP )NP l,h for O 2 , and r AlkP NP l,h for ALK, where r CP , r OCP , r ONP and r AlkP are (Redfield)  . Steady state, pre-industrial simulations (solid lines) compared to data (dots) of mean ocean profiles of (a) temperature (T ), (b) salinity (S) and (c) 14 C. Low-mid latitude and high latitude simulations are given in red and blue, respectively. Mean, databased profiles from the low-mid and high latitude sectors are given in black and grey, respectively. These profiles have been calculated from GEOSECS data as in Shaffer and Sarmiento (1995) and Shaffer (1996).
(O 2 ) N :P and ALK:P mole ratios, respectively. The subscripts C and N refer to a division of POM produced into "carbon" and "nutrient" parts, respectively, as explained below. Values for these stoichiometric ratios are taken to be 106, 118, 32 and 16 in the above order, whereby an implicit N:P ratio of 16 in new production has been assumed and the enhanced value of r OCP above the traditional Redfield et al. (1963) value of 106 reflects the influence of phytoplankton protein and lipids (Anderson, 1995), not considered in the original Redfield stoichiometry. The above is consistent G. Shaffer et al.: DCESS Earth System Model with an -O 2 :C assimilation ratio in new production of about 1.4 (Laws, 1991). For DI 13,14 C, the surface layer sinks due to new production and associated isotope fractionation are − i α l,h Org (DI i C/DI 12 C) l,h eu r CP NP l,h (i=13 and 14) where the subscript eu refers to euphotic zone (surface layer) values. We use 13 α l,h Org =1−{17 log (CO 2 (aq)) l,h eu +3.4}/1000 (Popp et al., 1989). The empirical relationship for 13 α Org assumes that (aqueous) CO 2 concentrations (in units mmol m −3 ) mainly control this fractionation during primary production in the ocean. For a warm to cold ocean CO 2 (aq) range from 7 to 24 mmol m −3 , this leads to a ‰ fractionation from about −18 to −27‰ in the organic carbon produced.
The surface layer production of biogenic calcite carbon is expressed as r l,h CalC r CP NP l,h where r l,h CalC are mole ("rain") ratios between new production of organic and calcite 12 C in the surface layer. The rain ratio is parameterized as (Maier-Reimer, 1993) where r CalC,m is a rain ratio upper limit, µ is a steepness factor and T ref is a reference temperature, taken to 10 • C. Both r CalC,m and µ are determined by model fit to ocean and ocean sediment data (see Sects. 3.1.3 and 3.2.1). Equation (20) yields lower rain ratios for lower temperatures as indicated by observations (Tsunogai and Noriki, 1991). Surface layer DI 12 C and ALK sinks from biogenic calcite production are −r l,h CalC r CP NP l,h and −2r l,h CalC r CP NP l,h , respectively. Surface layer DI 13,14 C sinks from this production and associated fractionation are − i α Cal (DI i C/DI 12 C) l,h eu r l,h CalC r CP NP l,h where i=13 and 14. Here 13 α Cal is taken to be 0.9988, corresponding to a fractionation of −1.2‰.
Particles are assumed to sink out of the surface layer with settling speeds high enough to neglect advection and diffusion of them and to take subsurface remineralization or dissolution of them to occur instantaneously. Following Shaffer (1996) and further motivated by the results of Shaffer et al. (1999), we assume an exponential law for the vertical distribution of remineralization of POM "carbon" and "nutrient" components, each with a distinct e-folding length, λ N and λ C , respectively. Likewise, we assume an exponential law for the dissolution of biogenic calcite particles with an e-folding length, λ Cal .
The vertical distributions of DI 13,14 C sources from remineralization and dissolution are where i=13 and 14. For simplicity, low-mid and high latitude values for λ N , λ C and λ Cal are assumed to be the same. Total sources/sinks of PO 4 , DI 12,13,14 C, O 2 and ALK from remineralization/dissolution at any depth z of each zone are calculated as the product of A l,h O (z) and remineralization/dissolution fluxes at z of each zone, as calculated from Eqs. (21) and (22). The fluxes of P and 12,13,14 C that fall in the form of POM and/or biogenic calcite particles on the model ocean sediment surface at any depth z of each zone are calculated as the product of dA l,h O (z)/dz there and the difference between the particulate fluxes falling out of the ocean surface layer and the remineralization/dissolution taking place down to z of each zone, as calculated by integrating Eqs. (21) and (22). Note that for λ N =λ C , C:P mole ratios of POM falling on the sediment surface vary with water depth.
Model calculations of air-sea exchange of carbon dioxide, carbon isotope fractionation during air-sea exchange and in ocean new production and dissolution of calcite in the ocean sediment require information on ocean distributions of CO 2 (aq) or CO 2− 3 . These distributions are calculated from the ocean carbonate chemistry equations, given pressure, model distributions of DIC, ALK, T and S and appropriate apparent dissociation constants for carbonic acid, boric acid and sea water as functions of T , S and pressure. Equations for these constants, and a relation between total borate concentration and S are from Millero (1995;his Eqs. 35, 36, 52, 53 and 62) but with corrections listed in http://cdiac.ornl.gov/oceans/co2rprt.html#preseff). Alkalinity includes hydroxide and hydrogen ion concentration but not minor bases. This nonlinear system is solved for all carbon species with the recursive formulation of Antoine and Morel (1995). The calculations also yield distributions of hydrogen ion concentration (including H + bound to SO 2− 4 and F − ) from which pH (seawater scale) is calculated. Profiles of CO 2− 3 saturation with respect to calcite are calculated as K CaCO 3 (z)/{(Ca 2+ ) 0 S(z)/35}, where K CaCO 3 is the apparent dissociation constant for calcite as a function of T , S and pressure (Mucci, 1983) and (Ca 2+ ) 0 is mean ocean Ca 2+ , 10.57 mol m −3 for present day.

Ocean sediment
For the DCESS model we have developed a new, fast, semianalytical module for addressing calcium carbonate dissolution and (oxic and anoxic) organic matter remineralization in ocean sediments. Module details are given in Appendix A. There is a sediment segment for each of the 2×55 ocean model layers. The area covered by each segment is determined by model topography (Fig. 2). Each segment is characterized by a bioturbated layer (BL), assumed to be 0.1 m Geosci. Model Dev., 1, 17-51, 2008 www.geosci-model-dev.net/1/17/2008/ thick. We do not consider dissolution and remineralization below the BL. The sediment is composed of calcite, CaCO 3 , non-calcite mineral, NCM, and reactive organic matter fractions. To a good approximation, CaCO 3 and NCM fractions are taken to be well mixed in the BL by bioturbation, D b , but the reactive organic matter fraction varies over the BL due to relatively rapid remineralization. We assume D b to be constant over the BL and to depend on organic carbon rain rates and ambient dissolved oxygen concentrations. For each sediment segment, NCM fluxes are prescribed based on data (see Eq. A30 in Appendix A) and calculated POM and calcite rain fluxes, O 2 , DIC, ALK, T , S and pressure are accepted from corresponding layers of the ocean module. Calculated fluxes of PO 4 , O 2 , DIC and ALK to/from the sediment are then removed from/returned to these ocean layers. As part of these calculations, the calcite dry weight/total sediment dry weight ratio, (CaCO 3 ) dwf , and the sedimentation velocity, w s , out of the BL are carried as prognostic variables. With this approach we can also deal with transient upward w s due to strong CaCO 3 dissolution, and the associated entrainment of buried CaCO 3 back up into the BL.
For each sediment segment, solutions are found for steady state profiles of reactive organic carbon (OrgC), pore-water O 2 and pore-water CO 2− 3 . From the latter, dissolution of CaCO 3 in the BL is calculated and used in steady state and time dependent calculations of calcite dry weight fractions. The OrgC and O 2 solutions are coupled and are solved in a semi-analytical, iterative approach, given (CaCO 3 ) dwf and w s . Explicit solutions are not sought for other species produced/consumed during anoxic respiration but the influence of these species is included via boundary conditions on O 2 . The CaCO 3 and CO 2− 3 solutions are coupled and are solved with a semi-analytical, iterative approach (steady state) or a semi-analytical, time stepping approach (time dependent). In principle, these latter solutions also depend upon the OrgC/O 2 solutions via the release of CO 2 during remineralization and carbonate chemistry (Archer, 1991). Here we do not treat this effect explicitly, as it needs careful treatment for the general oxic/anoxic case we consider (see Appendix A). However, we do include this effect implicitly to some degree in the "water column" dissolution of CaCO 3 (Sect. 2.4). Steady state solutions for pore-water concentrations hold on time scales of years while steady state solutions for reactive carbon hold on times scales increasing from tens of years on continental shelves to hundreds of years in the deep ocean.

Land biosphere
We consider a land biosphere model with carbon isotope reservoirs for leaves (M G ), wood (M W ), litter (M D ) and soil (M S ). Net primary production on land, NPP, takes up CO 2 from the atmosphere and depends on pCO 2 according to where NPP PI is the pre-industrial level of NPP and F CO 2 is a CO 2 fertilization factor. Following Siegenthaler and Oeschger (1987), we take NPP PI to be 60 Gt C yr −1 and the pre-industrial M G,W,D,S reservoirs to be 100, 500, 120 and 1500 Gt C, respectively. In a comparison of 11 coupledcarbon cycle models, Friedlingstein et al. (2006) found relatively strong CO 2 fertilization in all but one of the models (their Fig. 3a). We found a good agreement to an average of their results by using Eq. (23) with F CO 2 equal to 0.65, a value adopted below. For example, this leads to an increase of NPP from 60 to 87 Gt C yr −1 for a pCO 2 doubling from 278 to 556 µatm. Friedlingstein et al. (2006) found no model consensus on temperature dependence of NPP and very little dependence in five of the models of the intercomparison. On this basis we have neglected any such dependence in Eq. (23). Following Siegenthaler and Oeschger (1987), NPP is distributed between leaves and wood in the fixed ratio 35:25, all leaf loss goes to litter, wood loss is divided between litter and soil in the fixed ratio 20:5, litter loss is divided between the atmosphere (as CO 2 ) and the soil in the fixed ratio 45:10. Soil loss is to the atmosphere as CO 2 and, to a lesser extent, as CH 4 (see below). Organic burial on land is not considered.
Losses from all land reservoirs are taken to be proportional to reservoir size and, for litter and soil, to also depend upon mean atmospheric temperature according to where Q 10 is a (biotic) activity increase for a 10 degree increase ofT a . We chose a value for Q 10 of 2, a typical choice in carbon-climate models (Friedlingstein et al., 2006). For temperatures at or aboveT a,PI , simple relationships of this type approximate well the results of complex global vegetation models (Gerber et al., 2004).
Methane and nitrous oxide production occur in the soil, are proportional to the soil reservoir size and depend upon T a according to λ Q , again with a Q 10 of 2. A preindustrial balance is assumed between soil production and atmospheric consumption of CH 4 and N 2 O. Our assumptions on the climate dependence of methane and nitrous oxide production are simple but lead to results that are consistent with observed glacial-interglacial changes. A Last Glacial Maximum methane consumption of 0.150 Gt (CH 4 ) yr −1 follows from Eq. (14) and an LGM atmospheric content of 0.35 µatm (Jouzel et al., 1993). We find a matching methane production for an atmospheric temperature decrease of 5 • C (Scheider von Deimling et al., 2006) and a land biosphere carbon stock decrease of 300 Gt C from LGM reconstructions (Francois et al., 1999). A comparable calculation for nitrous oxide with an LGM atmospheric content of 0.185 µatm yields a consumption-production match for a Q 10 of 1.6 for the same temperature and biosphere carbon stock decreases as above (Leuenberger and Siegenthaler, 1992). For simplicity, we retain a Q 10 of 2 for this minor greenhouse gas. Inputs of 12,13,14 C to the atmospheric CH 4 pool from the soil are λ Q (λ CH 4 ,PI pCH 4,PI )(M S /M S,PI ) and 13,14 α M λ Q (λ CH 4 ,PI pCH 4,PI )( 13,14 M S /M S,PI ) where 13 α m is the 13 C fractionation factor for CH 4 production, taken to be 0.970 corresponding to a −30‰ fractionation and a δ 13 C value of about −55‰ for CH 4 released from the soil (Quay et al., 1988).
The conservation equations for the land biosphere reservoirs of 12 C are The conservation equations for the land biosphere reservoirs of 13 C are where 13 α L is the 13 C fractionation factor for land photosynthesis, taken to be 0.9819. This corresponds to a −18.1‰ fractionation, reflecting domination of C3 over C4 plant productivity, and a land biosphere δ 13 C value of about −25‰ (Joos and Bruno, 1998). Conservation equations for the land 14 C reservoirs are similar to Eqs. (28-31) but with additional radioactive sinks, −λ C14 14 M G,W,D,S . Inputs of 12,13,14 C to the atmospheric CO 2 pool from the land biosphere are .
With the above parameter and reservoir size choices, the preindustrial, steady state solutions for land biosphere 12,13,14 C are fully determined by prescribed pCO 2,PI and pCH 4,PI .

Rock weathering, volcanism and river input
Climate-dependent weathering of rocks containing phosphorus, W P , is taken to supply dissolved phosphorus for river input, R P , such that We assume that 80% and 20% of R P enter the low-mid and high latitude ocean surface layer, respectively, as shown by river runoff observations (Dai and Trenberth, 2002). For simplicity, we use the same Q 10 -based, climate dependency for weathering as for other model components above ) and again we take Q 10 to be 2. This gives a weathering dependence on global temperature very similar that from the function e (T a −T a,PI )/13.7 used in earlier work (Volk, 1987). Weathering rates may also depend on other factors like continental surface area and types of exposed bedrock (Munhoven, 2002). These could be included into the model if needed for a specific application.
Changes in the total phosphorus content of the ocean reflect imbalances of net inputs and outputs or where M P is the ocean phosphorus content and B OrgP is the total burial rate of phosphorus in organic matter, OrgP (i), from the ocean sediment module for each of the n bottom segments with separation in depth z for each of the ocean sectors. For simplicity, we assume a pre-industrial balance between total weathering and ocean burial or W P ,PI −B OrgP,PI =0. Multiplication of the above relationship for phosphate by r AlkP yields corresponding relationships for the part of the ALK fluxes associated with organic matter cycling.
The overall carbon balance of the model and the distribution of carbon among the different model components are influenced by climate-dependent weathering of carbonate and silicate rocks (W Cal and W Sil ), climate-dependent weathering of rocks containing old organic carbon (W OrgC ), and lithosphere outgassing (Vol). In simple terms, silicate weathering may be described by the left to right reactions in two reaction steps www.geosci-model-dev.net/1/17/2008/ and Thus there is an atmosphere sink of 2 moles of CO 2 per mole of silicate mineral weathered. Carbonate weathering is described by the left to right reaction in the second reaction step with an atmosphere sink of 1 mole of CO 2 per mole of carbonate mineral weathered. Both weathering types supply HCO − 3 for river input. Then river inputs of DIC and ALK are and where we assume the same 80%-20% river input partition as above and the same Q 10 -based, climate dependency with Q 10 =2 (note that we do not consider river input of organic carbon). In the assumed, pre-industrial steady state, just enough of the biogenic carbonate falling on the sediment surface is buried to satisfy where total, pre-industrial calcite burial rate, B Cal,PI , is calculated as for phosphorus burial above but with calcite burial fluxes. Likewise, pre-industrial CO 2 outgassing from the ocean to the atmosphere is equal to W Sil,PI +W Cal,PI . An analogous relation for river input of 13 C can be derived in a similar way, with proper account taken for the fact that all and half of the carbon involved in W Sil and W Cal , respectively, stem from atmospheric CO 2 , and with a proper choice for the 13 C content of carbonate rock weathered during W Cal (see below). W OrgC and Vol are the two external sources of atmospheric CO 2 in the model. We again adopt the same Q 10 -based, climate dependency and Q 10 =2 for weathering of old organic carbon such that W OrgC =λ Q W OrgC,PI . The sources of Vol are thermal breakdown of buried carbonate and organic carbon. Vol may either be taken constant and equal to its preindustrial value, Vol PI (see below) or may be prescribed as external forcing of the Earth system. Therefore, the total source of atmospheric CO 2 from lithospheric processes is W OrgC +Vol−2W Sil −W Cal (all in units of moles of carbon per unit time).
Changes in total carbon content of the combined atmosphere-ocean-land biosphere-ocean sediment system reflect imbalances of net inputs and outputs or where M C is the total carbon content of this system and B OrgC is the total organic carbon burial rate, calculated as for phosphorus burial above but with organic carbon burial fluxes. For the assumed, pre-industrial steady state, we have from Eqs. (36) and (37), With the additional assumption that silicate weathering takes place at a fixed ratio, γ Sil , to carbonate weathering, Eq. (36) gives W Sil,PI =γ Sil (1+γ Sil ) −1 B Cal,PI and Eq. (38) becomes From a detailed study of global weathering sources, we take γ Sil =0.85 (Lerman et al., 2007). An analogous, steady state equation for 13 C may be written as where the δ 13 C are defined in the usual way relative to the PDB standard, δ 13 C Vol,PI is taken to be −5‰ (Kump and Arthur, 1999) and δ 13 C Cal,PI and δ 13 C OrgC,PI are calculated from the steady state, pre-industrial model results (Sect. 3.2.3) as sector burial-weighted means of 13 C contents in organic carbon and biogenic carbonate produced in each ocean sector. Thus, OrgC,PI +B h OrgC,PI ) −1 . For simplicity, we have assumed that the 13 C content of old and new organic matter are the same in this steady state. We can now solve for W OrgC,PI and Vol PI using Eqs. (39) and (40): Changes in total oxygen content of the atmosphere and ocean reflect imbalances of net inputs and outputs whereby where M O is the total oxygen content of the atmosphere and ocean, M L is the total carbon content of the land biosphere (M G +M W +M D +M S ), r ONP , r OCP and r CP are defined in Sect. 2.4, and r OC,L is the -O 2 :C assimilation ratio in land new production, taken to be 1.1. The factor multiplying Vol gives the fraction of Vol from old organic matter. In Eq. (43), flows of phosphorus "stand in" for flows of nitrogen, not treated specifically here.

G. Shaffer et al.: DCESS Earth System Model
For the assumed, pre-industrial steady state, Eq. (43) becomes

Model solution, calibration and validation
3.1 Pre-industrial steady state solution for the atmosphere and ocean 3.1.1 Solution procedure Ocean tracer equations are discretized on a grid defined by the two meridional zones, with vertically-decreasing horizontal extents as set by observed ocean topography, and by constant vertical resolution of 100 m. Each of these 2×55 ocean boxes is associated with an ocean sediment segment with horizontal extent set by the observed topography. Model tracer boundary conditions account for the two-way exchange of heat and freshwater between the atmosphere and the ocean as well as the two-way exchange of gases between the atmosphere and the ocean, the land biosphere and the lithosphere. Other tracer boundary conditions account for particulate matter fluxes from the ocean to the ocean sediment, two-way exchange of dissolved substances between the ocean and the ocean sediment, and fluxes of dissolved and particulate matter from the lithosphere to the ocean. Prognostic equations for the atmosphere (including snow and sea ice cover), the land biosphere, the lithosphere and the ocean are solved simultaneously using a fourth order Runge Kutta algorithm with a two week time step. Prognostic equations for the ocean sediment are solved by simple time stepping with a one year time step. The complete, coupled model was written in Matlab and runs at a speed of about 10 kyr per hour on a contemporary PC.

Calibration procedure
We calibrated the parameters of the DCESS model by "trial and error" in four steps. In the first calibration step we considered the atmosphere module only with atmospheric pCH 4 and pN 2 O set at pre-industrial values of 0.72 µatm and 0.27 µatm, respectively, and with atmospheric pCO 2 set at its pre-industrial value of 278 µatm or at twice that value, 556 µatm (Etheridge et al., 1998a, b;Meure et al., 2006). We adjusted the four free parameters of this module (Table 1) to give a global mean atmospheric temperature of 15 • C and a climate sensitivity of 3 • C per doubling of atmospheric pCO 2 and poleward heat and water vapor transports in the atmosphere across the sector boundary consistent with observations (Trenberth and Caron, 2001;Dai and Trenberth, 2002).
In the second calibration step, we considered the ocean and land biosphere modules coupled to the preliminarily-calibrated, atmosphere module, with pre-industrial atmospheric δ 13 C of −6.4‰, observed atmospheric pO 2 of 0.2095 atm and observed ocean mean PO 4 , DIC and ALK of 2.12×10 −3 , 2.32 and 2.44 mol m −3 (Francey, 1999;Keeling et al., 1998;Shaffer, 1993Shaffer, , 1996. Atmospheric 14 C production was adjusted to maintain atmospheric 14 C at 0‰. In this step we assumed that all biogenic particles falling to the ocean bottom remineralize completely there. We then made initial guesses for the values of the 10 free parameters of the ocean module of which 4 are physical and 6 are biogeochemical (Table 1). These guesses were based in part on results from the HILDA model calibration (Shaffer, 1993(Shaffer, , 1996Shaffer and Sarmiento, 1995) and in part on literature values (Maier-Reimer, 1993). Note that values of all parameters in the land biosphere model were determined a priori.
The resulting atmosphere-ocean-land biosphere model was spun up from uniform atmosphere and ocean tracer distributions to a steady state solution after about 10 000 model years and results were compared with atmosphere and ocean data. Then parameter values were adjusted by expertiseguided "trial and error" to obtain steady state solutions that better satisfied a global mean atmospheric temperature of 15 • C, a climate sensitivity of 3 • C per doubling of atmospheric pCO 2 , an atmospheric pCO 2 value of 278 µatm, observed poleward heat and water vapor transports across the sector boundary (52 • latitude) and observed ocean distributions of T, 14 C, PO 4 , O 2 , DIC and ALK (Shaffer, 1993(Shaffer, , 1996Shaffer and Sarmiento, 1995). Note that atmospheric transport parameters were reduced in this step to obtain good model fit in the presence of meridional heat transports in the ocean.
In the third calibration step, we coupled our ocean sediment module to the model of the second step. In the resulting closed system, a tracer flux was added in appropriate form to the ocean surface layer of a sector at a rate equal to the total burial rate of that tracer in that sector from the sediment module. We then made initial guesses for the values of the 4 free parameters of the sediment module (Table 1), based on published results (Archer et al., , 2002. With these values and parameter values from the second calibration step, we solved for a new, pre-industrial steady state, taking advantage of the fast, steady state mode of the sediment module (see Appendix A). Then all model parameter values were adjusted by expertise-guided "trial and error" to obtain steady state solutions that better satisfied the data-based constraints of the second calibration step in addition to observed calcite and organic carbon distributions and inventories in the ocean sediment (Hedges and Keil, 1995;Archer, 1996a).
In the fourth and final calibration step, we coupled our lithosphere module to the model of the third step whereby river inputs of tracers were equated with tracer burial fluxes from the third step and tracer burial fluxes in this new open model leave the system. As outlined in Sect. 2.7, we also calculated weathering rates and lithosphere outgassing from the tracer burial fluxes and the convenient assumption of a  pre-anthropogenic steady state for P and 12,13 C. We then made slight final "trial and error" adjustments of model parameters values and of ocean mean PO 4 , DIC and ALK until the model steady state results again broadly satisfied the data-based constraints of the third calibration step. For this final calibration, atmospheric production of 14 C is 1.8752×10 4 atom m −2 s −1 and the ocean mean PO 4 , DIC and ALK are 2.089×10 −3 , 2.318 and 2.434 mol m −3 , respectively, that may be compared with the observed ocean inventories above.

Ocean tracer and biological production results
Our tuned parameter values for this pre-industrial, steady state calibration are listed in Table 1. The sea ice and snow lines for this solution are found at 63.5 and 55.8 • latitude, respectively. The total poleward heat transport across 52 • latitude in this steady state is 5.0 PW, with ocean and atmosphere contributions of 0.7 PW and 4.3 PW, respectively. Poleward water vapor transport in the atmosphere there is 0.36 Sv (1 Sv=10 6 m 3 s −1 ). All these transport estimates agree well with recent data and reanalysis based estimates (Trenberth and Caron, 2001;Dai and Trenberth, 2002). The atmospheric heat transport is divided between sensible and latent heat as 3.44 PW and 0.80 PW, respectively. The ocean heat transport consists of 0.23 PW in the deep upwelling circulation, V , and 0.50 PW in the wind-driven circulation and deep recirculation associated with K h . Equation (16) and the best fit estimate for the diffusion scale lead to an increase in vertical diffusion from 2×10 −5 m 2 s −1 to 10.2×10 −5 m 2 s −1 down through the lowmid latitude ocean. This agrees with observations of weak background mixing combined with bottom intensified mixing near rough topography (Ledwell et al., 1998;Polzin et al., 1997). However, our simple model does not capture the vertical component of ocean isopycnal mixing, an important component of upper ocean vertical exchange of tracers like 14 C and O 2 (Siegenthaler and Joos, 1992). Our simultaneous tuning to fit to such tracers and to temperature (which largely defines isopycnals and is therefore not mixed along them) is a tradeoff tending to overestimate the effective vertical exchange of heat but to underestimate the effective vertical exchange of the other tracers. But this is a useful tradeoff since it helps limit the number of model free parameters while still allowing good model agreement with observations.
Model ocean profiles of T , S and 14 C are shown in Fig. 3, together with data-based, sector mean profiles of these tracers (all mean tracer profiles in Figs. 3, 4, 5 and 7 have been calculated from GEOSECS data as in Shaffer and Sarmiento, 1995;Shaffer, 1996). In the 14 C comparison, . Low-mid latitude and high latitude simulations are given in red and blue, respectively. Mean, data-based profiles from the low-mid and high latitude sectors are given in black and grey, respectively. These profiles have been calculated from GEOSECS data as in Shaffer and Sarmiento (1995) and Shaffer (1996).
only ocean data from below 1000 m have been used as shallower depths are strongly affected by atomic bomb 14 C. With best fit parameters, the model achieves generally good fits to T and 14 C data. The high latitude temperature observations reflect deep water formation in geographically restricted sites not resolvable in our simple model. Model 14 C minimum for the low-mid latitude sector lies about 1 km deeper than in the observations and model 14 C values for the high latitude sector are a little high in the depth range 1000-2500 m. Model fit to the salinity data is not very good since salinity distributions in the real ocean are strongly controlled by vertically-structured, advective processes not captured in our simple model. In particular, the salinity minimum in the data at about 800 m depth reflects low saline, intermediate waters formed in the 50-60 • latitude band. The presence of these waters also helps maintain low-mid latitude, surface layer salinity relatively low.
Model ocean profiles of PO 4 , O 2 , DIC and ALK are shown in Fig. 4, together with data-based, mean sector profiles of these tracers. With best fit parameters, the model achieves good fits to PO 4 , O 2 , and DIC data in the low-mid latitude sector. High latitude sector differences in vertical structure between data and model simulations reflect geographically restricted deep water formation and vertically-structured, advective processes, as mentioned above. As in Shaffer (1996), simultaneous tuning to fit PO 4 and O 2 data reveal slower remineralization of the "carbon" component compared to the "nutrient" component of POM, as reflected by λ C >λ N in Table 1. This important property of POM remineralization in the ocean was also documented in an in-depth analysis of ocean tracer data (Shaffer et al., 1999). Model fit to ocean ALK data in the low-mid latitude sector is less impressive but still serves to help constrain the global biogenic calcite production and the calcite dissolution length scale, λ Cal . As for salinity, the relative model misfit to low-mid latitude ALK data can be traced to a relatively strong influence of vertically-structured, advective processes in the ocean. Model ocean profiles of CO 2− 3 are shown in Fig. 7a. The crossing point for CO 2− 3 and the CO 2− 3 saturation profiles is the calcite saturation depth (CSD). Model CSD's are 2928 m and 3186 m for the low-mid and high latitude zones, respectively. The model CSD's are 400-500 m shallower than the data-based estimates (Fig. 7a). This can be traced back to the ALK profile misfits discussed above. Low-mid latitude and high latitude simulations are given in red and blue, respectively. Mean, data-based profiles from the low-mid and high latitude sectors are given in black and grey, respectively. These profiles have been calculated from GEOSECS data as in Shaffer and Sarmiento (1995) and Shaffer (1996). Also shown are (a) simulated biogenic calcite δ 18 O (dashed lines), calculated from simulated ocean T and δ 18 O w using Eq. (18) and (b) simulated δ 13 C (dashed lines) for an alternative temperaturedependent fractionation (see Sect. 3.1.3). Note that ocean uptake of fossil fuel CO 2 has reduced near surface δ 13 C values by about 0.5‰ from pre-industrial levels (Sonnerup et al., 1999). and high latitude sectors, respectively. This global new production estimate is somewhat higher than the 4.6 Gt C yr −1 from Shaffer (1996) but still only about half as large as more recent estimates (cf. Falkowski et al., 2003). The tradeoff in tuning of vertical exchanges discussed above helps explain our relatively low result. New production in the tuned model is strongly constrained by high latitude surface layer PO 4 and ocean interior PO 4 and O 2 data. This leads to the relatively low value of 0.36 for the high latitude new production efficiency factor, L h f , indicative of strong light and/or iron limitation in this region. Global biogenic calcite production in our solution is 0.97 Gt C yr −1 , thereof 0.83 and 0.14 Gt C yr −1 in the low-mid and high latitude sectors, respectively. This global estimate lies well within the range 0.5-1.6 Gt C yr −1 of other such estimates (Berelson et al., 2007). Our model topography and calcite dissolution length . Steady state, pre-industrial, low-mid latitude simulations of (a) distributions over water depth of calcite dry weight fraction exported at the base of the sediment bioturbated layer for three different sediment model versions: 1) with "exact" carbon pore water chemistry but neglect of the organic carbon fraction in sediment bulk and density (blue line), 2) with "approximate" carbon pore water chemistry and neglect of the organic carbon fraction in sediment bulk and density (green line), and 3) with "approximate" carbon pore water chemistry but with consideration of the organic carbon fraction in sediment bulk and density (red line), (b) profiles across the sediment bioturbated layer at 5000 m water depth of carbonate ion for versions 1 and 3 (blue and red line, respectively) and of dis-  (Table 1) imply that 64% of this calcite production dissolves in the water column; the rest falls on the sediment surface. Model calcite production is constrained strongly by ALK data but also by ocean sediment data (Sect. 3.2.1). This has led to a rather high value of 0.36 for r CalC,m , the rain ratio upper limit parameter and a rather strong temperature dependency of the rain ratio, as expressed by the value of 0.18 for µ. Still, the above low-mid latitude results lead to calcite carbon to organic carbon flux ratios of 0.6, 1.1 and 2.0 at 1000, 2000 and 3000 m depths, respectively, in good agreement with ocean sediment trap results reviewed by Berelson et al. (2007).

Ocean isotope results
Model ocean profiles of δ 18 O w and δ 13 C are shown in Fig. 5, together with mean, data-based sector profiles. Atmospheric processes coupled to evaporation/precipitation force δ 18 O w . Therefore, the distribution of this tracer (Fig. 5a, solid lines) mirrors that of salinity, as can be seen by comparing model and the data based profiles (Figs. 3b and 5a). Our treatment of δ 18 O w yields correct δ 18 O w :S ratios (see also Olsen et al., 2005)  Model δ 13 C results capture much of the vertical structure of ocean observations but exhibit an offset of about 1‰ toward higher values (Fig. 5b). Large model offsets toward higher δ 13 C have also been found in earlier work (Maier-Reimer, 1993). The reason for this offset is not clear. Much of the deep water formation and deep recirculation near Antarctica may occur under the sea ice cover, limiting air-sea gas exchange there and thereby reducing deep ocean δ 13 C. The seasonal cycle in high latitude surface temperatures and sea ice coverage would have a similar effect. Indeed, we were able to achieve a good fit to ocean δ 13 C data (for the fixed atmospheric value of −6.4‰) by sufficiently reducing the gas transfer velocity in our high latitude box. However, we found that such a reduction seriously degraded Geosci. Model Dev., 1, 17-51, 2008 www.geosci-model-dev.net/1/17/2008/ our model fit to ocean 14 C data. Excessive air-sea exchange fractionation at low temperatures is another possible explanation. To illustrate this, we altered the formulation of the Zhang et al. (1995) temperature-dependent, fractionation factor 13 α HCO 3 , such as to yield the same value at 25 • C as the original 13 α HCO 3 but with the weaker, temperature dependency slope of 13 α CO 3 . The results with the altered 13 α HCO 3 show much better agreement with data, including higher surface layer values at low-mid latitudes than at high latitudes (Fig. 5b). We note that the Zhang et al. (1995) results are not based on any measurements at temperatures below 5 • C but we have no other reason to doubt these results. A simulation with a 50% increase in fractionation during new production also yields a considerably better model agreement with ocean mean δ 13 C. However, our results from such a simulation overestimate observed surface layer values and are not consistent with observed δ 13 C in ocean particulate organic matter (Hofmann et al., 2000).
3.2 Pre industrial, steady state solution for the ocean sediment and the lithosphere

Sediment inventories
Values for the oxic remineralization rate scale, λ 0 ox , the anoxic-oxic remineralization rate ratio, β 0 , and the organic rain dependence power, γ , have been chosen to yield model results that satisfy two conditions, given rain rates from the ocean model and the prescribed rain of non-calcite minerals. The conditions are 1) the ocean mean burial fraction for organic matter falling on the sediment surface should be about 0.1 and, 2) the organic matter burial at depths of 1000 m or less should be a fraction of 0.8-0.9 of total ocean organic matter burial (Hedges and Keil, 1995). From our tuning of these parameters (Table 1), the pre-industrial, steady state values for these two fractions are 0.093 and 0.897, respectively. Burial fractions for the low-mid and high latitude sectors are 0.090 and 0.094 and the total organic carbon burial rate, B OrgC,PI , is 0.073 Gt C yr −1 . Model global inventories of erodible and bioturbated layer organic carbon are 130 and 92 Gt C. The best fit value for λ 0 ox agrees with sediment observations for moderate organic carbon rain rates (Emerson, 1985). The best fit values for β 0 and γ lead to much reduced remineralization rates in the sediment under anoxic conditions, as indicated by data (Archer et al., 2002). As an illustration of model sensitivity to anoxic remineralization rate, ocean mean burial fractions are 0.316 and 0.004 when this rate is set to 0 and to the oxic rate.
The value for the calcite dissolution rate constant, k Cal , in Table 1 was chosen to approximate a global inventory of erodible calcite in ocean sediments of about 1600 Gt C (Archer, 1996a) and sublysocline transition layer thicknesses around 1500-2000 m, given biogenic rain rates from the ocean model and the prescribed rain of non-calcite minerals on the ocean surface. The best fit value for k Cal is in the range for which Archer et al. (1998) found good agreement among model results based on linear and non-linear kinetics for calcite dissolution. With this value, model global inventories of erodible and bioturbated layer calcite are 1603 and 1010 Gt C. The model mean calcite dry weight fraction (dwf) is 0.360, close to a data-based estimate of 0.34 (Archer 1996a). The calcite burial rate for the pre-industrial, steady state solution, B Cal,PI , is 0.20 Gt C yr −1 of which 0.13 Gt C yr −1 takes place at water depths greater than 1000 m. Results above give an overall calcite-C to organic carbon-C burial ratio of less than 3 while the corresponding overall sediment inventory ratio is greater than 10. This contrast is explained by the fact that most of the sediment organic carbon is found at shallow depths where the sedimentation velocity is much greater than deeper down.
To test the influence of porosity formulations on our results, we recalculated steady state, sediment calcite carbon and organic carbon inventories using the Zeebe and Zachos (2007) function for the limiting porosity at depth, φ min , together with the Archer (1996a) function for porosity change across the bioturbated sediment layer (see Sect. A1 of Appendix A). This is the same approach as was used by Ridgwell (2007). With this formulation, model global inventories of erodible and bioturbated layer organic carbon are 159 and 80 Gt C and model global inventories of erodible and bioturbated layer calcite are 1659 and 754 Gt C. These results are similar to those reported above from our standard porosity formulation but show slightly lower bioturbated layer inventories. The lower inventory for bioturbated layer calcite agrees somewhat better with a data-based estimate of about 800 Gt C (Archer, 1996a). We also recalculated these inventories for assumed constant sediment porosity and obtained quite different results with considerably lower calcite inventories for example. These findings underline the importance of using an appropriate depth-and composition-dependent porosity formulation in an ocean sediment module.

Bioturbated sediment layer distributions
We compare in Fig. 6 results for three different sediment model versions: 1) With "exact" carbon pore water chemistry (Eqs. A18, A20 and A21 in Appendix A) and neglect of the organic carbon fraction in sediment bulk and density, 2) As version 1 but with "approximate" carbon pore water chemistry (Eq. A23 in Appendix A) and 3) As version 2 but accounting for the organic carbon fraction in sediment bulk and density. Figure 6b shows BL profiles of CO 2− 3 and dissolved CO 2 at 5000 m water depth for version 1 and Fig. 6c shows corresponding HCO − 3 and DIC profiles. The "approximate" CO 2− 3 profile for version 3 only differs slightly from the version 1 profile (Fig. 6b), implying that calcite dissolution rates are very similar in both versions but are slightly enhanced in version 1. This explains the slightly broader and deeper, low-mid latitude sublysocline transition layer in version 1 as compared to those of versions 2 and 3 (Fig. 6a).
www.geosci-model-dev.net/1/17/2008/ Geosci. Model Dev., 1, 17-51, 2008 Lower calcite dry weight fractions in the low-mid latitude, upper ocean sediments for version 3 compared to versions 1 and 2 are explained by the "placetaking" of organic carbon in the version 3 solution (Fig. 6a). These results show that model calcite inventories are more sensitive to the neglect of organic carbon "placetaking" than to the use of "approximate" carbon pore water chemistry, as is also reflected by the low-mid latitude calcite inventories for versions 1, 2 and 3 of 1390, 1397 and 1335 Gt C, respectively. Based on these results and on computation times needed for each of the versions, we choose intermediate complexity, version 3 as our standard sediment model below. Note that this version was also used in the above sediment model calibration. Figure 7b and c shows standard case model distributions with water depth of calcite and organic carbon dwf raining onto and sedimented down out of the BL. Relatively low calcite dwf and the high sedimentation rates (Fig. 7d) at the shallowest water depths reflect our prescribed high rain rate of non-calcite minerals and the relatively high organic matter rain rates at such depths (Eq. A29 in Appendix A). These rains combine to flush the BL relatively rapidly at the shallowest water depths, favoring relatively high organic carbon dwf and organic carbon burial in the model there, as in the real ocean (Emerson, 1985). Relatively high calcite dwf, rapidly decreasing organic carbon dwf and moderate sedimentation rates at intermediate depths above the model CSD's reflect much lower non-calcite mineral and organic matter rain rates, rather constant calcite rain rates, and slower BL flushing, allowing more complete organic carbon remineralization. Rapidly-decreasing calcite dwf, very low organic carbon dwf and low sedimentation rates below the model CSD's reflect calcite dissolution combined with constant non-calcite rain rates and still lower organic carbon rain rates. For the CSD as an upper boundary and the depth where calcite dwf equals 0.1 (the calcite compensation depth) as the lower boundary, sublysocline transition layer thicknesses are 1696 and 1182 m for the low-mid and high latitude zones, respectively. These results agree with the compilations of Archer (1996a), including the sharper transition layer at high latitudes.
Standard case model profiles over the BL for different water depths in the low-mid latitude zone are shown in Fig. 8. The minimum of sediment porosity at intermediate water depths follows from the broad maximum of calcite dwf there and Eq. (A1) of Appendix A (Figs. 8a and 7b). The profiles of calcite undersaturation show that calcite dissolution occurs in the upper few centimeters of the BL which is assumed to be well-mixed in calcite by bioturbation (Fig. 8b). The rapid decrease of organic carbon dwf near the top of the BL for shallow and intermediate water depths reflects rapid, oxic remineralization as compared to slower, anoxic remineralization deeper down in the BL (Fig. 8c). In contrast, increased vertical structure of the organic carbon dwf at deeper water depths can be traced to porosity increases at these depths (Fig. 8a) and the upper boundary condition for the organic carbon solution in the BL (Eq. A7 in Appendix A). Note, however, that organic carbon dwf at the bottom of the BL decreases monotonically with water depth, consistent with the results in Fig. 7c. The upper oxic part of the BL in the lowmid latitude zone is thinnest (4-5 mm thick) at water depths of 500-1000 m where the water column oxygen minimum is found (Figs. 8d and 4b). Below a water depth of about 3000 m, the entire BL is oxygenated. In the high latitude zone where model organic carbon rain rates are larger by about a factor of 5, the oxic part of the BL is even thinner at shallow water depths and anoxic remineralization is found down to depths of 4200 m (not shown). Table 2 lists model weathering rates calculated from overall 12 C and 13 C balances for an assumed pre-industrial steady state (see Sect. 2.7). The lithosphere outgassing estimate is about 1.6% of the present day carbon source from fossil fuel burning and is consistent with other, data-based, lithosphere outgassing estimates (Mörner and Etiope, 2002). From the factor multiplying Vol in Eq. (10) and with calculated values for δ 13 C Cal,PI and δ 13 C OrgC,PI of 1.15 and −23.17‰, respectively, we find that 25.3% and 74.7% of the lithosphere outgassing derives from old organic carbon and carbonate, respectively. The model estimate for total river inflow of inorganic carbon, R C,PI (=2(W Carb,PI +W Sil,PI )=2B Cal,PI ), is 0.40 Gt C yr −1 , in agreement with recent, data-based estimates (Lerman et al., 2007).

Lithosphere results
We find that the pre-industrial carbon sink due to silicate weathering is about 1.3 times as large as the carbon sink due to organic carbon burial and about 2.5 times as large as the carbon source associated with weathering of rocks containing old organic carbon. Pre-industrial, ocean outgassing of CO 2 , R C,PI −B Cal,PI −B OrgC,PI , is 0.13 Gt C yr −1 for the results above. In this steady state solution, ocean outgassing is balanced by net uptake of atmospheric CO 2 , 2W Sil,PI +W Cal,PI −Vol PI −W OrgC,PI .

Greenhouse gas evolutions and warming
To test model performance on decade to century time scales, we made a simulation from 1765 to 2000 AD. Forcing and simulation results as well as comparisons with observations are shown in Fig. 9. Initial conditions were taken from the pre-industrial steady state solution described above. The simulation was forced by prescribed anthropogenic emissions of CO 2 , CH 4 , and N 2 O (Fig. 9b) and by prescribed radiative forcing changes from aerosols, volcanos, variations in solar radiation and the rest of the greenhouse gases for this period (green line in Fig. 9a). For simplicity and due to lack of data, anthropogenic emissions of N 2 O are taken to be proportional to those of CH 4 (in units mol s −1 ) with a proportionality constant of 0.007, chosen for a good model fit to observations (Fig. 9f).
Model evolutions of pCO 2 , pCH 4 , and pN 2 O from 1765 to 2000 agree well with observations. From 1880 to 1950, pCO 2 is slightly underpredicted and during the last decade of the simulation, pCO 2 and, in particular, pCH 4 are slightly overpredicted (Fig. 9e, f). Figure 9a shows the predicted radiative forcing anomalies from the sum of these three greenhouse gases as well as the total radiative forcing change. Modelled atmospheric warming shows polar amplification whereby the high latitude temperature increase is about 50% greater than the low-mid latitude increase (Fig. 9c). In the model this is due largely to poleward retreat of snow and ice lines and the associated ice albedo feedback (Fig. 9d). Although model polar amplification leads to a weaker meridional atmospheric temperature gradient, and thereby a weaker sensible heat transport, the latent heat transport increases slightly  , 2000) and the total radiative forcing changes as the sum of these two (blue), (b) anthropogenic CO 2 emissions from fossil fuel burning (blue; Marland et al., 2007) and from land use change (red) as well as anthropogenic CH 4 emissions (green), (c) simulated low-mid latitude, high latitude and global mean atmospheric temperature change (red, blue and black, respectively), observed global mean atmospheric temperature change (green; Jones et al., 2006) and simulated global mean ocean temperature change (dashed blue), (d) simulated ice line and snow line changes (blue and green, respectively), (e) simulated and observed atmospheric pCO 2 (black and green, respectively) and (f) simulated and observed atmospheric pCH 4 (blue and green, respectively) and simulated and observed atmospheric pN 2 O (red and green, respectively). Observed atmospheric pCO 2 and pCH 4 before 1850 are from Etheridge et al. (1998a) and Etheridge et al. (1998b), respectively, and after 1850 from Hansen and Sato (2007). Observed atmospheric pN 2 O before 1979 are from Meure et al. (2006) and after 1979from IPCC (2007. Land use change CO 2 emissions for 1860 to 1950 are from Houghton (2002). For 1980, Houghton (2002 values were multiplied by a factor of 0.73 to give approximate mean values of 1.4 and 1.6 Gt C yr −1 for the 1980s and 1990s, respectively (IPCC, 2007). Values for 1765 to 1860 were calculated by a linear interpolation from 0 to the 1860 value. Values for 1950-1980 were obtained by multiplying Houghton (2002) values by a factor varying linearly from 1 to 0.73 over this period. Anthropogenic CH 4 emissions from 1860 to 1994 are from Stern and Kaufmann (1998) minus a constant value of 0.04 Gt(CH 4 ) yr −1 , the anthropogenic, pre-industrial value estimate included in the land biosphere emissions (see Sect. 2.2). Values for 1765 to 1860 were calculated by a linear interpolation from 0 to the "corrected" 1860 value. Values for 1995 to 2000 were assigned the "corrected" 1994 value.
due to the greater moisture carrying capacity of warmer air. Likewise, the modeled atmospheric water vapor transport from the low-mid to the high latitude sector increases, from 0.355 Sv in 1765 to 0.371 Sv in 2000.
Modelled mean atmospheric temperature change for the period agrees well with observations and the mean atmospheric warming from 1765 to 2000 is 1.015 • C (Fig. 9c). However, a warming-cooling cycle from about 1920 to 1950 Geosci. Model Dev., 1, 17-51, 2008 www.geosci-model-dev.net/1/17/2008/ is not captured. This cycle may be associated with changes in the Atlantic thermohaline circulation (Zhang and Delworth, 2005). Our ocean circulation and mixing are held fixed to pre-industrial calibrations but even coupled climate models with a dynamic ocean have not simulated the 1920-1950 warming-cooling cycle well with natural and anthropogenic radiative forcing (IPCC, 2007). General model agreement with observed greenhouse gas evolution and global warming speaks well for the model design and calibration, including our choice of a central, 3 • C climate sensitivity. Mean ocean temperature increases much slower than atmospheric temperature due to the slow ocean exchange and large ocean heat capacity (Fig. 9c). Our calculated mean temperature increase for 0-3000 m depths during 1955-1998 of 0.062 • C is more than 50% greater than the Levitus et al. (2005) estimate of 0.037 • C. This may reflect, in part, model overestimate of ocean heat transport for our model calibration (see above) and, in part, a data-based underestimate of ocean heating due to systematic observation errors in some instruments (Gouretski and Koltermann, 2007). After correction for such errors, these authors find an ocean heat content increase (0-3000 m) between 1957-1966 and 1987-1996 that corresponds to a temperature increase of 0.033±0.020 • C. Our result for this period, 0.040 • C, is well within this latest estimate. Figure 10a and c shows anthropogenic CO 2 uptake rates and carbon inventory changes, respectively, for the atmosphere, the ocean and the land biosphere in our 1765 to 2000 simulation. By year 2000, model uptake rates have increased to 4.01, 2.15 and 2.55 Gt C yr −1 for the atmosphere, ocean and land biosphere, respectively. Together with our observationbased estimate for 2000 of emissions from land use change, 1.52 Gt C yr −1 , the net land uptake rate is 1.03 Gt C yr −1 . All these uptake rates agree very well with the latest IPCC consensus (Table 7.1 in IPCC, 2007). Total model carbon inventory increase from 1765 to 2000 is 464.5 Gt C, the sum of increases of 205.1, 118.9 and 140.5 Gt C for the atmosphere, ocean and land biosphere, respectively. In response to the overall 1 • C warming over this 235 year period, the model atmospheric CO 2 sink from weathering increases by 0.018 Gt C yr −1 while the model net oceanic carbon source from river inflow minus burial increases by 0.013 Gt C yr −1 .

CO 2 uptake rates and atmosphere tracer evolutions
Net primary production on land, NPP, increases from from 60 to 71.5 Gt C yr −1 in our 1765 to 2000 simulation (Fig. 10b). This increase is due to the CO 2 fertilization effect (see Sect. 2.6). New production in the ocean, NP, increases from 5.37 to 5.77 Gt C yr −1 over the simulation (Fig. 10b). This increase is due almost entirely to a high latitude sector increase associated with sea ice retreat and more open water available for planktonic production. This physicalbiogeochemical interaction leads to a pCO 2 drawdown in the high latitude surface layer via a net downward transport of in-organic carbon, acting as negative feedback on atmospheric pCO 2 and global warming. However, this feedback is rather weak: when compared to a model simulation for the ice edge free to respond but high latitude new production held constant, the model simulation with this feedback only led to a 0.09 Gt C yr −1 increase in ocean CO 2 uptake and a 1.3 µatm decrease in atmospheric pCO 2 by year 2000.
Model evolutions of atmospheric δ 13 C and 14 C in our simulation agree quite well with corresponding data (Fig. 10d and e). The data-model agreement for atmospheric δ 13 C is even better in recent decades than shown in Fig. 10d: The atmospheric δ 13 C data in the figure is from the high latitude Southern Hemisphere and is less negative in recent decades than the global average, due to Northern Hemisphere fossil fuel sources. For example, mean values for year 2000 are −7.97‰ at Cape Grim (41 • S, 145 • E) and −8.05‰ for an average from four Southern Hemisphere and four Northern Hemisphere stations (Allison et al., 2003). The model value for that year is −8.07‰. The atmospheric 14 C simulation captures well the observed general decrease in 14 C until about 1950, forced by the burning of fossil fuel devoid of 14 C (Suess effect). This result and the excellent modeldata agreement for atmospheric δ 13 C provide more support for model calibrations of air-sea exchange and ocean circulation and mixing. On the other hand, the model assumes constant ocean circulation and mixing and the constant atmospheric 14 C production of the pre-industrial, steady state solution and, therefore, lacks the means to explain the significant, decade scale variability in the atmospheric 14 C observations. Figure 10f shows the simulated atmospheric pO 2 evolution from 1765 to 2000 driven by a large sink from fossil fuel burning and smaller sources from O 2 ocean outgassing and net O 2 production associated with net land biotic carbon sinks (green curve in Fig. 10a minus the red curve in Fig. 9b). In the model simulation, pO 2 decreased by 26.4 µatm from 1993 to 2000 in good agreement with a decrease for that period of 25.5 µatm, as calculated from average observed O 2 /N 2 ratios and estimated N 2 outgassing rates reported by Manning and Keeling (2006). Figure 11 shows the vertical distributions of modeled changes between years 1765 and 2000 of some ocean properties. Although atmospheric temperature warms most at high latitudes, ocean surface layer warming is greater at low-mid latitudes since sea ice cover shields much of the high latitude surface layer from heating and since intense vertical mixing there resists surface layer change (Fig. 11a). However, the intense mixing also heats deeper layers faster at high latitudes. The increased atmospheric water vapor transport forces a saltier and a fresher surface layer in the low-mid and high latitude sectors, respectively (Fig. 11b). The DIC increase in Fig. 11c shows that the high latitude ocean also takes up anthropogenic CO 2 faster. By year 2000, 35.2% (a) simulated CO 2 uptakes by the atmosphere, the land biosphere and the ocean (red, green and blue lines, respectively), (b) simulated net production on land, NPP, and new production in the ocean, NP (green and blue lines, respectively), (c) simulated cumulative CO 2 inventory changes in the atmosphere, on land and in the ocean (red, green and blue lines, respectively), (d) simulated and observed changes in atmospheric δ 13 C (black and green lines, respectively), (e) simulated and observed changes in atmospheric 14 C until 1954 (black and green lines, respectively) and (f) simulated pO 2 changes (black line). Fossil fuel δ 13 C values from 1860 to 1992 for the simulation were taken from Andres et al. (1996). Values before 1860 were set to the 1860 value and values after 1992 were set to the 1992 value. Atmospheric δ 13 C observations were taken from Francey et al. (1999) and atmospheric 14 C observations were taken from Stuiver and Quay (1981). The fossil fuel -O 2 :C mole ratio was taken to be 1.391 in the pO 2 simulation (Keeling et al., 1998). of total ocean DIC increase is found in this sector with only 13.4% of the model ocean volume. This CO 2 uptake forces decreased CO 2− 3 concentrations and pH via ocean carbonate chemistry ( Fig. 11c and f). The surface layer pH decrease over the period is about 0.1 in both sectors, in agreement with other studies (Caldeira and Wickett, 2003) and model CSD's shoal by 43 m and 294 m in the low-mid and high latitude sectors, respectively. However, the shoaling is too little and the period too short to lead to any significant extra sediment calcite dissolution by year 2000. Dissolved oxygen decreases over the period due to less solubility for warmer tempera-tures (Fig. 9d). The vertical structure of this decrease in the high latitude sector also reflects increased new production there and associated increased remineralization. The δ 13 C decrease over the period is forced mainly by exchange with the atmosphere and is greatest in the low-mid latitude surface layer (−1.13‰ vs. −1.68‰ in the atmosphere; Fig. 11c).

Changes in the ocean interior
Geosci. Model Dev., 1, 17-51, 2008 www.geosci-model-dev We present here the results from several long, forced simulations designed to illustrate the workings and behavior of the ocean sediment and lithosphere modules and the interaction with other modules of the DCESS model. All these simulations start from our pre-industrial, steady state solution. In a first group of simulations over 100 000 years, we forced the model by injecting 5000 Gt C CO 2 into the atmosphere over a ∼5000 year time scale at the beginning of the simulation with a maximum CO 2 input rate of 1.11 Gt C yr −1 at simulation year 4000 (grey line in Fig. 12a). Injections of similar sizes and timescales may have occurred in connection with past warming events on Earth (Pagani et al., 2006). In the final simulation over 1.5 million years, we forced the model by doubling the pre-industrial, model lithosphere outgassing (Table 2) to 0.247 Gt C yr −1 at the start of the simulation. For simplicity and relative ease in interpretation, the non-calcite mineral input to the ocean remains unchanged and the rain ratio remains only a function of temperature for all these simulations.
Results from the 100 000 year simulations are shown in Fig. 12. For our standard case simulation, atmospheric pCO 2 rises to 718 µatm at year 7700, drops relatively rapidly over the next 30 000 years or so and decreases more slowly to a pCO 2 of 356 µatm at the end (Fig. 12a). The difference (d) Fig. 12. Model results for simulations over 100 000 years in response to an idealized, "slow" 5000 Gt C injection of CO 2 to the atmosphere. (a) pCO 2 from the standard case simulation (solid line), from a simulation with constant, preanthropogenic weathering (dashed line), from a simulation with constant, preanthropogenic weathering and constant, pre-anthropogenic land biosphere size (dashed-dotted line), and from a "closed" simulation with no interaction with the ocean sediment nor the lithosphere (dotted line). Also shown is the CO 2 source (grey solid line), (b) Standard case simulation, net CO 2− 3 -C fluxes to the atmosphere-ocean system from changes in the bioturbated sediment layer (BL) calcite inventory (solid line) and from the difference between weathering (carbonate plus silicate) and calcite burial down out of the BL (dashed line) as well as net CO 2 -C fluxes to the atmosphereocean system from changes in BL organic carbon inventory (dotted line close to zero line) and from the difference between lithosphere outgassing, weathering of old organic carbon, silicate weathering and organic carbon burial down out of the BL (dashed-dotted line) (c) Standard case simulation changes in carbon inventories for the atmosphere (solid line), the ocean (dashed line), the bioturbated sediment layer (dashed-dotted line) and the land biosphere (dotted line). (d) Standard case simulation, CO 2− 3 -C fluxes from weathering (carbonate plus silicate; upper solid line) and from calcite burial (upper dashed line) as well as CO 2 -C fluxes from lithosphere outgassing and weathering (old organic carbon minus silicate; lower solid line) and as organic carbon burial (lower dashed line). Both burials in (d) have been multiplied by −1 to facilitate comparison. of this value from the pre-industrial 278 µatm represents an "airbourne" fraction of 0.032 of the total CO 2 injection. In this case with climate-dependent silicate weathering, the airbourne fraction will fall to zero (pCO 2 =278 µatm) over an order of magnitude longer time scale (see below). For a simulation with weathering at constant, pre-anthropogenic levels, atmospheric pCO 2 rises to 779 µatm at model year 8450, drops less rapidly over the next 50 000 years or so and then slowly approaches a near steady state by the end with a pCO 2 of 457 µatm and an airbourne fraction of 0.075. Without weathering feedbacks, there is slower neutralization of CO 2 invading the ocean but the system reaches a final (warmer) steady state much sooner. At 10 000 years after the forcing maximum, airbourne fractions for the simulations with and without climate-dependent weathering were 0.13 and 0.18, respectively. This is in the lower end of the range of 0.10-0.30 found after 10 000 years in other climate -carbon cycle models with similar forcing (Archer and Brovkin, 2008).
For a simulation with constant weathering and constant, pre-industrial land biosphere, pCO 2 rises to 911 µatm at model year 8200, drops somewhat more rapidly over the next 50 000 years or so and then also approaches a (slightly warmer) steady state with a slightly higher pCO 2 and airbourne fraction. The difference between the latter two simulations is due to CO 2 taken up and released by the land biosphere. Finally for a simulation with internal cycling only, i.e. complete remineralization/dissolution of all biogenic particles falling on the ocean floor and lithosphere outgassing and weathering set to zero, atmospheric pCO 2 rises to a constant level of 1006 µatm a few thousand years after the CO 2 injection. This leads to a total airbourne carbon fraction of 0.303 (including the small contribution from the simulated atmospheric pCH 4 increase from 0.72 to 1.54 µatm). Corresponding fractions for the ocean and land biosphere are 0.514 and 0.183, respectively.
In the standard case simulation, the increase of CO 2− 3 -C flux to the atmosphere-ocean system from sediment calcite dissolution peaks at model year 5200, about 1200 years after maximum CO 2 injection (Fig. 12b). This dissolution is driven by a decreasing carbonate ion concentrations at middepths due to carbonate ion depletion via reaction with the extra CO 2 invading the ocean (Fig. 13a). This process drives shoaling of the CSD to a minimum of 617 m by model year 7200. The associated increase in CO 2− 3 -C flux to the ocean from enhanced sediment calcite dissolution is first matched and then exceeded in strength and especially in duration by the increase in CO 2− 3 -C flux from the excess of carbonate and silicate weathering over calcite burial (Fig. 12b). Note that this flux increase is driven as much or more by a calcite burial decrease than by a weathering increase (Fig. 12d). For the case above with constant weathering, the calcite burial decrease accounts for all the increase of this flux in the "terrestrial neutralization" process Ridgwell and Hargreaves, 2007). Model calcite burial decrease is accompanied by decreasing sedimentation velocities (Fig. 13c) and decreasing sediment calcite content (Fig. 13b) at mid-depths, both driven by enhanced sediment calcite dissolution. Thus, the impact of this dissolution on neutralization of CO 2 invading the ocean, and thereby on the drawdown of atmospheric pCO 2 , is much greater than would be gauged solely from the decrease in sedimentary calcite inventory.
The greatest decrease of calcite burial at model year 5800 is directly associated with a sedimentation velocity minimum at this time (Figs. 12d and 13c). The sedimentation velocity is actually directed upward for a period of about 3000 years over a depth range exceeding 1000 m. During this event, sediment is being "mined" from below the bioturbated layer. However, the mining rate is so slow in this case that less than a total of 1 Gt C was mined. In response to the CO 2 injection, the carbon inventory in the BL decreases from 1100 Gt C to a minimum of 470 Gt C at model year 12 250 in the standard case simulation (Fig. 12c). This total decrease is composed of a decrease in BL calcite inventory of 645 Gt C and a simultaneous increase of organic carbon inventory in the BL of 15 Gt C. Subsequently, the BL carbon inventory builds up again over the next 40-50 kyr to exceed its original storage as the CSD falls somewhat below its original depth (Fig. 13b).
There is a decrease in net CO 2 -C flux to the atmosphereocean system in the standard case simulation in response to Low-mid latitude model results as functions of water depth for the standard case simulation over 100 000 years in response to the 5000 Gt C injection of CO 2 to the atmosphere. (a) the deviation of carbonate ion from its saturation value (CO 2− 3 -CO 2− 3 (sat)), in mol m −3 , (b) calcite dry weight fraction ((CaCO 3 ) dwf ) in the bioturbated sediment layer and (c) sedimentation velocity, w s , at the base of the bioturbated sediment layer, in cm kyr −1 . the CO 2 injection (Fig. 12b). This is explained by warming and the dominance of silicate weathering over the weathering of old organic carbon in our calibration (Table 2), paired with an increase in organic carbon burial (Fig. 12d). The initial stage of this burial increase stems from increased high latitude new production as sea ice recedes poleward in response to warming, as discussed above. A subsequent, weaker burial increase stems from a switch toward more anoxic remineralization in the sediment in response to this production increase and to decreasing O 2 levels in the ocean, as driven by this production increase and by decreasing O 2 solubility in the warming ocean surface layers. . Standard case model simulation over 1.5 million years in response to "instantaneous" doubling of lithosphere CO 2 outgassing. (a) low-mid latitude, high latitude and global mean atmospheric temperatures (solid red, blue and black lines), (b) atmospheric pCO 2 , (c) phosphorous source from weathering, phosphorus sink from organic matter burial, and the sum of these two (solid blue, green and black lines respectively), (d) carbon sources from lithosphere outgassing and "old" organic carbon weathering (solid violet and blue lines, respectively), carbon sinks from silicate weathering and organic matter burial (solid red and green lines, respectively) and the sum of these sinks and sources (solid black line), (e) the fraction of ocean new production that is buried in the ocean sediment, in terms of carbon and phosphorus (solid blue and red lines, respectively) and (f) atmospheric pO 2 .
The neutralization sink of CO 2 from the CO 2− 3 -C flux increase exceeds the direct sink from the decrease in CO 2 -C flux by almost an order of magnitude directly following the CO 2 injection and subsequently decreases more rapidly than the direct sink (Fig. 12b). The ∼15 kyr, e-folding time scale of this decrease is set mainly by the time needed to replenish the sediment calcite inventory. By model year 37 200, external carbon inputs and outputs balance and the total extra carbon inventory of the combined atmosphere-land biosphere-ocean-ocean sediment system reaches its maximum of 6991 Gt C of which 93.9% resides in the ocean (Fig. 12c). During the last part of the simulation, carbonate burial exceeds the sum of carbonate and silicate weathering and continued drawdown of atmospheric pCO 2 is due solely to enhanced silicate weathering and organic carbon burial (Fig. 12d). Figure 14 shows results from the 1.5 million year simulation. Since this simulation was forced by a doubling of lithosphere CO 2 outgassing, a model steady state can only be reached when the net CO 2 sink from the sum of silicate weathering, weathering of old organic carbon and organic carbon burial also doubles. We adopted a simple Q 10 -based, climate dependency with Q 10 =2 for both types of weathering. If such a dependency also held for organic carbon burial, the final steady state, global mean temperature would be 25 • C, exactly 10 • C greater than our pre-industrial value of 15 • C. On the other hand if organic carbon burial remained constant over this simulation, a global mean temperature of 32.8 • C (and a pCO 2 well over 10 000 µatm) would be required for steady state, given the results in Table 2 and above for the pre-industrial, steady state calibration.

A 1.5 million year simulation
By the end of the simulation, a new steady state has nearly been reached with global mean temperature of 24.5 • C and an atmospheric pCO 2 of 2636 µatm ( Fig. 14a and b). Organic carbon burial does increase by almost a factor of two over the simulation but the size and structure of this increase are not simple functions of the global warming but rather reflect changes in biogeochemical cycling in the model ocean and ocean sediment that accompany this warming (Fig. 14c, d and e). There is an increase in organic carbon and phosphorus burial during the first ∼30 kyr of the simulation due to increased high latitude new production as sea ice recedes poleward and disappears in response to initial warming. As global temperatures continue to increase, weathering input of phosphorus exceeds phosphorus loss through burial with a maximum of this imbalance centered at about model year 100 000. This leads to increased ocean phosphate inventories, new production and burial. There are further burial increases over the rest of the simulation due to the switch discussed above toward more anoxic remineralization in the sediment, above all in response to decreasing O 2 levels in the ocean, driven by higher new production and by decreasing O 2 solubility for warmer conditions. This O 2 feedback drives an increase in model carbon burial fraction from 1.36% to 2.01% over the simulation (Fig. 14e) and thereby explains most of the simulated organic carbon burial increase. This feedback and associated burial increase also act to inhibit anoxia in the mid-depth ocean, by limiting oxygen consumption there and by limiting ocean phosphate inventories and thereby ocean new production.
The forcing and model response over this simulation have a small but well-defined effect on atmospheric O 2 concentrations (Fig. 14f). Initally, the oxidation of the reduced part (25.3%) of the increased lithosphere outgassing forces decreasing atmospheric pO 2 until about model year 370 000. Subsequently, pO 2 rises as the carbon burial Geosci. Model Dev., 1, 2008 www.geosci-model-dev.net/1/17/2008/ fraction increases, leaving behind an increasing surplus of O 2 produced in new production over that consumed in remineralization and reaction with reduced gases in the BL. This modeled increase in atmospheric pO 2 would continue beyond 1.5 million years until the pO 2 is large enough to lead to high enough O 2 concentrations in the ocean to reduce the carbon burial fraction enough to restore O 2 balance. However, the carbon and phosphorus cycling would also be affected leading to a further climate drift. A proper treatment of the coupled carbon, nutrient, oxygen and climate system over such long time scales is beyond the scope of the present model and would require, for example, a treatment of sulfur cycling (Berner, 2006).

Discussion and conclusions
We have put considerable emphasis on and much effort into calibrating the DCESS model to pre-industrial conditions by fitting to available data. We believe that simple models should be calibrated to observations to the greatest extent possible and to the results of complex models to the least extent necessary. After all, complex models may do rather poorly in tests against observations. For example, a number of ocean carbon cycle models fail to simulate well observed ocean 14 C distributions and, therefore, present day ocean circulation and mixing (Matsumoto et al., 2004). Our calibrations to complex model results have been limited to the choice of a CO 2 fertilization factor for the land biosphere and to the dependence of methane consumption on methane concentration in the atmosphere. Neither of these important items can be constrained sufficiently well by available observations. For simplicity, the ocean module of the DCESS model has no dynamics and ocean circulation and mixing are prescribed at well-calibrated, pre-industrial calibration levels. This approach serves to keep down the degrees of freedom in our low-order model but at the expense of being able to deal with climate change associated with simulated changes in ocean circulation. Most coupled climate models exhibit weakened meridional overturning circulations for global warming (Gregory et al., 2005). On the other hand, reconstructions of past conditions show weaker overturning for cooler conditions and data from high-resolution climate archives like ice cores indicate that there is considerably more climate variability due to ocean circulation changes when climate is colder than during warmer interglacial periods (Toggweiler and Russell, 2008;Grootes and Stuiver, 1997). To deal well with colder climates, a low-order Earth System Model should include (at least) a continental ice sheet module, an improved description of the land biosphere model and a more sophisticated sea ice formulation.
We put particular effort into the development and calibration of a new, semi-analytical ocean sediment model, described in detail in Appendix A. Our goal was a model suffi-ciently simple for rapid simulations while retaining sufficient complexity and flexibility to deal with organic and inorganic sediment fractions at all ocean depths. As opposed to most coupled model work up to now (Heinze et al., 1999;Ridgwell, 2007), we consider oxic and suboxic remineralization of organic matter, the dependence of bioturbation rate on organic carbon fluxes to the sediment and dissolved oxygen concentrations and the dependence of sediment remineralization rates on bioturbation rates.
We have also included the dependency of porosity on sediment composition. In times of changing sediment composition, this effect feeds back upon sediment composition and sediment burial rate via effective pore water diffusion and via sedimentation velocities at the base of the bioturbated sediment layer. However, we have not modeled the effect of organic carbon remineralization on calcite dissolution in the sediment. A proper treatment of this would require explicit consideration of chemical species involved in anoxic remineralization that is beyond the scope of the present work (but see Ridgwell, 2007). However, as in many other coupled models, dissolution above the calcite saturation depth of biogenic calcite particles sinking out of the surface layer is included as a simple function of depth. This may represent dissolution in the sediment as well as in the water column. Our ocean sediment module would also be well suited for use in more complex Earth System Models.
In conclusion, we developed, calibrated and tested against data a new, low-order Earth System Model designed to be comprehensive, flexible and fast. This DCESS model should serve as a useful tool for studies of past, present and future global change on time scales of years to millions of years, in particular for climates as warm as or warmer than the present one.

Appendix A
The ocean sediment module

A1 General features
The sediment module addresses calcium carbonate dissolution and (oxic and anoxic) organic matter remineralization in 0.1 m thick bioturbated layers (BL) on sediment segments for each of the 2×55 ocean model layers. The sediment is composed of calcite, CaCO 3 , non-calcite mineral, NCM, and reactive organic matter fractions. To a good approximation, CaCO 3 and NCM fractions are taken to be well mixed in the BL by bioturbation, D b , but the reactive organic matter fraction varies over the BL due to relatively rapid remineralization. For each sediment segment, NCM fluxes are prescribed based on data and calculated particulate organic matter ( the sediment are then removed from/returned to these ocean layers. For each sediment segment, solutions are found for steady state profiles of reactive organic carbon (OrgC), porewater O 2 and pore-water CO 2− 3 . From the latter, dissolution of CaCO 3 in the BL is calculated and used in steady state and time dependent calculations of calcite dry weight fraction, (CaCO 3 ) dwf .
Sediment porosity, φ, defined as the ratio of pore volume to total volume in the sediment, is a key property in the sediment module. Empirical data show φ to be a function of calcite dry weight fraction, as represented here by (Archer, 1996a), where the vertical coordinate ζ is taken positive downward from the sediment surface (ζ =0), φ min =1−(0.483+0.45 (CaCO 3 ) dwf )/2.5 and α=0.25 (CaCO 3 ) dwf +3 1− (CaCO 3 ) dwf in centimeters.
Another recent estimate for the limiting porosity at depth in the deep sea, φ min , is φ min =(φ min,n + (CaCO 3 ) dwf )/(1+ (CaCO 3 ) dwf ) where =(φ min,c −φ min,n )/(1−φ min,c ) and φ min,n and φ min,c , the limiting porosities for pure non-calcite and pure calcite sediment, are taken to be 0.88 and 0.61 (Zeebe and Zachos, 2007). Here we use the Archer (1996a) formulation for φ min as standard. Another key sediment property here is the sediment formation factor, F s (ζ ), needed to calculate bulk sediment diffusion coefficients of pore water solutes. These coefficients are reduced from free water molecular diffusion values by complex sediment structure (tortuosity).

A2 Organic carbon
The general, steady state equation governing bioturbated layer profiles of solid reactive organic carbon (in moles cm −3 of solid sediment) is where λ ox and λ anox are the oxic and anoxic remineralization rates for the parts of the BL above and below the depth ζ o where O 2 goes to zero. To proceed in the solution of Eq. (A2), one could substitute for φ(ζ ) using Eq. (A1) and seek analytical solutions to the resulting nonhomogenious, second order differential equation for the regimes above and below ζ o . However we choose to take another approach and to divide the BL into k sublayers with assumed constant φ for each sublayer i. From Eq. (A2) this leads to the simpler governing equation for each solution layer j (note that, in general, there will be one more solution layer than there are sublayers since there will be an oxic and an anoxic solution above and below ζ o in the sublayer where ζ o is found). For each BL sublayer, mean φ and F s are calculated by taking averages over sublayer thicknesses, i.e.
ζ =ζ i φ(ζ ), (φ(ζ )) −3 dζ . We found that sufficient resolution of φ and F s could be achieved for all (CaCO 3 ) dwf with seven sublayers bounded at ζ =0, 0.2, 0.5, 1, 1.8, 3.2, 6 and 10 cm. Corresponding φ i and F s,i are plotted in Fig. A1. The general solutions of Eq. (A3) are and (OrgC) anox,j (ζ ≥ ζ 0 )=A 2,j exp(s 3 ζ ) + B 2,j exp(s 4 ζ ) (A5) where s 1 , s 2 = 0.5 w s /D b ± (w s /D b ) 2 + 4λ ox /D b where F OrgC is the rate of organic carbon rain (in mol cm −2 s −1 ) at the sediment surface from the ocean module. The boundary condition at the bottom of the BL reduces to a vanishing OrgC gradient there, due to the assumption of no bioturbation below ζ b . All this leads to 16 equations (two of them non-linear) in the 17 unknowns A 1,j , B 1,j, A 2,j, B 2,j and ζ o for our seven sublayer case and the complete solution must await simultaneous solutions of the coupled organic carbon-dissolved oxygen problem given below. For a completely oxygenated BL, corresponding to weak organic carbon fluxes to the sediment, Eq. (A4) is the general solution and the complete organic carbon solution may be obtained at once by solving for A 1,j and B 1,j (now j =1, . . . 7) from the above condition at ζ =0, matching conditions at the BL sublayer boundaries, and from d(OrgC) ox,7 /dζ =0 at ζ =ζ b . We can calculate profiles of sediment density, ρ s , and organic matter dry weight fraction, (Org) dwf , from the organic carbon solution with the assumption that the densities of the calcite and non-calcite mineral fractions, ρ min are the same (2.7 g cm −3 ) such that where M C is the molecular weight of carbon, WR Org is the total weight/carbon weight ratio in organic matter (2.7) and ρ Org is the organic matter density (1.1 g cm −3 ).

A3 Remineralization and dissolved oxygen
The total organic carbon remineralization over the BL, RM OrgC (in moles cm −2 s −1 ), is where depending on where ζ o is found. RM OrgC also gives the DIC flux from organic carbon remineralization that is fed into the appropriate ocean module layer. The organic carbon burial flux, BF OrgC , is equal to (1−φ k )w s (OrgC) at ζ =ζ b and, for these steady state solutions, is also equal to F OrgC −RM OrgC . Finally, the organic carbon burial fraction is BF OrgC /F OrgC . For simplicity we also apply the above model with the same sediment remineralization rates to organic phosphorus raining on the sediment surface to obtain RM OrgP and BF OrgP . This leads to identical sediment remineralization and burial fractions for organic P as for organic C (recall however that, in the model, C:P ratios in POM rain to the sediment surface are different from C:P ratios in new production due different water column remineralization scales for "nutrient" and "carbon" fractions of POM). The sink of ALK related to organic matter remineralization in the sediment is r AlkP RM OrgP and this sink is subtracted from the appropriate ocean module layer. Also, for simplicity at this stage in model development, 13,14 C contents in DIC leaving the sediment at specific locations and times are coupled directly to 13,14 C contents of organic carbon and calcite particles raining on the sediment surface at those locations and times. In future work we will solve explicitly for sediment distributions of these isotopes (as well as for pelagic and benthic 18 O). This will allow explicit calculation of 13,14 C contents in DIC leaving the sediment.
For the general case of λ ox =λ anox and ζ o <ζ b , the above solutions for sediment remineralization and burial depend upon ζ o from the solution for pore water O 2 (in moles m −3 of pore water). The general, steady state governing equation for pore water O 2 is where D O 2 is the temperature-dependent, free solution, molecular diffusion coefficient of O 2 (in cm 2 s −1 ), the factor (φF s (ζ )) −1 describes the attenuation reduction of this diffusion coefficient as discussed above, r OC,S is the sediment mole ratio of oxygen consumed per carbon remineralized, taken to be 1.4 for consistency with ocean module choices. For our BL sublayer solution approach, Eq. (A12) reduces to and the general solution to Eq. (A13) is with s 1 , s 2 from Eq. (A6) and A 1 and B 1 from the (OrgC) ox solution above. Specific solutions are obtained for all ζ ≤ζ o by applying the boundary/ matching conditions. These are O 2 =O 2,ocean , the O 2 concentration from the appropriate ocean module layer, at ζ =0, O 2 =0 at ζ =ζ o and continuous concentrations and fluxes at BL sublayer boundaries. For our approach, the flux conditions are A final boundary condition is based on the following: During anoxic remineralization below ζ o reduced species are produced (like H 2 S). Some of these species will precipitate with available metals (like Fe) at a rate limited by the rain rate of such metals to the sediment. However, in general in the steady state, by far most of the reduced species produced diffuse upward to be oxidized by O 2 near ζ =ζ o . Thus, to a good approximation, this total extra oxygen demand can be equated with r OC,S times the total anoxic remineralization of organic carbon in the BL below ζ 0 and can be taken to be a line sink at ζ =ζ o . In the model, this extra oxygen demand is supplied by downward diffusion. Thus, at ζ =ζ o , For a partially-anoxic BL, this leads to n algebraic equations (four of them non-linear) in n unknowns A 1,j , B 1,j , A 2,j , B 2,j , C j , D j and ζ o . Solutions for these equations, and thus for the coupled OrgC-O 2 problem, are obtained by iterating toward the solution value of ζ o (we use the Matlab function fzero for this). For a completely oxygenated BL, the O 2 problem for our seven layer case reduces to 14 linear algebraic equations in C 1,2,..,7 , and D 1,2,..,7 , given A 1 and B 1 from the OrgC solution for an oxygenated BL. These equations derive from the sediment surface boundary condition above, 12 concentration and flux matching conditions at the sublayer boundaries and a no flux condition at ζ =ζ b . In either case, total O 2 consumption in the BL is r OC,S RM OrgC , equal to the O 2 flux to the BL (φ 1 D O 2 (F s,1 ) −1 d(O 2 ) 1 /dζ at ζ =0). This is then the O 2 flux subtracted from the appropriate ocean module layer. Organic carbon rain to the sediment surface provides sustenance for the benthic fauna. Thus, bioturbation rates due to the actions of this fauna should depend upon organic carbon rain rates but should also be attenuated as dissolved oxygen concentrations become very low. Our parameterization for these relationships is based on the bioturbation estimates and approach of Archer et al. (2002) but with D b constant over the BL for simplicity: where O 2,ocean is the ocean O 2 concentration at the sediment surface and the bioturbation rate scale, D 0 b , the organic carbon rain rate scale, F 0 OrgC , and O 2,low are taken to be 1×10 −8 cm 2 s −1 , 1×10 −12 mol cm −2 s −1 , and 20×10 −3 mol m −3 , respectively (Archer et al., 2002). Furthermore, we assume that oxygen remineralization rates in the BL scale as bioturbation rates (and thereby as organic carbon rain rates; Archer et al., 2002), such that λ ox =λ 0 ox D b /D 0 b , where λ 0 ox is an oxic remineralization rate scale.
The anoxic remineralization rate will depend upon the specific remineralization reactions involved. For example, denitrification will occur below the oxic layer in the BL but above the layer where sulfate reduction occurs at a slower rate than denitrification. Furthermore, more organic rain would be associated with a more anoxic BL and a shift toward sulfate reduction. Therefore we assume here that λ anox =βλ ox whereby β is taken to decrease for increasing organic carbon rain rate such that β=β o (F OrgC /F 0 OrgC ) γ . In Sect. 3.2, λ 0 ox , β o and γ are constrained by organic carbon burial observations.

A4 Calcite
The dissolution rate of calcite in the BL depends on the carbonate ion concentration in the pore water and the calcite concentration in the solid phase. In turn, the carbonate ion distribution in the BL is related to the distributions of the other inorganic carbon species there via carbonate chemistry. Steady state equations describing this system (Boudreau, 1987;Archer, 1996b) can be reduced for each of our BL sublayers to three coupled equations for CO 3 , HCO 3 and CO 2 (for convenience we drop charges here). The first of these equations may be obtained from subtracting the resulting DIC equation from the resulting ALK equation, under neglect of a small borate contribution (Archer, 1996b), to yield, where D CO 2 and D CO 3 are temperature-dependent, molecular diffusion coefficients for CO 2 and CO 3 , respectively, (CO 3 ) sat is CO 3 saturation with respect to calcite, calculated as in Sect. 2.4, for ambient T , S and pressure from appropriate ocean module layers, and where k Cal is a calcite dissolution rate constant, and M Cal is the molecular weight of calcite (100 g mole −1 ). Note that for simplicity and to a good approximation, the mean BL sediment density, ρ sm , is used in Eq. (A19). The term in Eq. (A18) multiplied by Cal,i describes the effect of calcite dissolution with linear dissolution kinetics and is nonzero only when CO 3 <(CO 3 ) sat . Archer et al. (1998) found that, for a proper value for k Cal , calcite dissolution can be equally well described with simpler linear kinetics as with often-used, non-linear kinetics.
A second BL sublayer equation may be obtained by subtracting the resulting ALK equation from twice the resulting DIC equation (again neglecting the borate correction) to yield, d 2 (HCO 3 ) i /dζ 2 −2(D CO 2 /D HCO 3 )d 2 (CO 2 ) i /dζ 2 =0 (A20) where D HCO 3 is the temperature-dependent, free solution, molecular diffusion coefficient for HCO 3 . The third BL sublayer equation results from carbonate chemistry yielding where K 1 and K 2 are the first and second apparent dissociation constants for carbonic acid as functions of ambient T , S and pressure from appropriate ocean module layers.
To be complete, the coupled system described by Eqs. (A18), (A20) and (A21) should also include the effects of oxic and anoxic respiration. Oxic respiration adds CO 2 to this system, reducing CO 3 concentrations and tending to enhance calcite dissolution. The effect of anoxic respiration is more subtle (Boudreau, 1987;Archer, 1996b) but adds mostly HCO 3 to the system, increasing CO 3 concentrations and tending to reduce calcite dissolution. A correct treatment of this complex system would require detailed treatment of species involved in anoxic respiration and is beyond the scope of the present work. Thus we have chosen to neglect these respiration effects in the present treatment. But to some extent these effects are included in the water column dissolution of calcite as parameterized in Sect. 2.4.
The CO 3 concentration in the model BL equals the adjacent ocean CO 3 concentration, and there is no calcite dissolution in the model, if CO 3,ocean ≥(CO 3 ) sat (as calculated from ambient DIC, ALK, T , S and pressure). If CO 3,ocean <(CO 3 ) sat , we seek solution of Eqs. (A18), (A20) and (A21) as follows: First, guided by the mathematical nature of the problem, we assume solutions for BL sublayer CO 2 of the form (CO 2 ) i =A 3,i exp(s 5,i ζ ) + B 3,i exp(−s 5,i ζ ) Second, we obtain general solutions for CO 3,i and HCO 3,i from Eqs. (A18) and (A20). Third, we obtain specific solutions for all three carbon species by applying boundary and matching conditions: concentrations at ζ =0 from the ocean module, vanishing gradients at ζ =ζ b and matching concentrations and fluxes at sublayer boundaries. Fourth, we use these specific solutions to check if the solutions satisfy Eq. (A21) everywhere in the BL. Fifth, we repeat steps 3 and 4 for different choices of s 5,i until the solutions satisfy Eq. (A21) well in the BL, giving the steady state, pore water solutions for a specified (CaCO 3 ) dwf .
A particularly simple solution for Eq. (A18) is obtained if the second term in that equation can be neglected relative to terms one and three, yielding This is equivalent to neglecting detailed carbonate chemistry within the BL. The general solution to Eq. (A23) is (CO 3 ) i =(CO 3 ) sat + A 4,i exp(s 6,i ζ ) + B 4,i exp(−s 6,i ζ )(A24) where s 6,i =( Cal,i ) 0.5 . Specific solutions are then found by applying the corresponding boundary/matching conditions to those listed above. In Sect. 3.2.2, we compare the results of this simplified solution to those from the complete carbonate chemistry solutions above. For a calcite steady-state in the BL, the calcite flux to the sediment surface, F Cal , must be balanced by the sum of calcite dissolution within the BL, DIS Cal , and the calcite flux down (or up) through ζ =ζ b , BF Cal , (all in mol m −2 s −1 ) such that From the above we have where where again ρ sm has been used for simplicity and to a good approximation. DIC and ALK fluxes from calcite dissolution that are fed into the appropriate ocean module layer are given by DIS Cal and 2DIS Cal , respectively. Furthermore, BF Cal = (1 − φ k )w s ρ sm (CaCO 3 ) dwf /M Cal (A28) The steady state sedimentation velocity at ζ =ζ b follows from overall mass balance, where F NCM is the flux of non-calcite mineral to the sediment surface (g cm −2 s −1 ). This w s is then used in the organic carbon problem above and in the steady state calcite problem below. The sources of F NCM are atmospheric dust input, river input of terrigenous material and non-calcite, biogenic minerals, in particular opal produced mainly by diatoms. To capture a decrease in F NCM from the coast toward open ocean background values, as would be expected from the sum of the above sources, we take where NCF is the open ocean, non-calcite flux, CAF is the amplification factor at the coast (i.e. at z=0) and λ slope is an e-folding, water depth scale. In our simple model, λ slope "stands in" for a distance from the coast, given typical continental slope topography. We take λ slope to be 300 m such that at 2000 m depth at the outer edge of the slope, the "near shore" component of F NCM has been reduced to only about 1 per mil of its value at the coast. The value for NCF is taken to be 0.95×10 −7 g m −2 s −1 (0.3 g cm −2 kyr −1 ), estimated as a mean, open ocean value from Fig. 3a in Archer (1996b). The value for CAF is taken to be 20 to yield realistic, model "shelf" sedimentation rates of about 20 cm kyr −1 (see also Sect. 3.2.1).

A5 Coupled solutions
Steady state solutions for sediment calcite concentrations can now be sought as follows: First, initial guesses are made for (CaCO 3 ) dwf , ρ sm and w s and the porosity profile is calculated from Eq. (A1) using (CaCO 3 ) dwf . Second, based on the above formulations and flux and concentration forcing from the ocean module, organic carbon remineralization and calcite dissolution are calculated. Third, based on these results, ρ sm , w s and, subsequently, calcite burial are calculated. Fourth, the overall calcite balance in the BL (Eq. A25) is now checked and the system is iterated with new choices of (CaCO 3 ) dwf until this balance is fulfilled, yielding the steady state value of (CaCO 3 ) dwf and, subsequently, solutions for all model sediment components. Time-dependent solutions, forced by time-varying fluxes and concentrations from the ocean module, are obtained from time-dependent, mass, calcite and organic carbon inventory balances, where the overbars indicate means taken over the BL. Note that for simplicity and to a good approximation, the (OrgC) dwf in Eq. (A33) and below is used in a BL-mean sense.
From the above and mass balance (Eq. A31), timedependent ρ sm and w s are ρ sm = (OrgC) dwf W Org ρ Org +(1−(OrgC) dwf W Org )ρ min (A34) and w s ={(F Cal −DIS Cal )M Cal +F NCM +(F OrgC −RM) M C W Org +ρ sm z b dφ/dt}/{ρ sm (1−φ 7 )} Note the "extra" term in Eq. (A35) involving the time rate of change of porosity, a term that can be important during periods of rapidly changing (CaCO 3 ) dwf (cf. Eq. A1). A corresponding term for the time rate of change of ρ sm has been neglected since non-carbonate and carbonate mineral fractions are assumed to have the same density and since these fractions are typically much larger that the organic matter fraction.
From calcite and organic carbon balances (Eqs. A32 and A33), we have where ρ sm and w s are taken from Eqs. (A34) and (A35). We also solved a simpler version of the calcite part of our sediment model that considers only the calcite and noncalcite mineral sediment fractions, thereby neglecting terms in Eqs. (A29), (A31), (A34) and (A35) associated with the relatively-small, organic sediment fraction. In Sect. 3.2.2, we compare results from this simpler version with the more complete version above. Equation (A35) describes the time evolution of the advection of sediment across the bottom boundary of the BL. This advection will of course vary for changing fluxes to the sediment surface but can also vary significantly for changing dissolution rates of calcite (and to a lesser extent for changing remineralization rates of organic carbon). For large enough dissolution, w s can reverse signs and sediment buried earlier will be reintroduced into the BL. However, large dissolution will also lead to large reductions in calcite dry weight fraction and, thereby, lead to significant porosity increases (cf. Eq. A1). Thus total sediment solid phase inventory in the BL will decrease as the calcite inventory decreases, reducing the need to reintroduce sediment from below to satisfy mass balance. When such chemical erosion occurs, Eqs. (A36) and (A37) are modified such that ρ sm , (CaCO 3 ) dwf , (OrgC) dwf and φ 7 in the second term on the right hands side of these equations are assigned values for the sediment being reintroduced from below the BL, values taken from the appropriate previous solutions.
The time-dependent problem is solved by stepping Eqs. (A36) and (A37) forward in time from given initial conditions (for example from steady state solutions), for specified (time-dependent) ocean forcing and boundary conditions and with the use of Eqs. (A34) and (A35). At each time step, steady state solutions for organic carbon/dissolved oxygen and pore water carbonate are calculated, based on the ocean forcing and boundary conditions, calcite fraction (and associated porosity), organic carbon fraction (and associated sediment density) and sedimentation velocity. Calcite dissolution and organic carbon remineralization that result are used Geosci. Model Dev., 1, 17-51, 2008 www.geosci-model-dev.net/1/17/2008/ in calculating changes of calcite and organic carbon fractions and thereafter in updating these fractions, sediment density and sedimentation velocity.