LOSCAR : Long-term Ocean-atmosphere-Sediment CArbon cycle Reservoir Model

Introduction Conclusions References


Introduction
Various carbon cycle models that are computationally inexpensive have been developed in the past, in particular box models of the ocean's carbon cycle (e.g.Sarmiento and Toggweiler, 1984;Walker and Kasting, 1992;Toggweiler, 1999;Munhoven and Francois, 1996;K öhler et al., 2005).However, only a few studies have considered sediments (e.g.Sundquist, 1986;Sigman et al., 1998;Ridgwell, 2001)  long-term carbon cycle fluxes such as carbonate and silicate rock weathering (e.g.Shaffer et al., 2008).The LOSCAR model (Long-term Ocean-atmosphere-Sediment CArbon cycle Reservoir model) closes this gap.LOSCAR is primarily designed to efficiently compute the partitioning of carbon between ocean, atmosphere, and sediments on time scales ranging from centuries to millions of years.LOSCAR includes various biogeochemical tracers such as total dissolved inorganic carbon (TCO 2 ), total alkalinity (TA), phosphate (PO 4 ), oxygen (O 2 ), stable carbon isotopes (δ 13 C), and %CaCO 3 in sediments.Based on the predicted tracer distributions, different variables are computed including atmospheric CO 2 , ocean pH, calcite and aragonite saturation state, calcite compensation depth (CCD) and more.LOSCAR also allows for changes in the major ion composition of seawater, including the seawater Mg/Ca ratio, which is critical for paleo-applications.The major ion seawater composition affects thermodynamic quantities such as equilibrium constants and solubility products, which in turn affect the predicted ocean carbonate chemistry and atmospheric CO 2 .
We have previously published several applications of LOSCAR dealing, for instance, with future projections of ocean chemistry and weathering, pCO 2 sensitivity to carbon cycle perturbations throughout the Cenozoic, and carbon/calcium cycling during the Paleocene-Eocene Thermal Maximum (PETM) (Zeebe et al., 2008(Zeebe et al., , 2009;;Uchikawa and Zeebe, 2008;Stuecker and Zeebe, 2010;Uchikawa and Zeebe, 2010;Komar and Zeebe, 2011;Zeebe and Ridgwell, 2011).The subject of the present contribution is the detailed description of the model including numerical architecture, processes and parameterizations, tuning, and examples of input and output.It may appear that publishing model applications before a detailed model description is putting the cart before the horse.One of the reasons for this is that the journals interested in publishing the model applications have little or no interest in publishing a detailed model description.
Journals that provide a forum for technical model descriptions are rare, and so the recent appearance of Geoscientific Model Development has encouraged me to provide a coherent model description of LOSCAR that will hopefully be useful for the readership of the journal, as well as the users of the model.On the other hand, publishing a Figures few model applications before the detailed model description also has an advantage.LOSCAR, for example, has been extensively tested by now and several bugs and scientific/numerical issues have already been fixed.LOSCAR's main components include ocean, atmosphere, and marine sediments.The model architecture, main components, model variables, and process parameterizations will be described in the following.Finally, two input/output examples will be presented, one dealing with anthropogenic fossil fuel emissions, the other with carbon release during the PETM.

Architecture
The basic architecture of the model is fairly simple.For all model variables y i , i.e. all tracers in all compartments (atmosphere, ocean boxes, and sediment boxes), a system of coupled, first-order ordinary differential equations is solved: where t is time, NEQ is the total number of equations, i = 1,2,...,NEQ, and F s are known functions.Note that for most applications, the derivatives (right-hand side of Eq. 1) do not explicitly depend on the independent variable t.For given start (initial) conditions y 0 at t start , the equations are then integrated forward in time over the interval from t start to t final .Standard numerical procedures (solvers) for this sort of problem are available.One thing to keep in mind is that the equations solved in LOSCAR are typically stiff and involve different time scales, which requires a solver for stiff problems with adaptive stepsize control.The solver implemented in the C version of LOSCAR is a fourth-order Rosenbrock method with automatic stepsize adjustment (Press et al., 1992).
Once the initial conditions y 0 and derivatives F 's have been supplied, the solution of the problem is usually straightforward.However, setting up y 0 and F requires some Figures

Back Close
Full work.In the following, the individual model components will be described and expressions will be given for individual F 's that enter Eq. ( 1).The current setup includes two different model versions: a "modern" version and a Paleocene/Eocene version ("P/E"version for short).

Geometry
The global ocean is geometrically divided in LOSCAR into separate ocean basins representing Atlantic, Indian, and Pacific Ocean (plus Tethys in the P/E-version).In turn, each ocean basin is subdivided into surface, intermediate, and deep ocean (Fig. 1).
In addition, the model includes a generic high-latitude box (H-box), representing cold surface waters without reference to a specific location (cf.Walker and Kasting, 1992;Toggweiler, 1999).As a result, the total number of ocean boxes is NB = 10 in the modern version and NB = 13 in the P/E-version.Box areas and volumes are given in Table 1.The modern ocean geometry in LOSCAR is not unlike the one used by Walker and Kasting (1992).However, Walker and Kasting (1992) combined the warm surface and thermocline waters each into a single reservoir for a total of 6 boxes to represent the global modern ocean.
The modern and Paleocene/Eocene ocean bathymetry in LOSCAR is based on Menard and Smith (1966) and Bice and Marotzke (2002), respectively.The bathymetry determines the surface area and volume of ocean boxes (Table 1) and the surface area-depth level relationship of the sediment boxes (Sect.6).

Ocean tracer equations
Let y k be a subset of y (Eq.1), representing ocean tracer variables including TCO 2 , TA, PO 4 , etc. (in this particular order).2,...,2NB for TA, k = 2NB + 1,2NB + 2,...,3NB for PO 4 , and so on.If the total number of ocean tracers is NOCT, then the total number of equations for all ocean tracers and boxes is NOCT × NB.The differential equation for an ocean tracer y k may be written in the general form: where V k is the volume of box k and F 's are fluxes due to (thermohaline) circulation and mixing, air-sea gas exchange (e.g. in case of TCO 2 ), biological uptake and remineralization, riverine/weathering input, and sediment fluxes.The first three flux terms on the right-hand side of Eq. ( 2) will be explained in the following subsections; the riverine/weathering and sediment flux terms will be explained in Sects.4 and 6.

Circulation, mixing, and air-sea gas exchange
Given a prescribed ocean circulation-and mixing scheme, F thm is of the form: where T is the volume transport of the conveyor circulation and m l k are mixing coefficients between boxes l and k (Fig. 2, Table 2).The box indices j and l are set by the prescribed circulation/mixing scheme (Fig. 2).The coefficients m l k represent bidirectional mixing, hence m l k = m kl .The air-sea gas exchange term reads: Figures

Back Close
Full where κ as is the air-sea gas exchange coefficient for CO 2 and A k is the area of surface box k; pCO 2 a and P CO 2 k is the atmospheric pCO 2 and the pCO 2 in equilibrium with dissolved CO 2 in surface box k, respectively.The index k runs over all surface boxes and for tracers such as TCO 2 and T 13 CO 2 .

Biological pump
The biological uptake and recycling of tracers is parameterized based on phosphate (PO 4 for short).For instance, net uptake of PO 4 in the low-latitude surface ocean (equivalent to particle export flux from the mixed layer) is calculated as: where the parameter f epl describes the efficiency for PO 4 uptake in the low-latitude surface boxes, m j k × [PO 4 ] j is the flux of PO 4 supplied by upwelling/mixing from the underlying intermediate box j into the surface box k. (Note that in the model, the conveyor transport T does not directly supply nutrients to the warm surface waters; it does so, however, to the cold surface waters, see Fig. 2).If f epl were to approach 1.0 (100 % efficiency), all upwelled PO 4 would be converted to sinking particles and the phosphate concentration of surface box k would be zero.In the model, as well as in reality, f epl is usually less than 1.0 (Table 2).influx entering the H-box, simple Michaelis-Menten kinetics prevent PO 4 from becoming negative.Caution is therefore advised when increasing F pph because the actual high-latitude export flux may be less than the value assigned to F pph .The high-latitude export flux is remineralized in the deep boxes.The fluxes of TCO 2 and TA due to biological uptake and recycling are computed based on PO 4 using a given Redfield-and rain ratio (Table 2).Note that there is a small contribution to alkalinity changes from organic carbon production and respiration as a result of nitrate uptake and release (e.g.Zeebe and Wolf-Gladrow, 2001).The major contribution to alkalinity changes in the model is associated with CaCO 3 fluxes.Of the total CaCO 3 export flux, the larger fraction is destined for accumulation or dissolution in sediments, the latter of which returns total carbon and alkalinity to the ocean (Sect.6).
A smaller fraction of the CaCO 3 export is assumed to dissolve in the water column (Table 2).This assumption yields better agreement with observed TA fields and is consistent with observations and modeling studies indicating substantial water column dissolution above the lysocline (e.g.Archer et al., 1998;Milliman et al., 1999;Feely et al., 2002).

Carbonate and silicate weathering
Weathering of carbonate rocks on the continents takes up atmospheric CO 2 and supplies calcium and bicarbonate ions to the ocean: Hence two moles of carbon and one mole of Ca 2+ enter the ocean for each mole of CaCO 3 weathered, raising ocean TCO 2 and TA by two units each (Fig. 3).If the CaCO 3 riverine/weathering influx is denoted by F cc (in units of mol CaCO 3 y −1 , see where k = 1,...,NOC runs over all low-latitude surface boxes and NOC is the number of corresponding ocean basins.Note that in steady state, subsequent precipitation of CaCO 3 in the ocean (Reaction 6 backwards) releases the same amount of CO 2 back into the atmosphere as was taken up during weathering.In other words, the CO 2 for carbonate weathering essentially originates from the ocean (Fig. 3).As a result, although the addition of Ca 2+ and 2 HCO − 3 increases ocean TCO 2 : TA in a 2:2 ratio, on a net basis CaCO 3 weathering increases ocean TCO 2 : TA in a 1:2 ratio because one mole of CO 2 returns to the atmosphere.If influx equals burial, carbonate weathering thus represents a zero net balance for atmospheric CO 2 .The steady-state balance is restored after a perturbation on a time scale of 5 to 10 ky and is referred to as "calcite compensation" (Broecker and Peng, 1987;Zeebe and Westbroek, 2003).
Weathering of silicate rocks and simultaneous uptake of atmospheric CO 2 may be described by: If the CaSiO 3 riverine/weathering influx is denoted by F si (in units of mol CaSiO 3 y −1 , see Table 3), then: Note that silicate weathering removes 2 moles of CO 2 from the atmosphere for each mole of CaSiO 3 weathered.Subsequent precipitation and burial of CaCO 3 (Reaction 6 backwards) releases one mole of CO 2 back to the atmosphere, the other mole is buried in the form of CaCO 3 in sediments (Fig. 3).In steady state, the balance is closed by long-term CO 2 input to the atmosphere from volcanic degassing.Putting it the other Figures

Back Close
Full way, the CO 2 released by volcanoes is balanced by silicate weathering and subsequent carbonate burial in the ocean (Fig. 3).The net reaction is: The steady-state balance for silicate weathering is restored after a perturbation on a time scale of 10 5 to 10 6 yr.This process also restores the partial pressure of atmospheric CO 2 in order to maintain a mass balance of long-term carbon cycle fluxes (e.g.Berner et al., 1983;Zeebe and Caldeira, 2008).
The restoring time scale for silicate weathering is much longer than for carbonate weathering for two reasons.First, silicate weathering requires whole-ocean TCO 2 to adjust, whereas carbonate weathering only requires the ocean's carbonate ion concentration to adjust.On average, the modern TCO 2 inventory is about 20 times larger than mean-ocean [CO 2− 3 ].Second, carbonate weathering fluxes have been estimated to be about 2.5-times larger than silicate weathering fluxes (Table 3).Combined, this gives a factor of about 50, which, multiplied by the calcite compensation time scale of 10 ky, gives 500 ky, which is about right.

Weathering feedback
The feedback between atmospheric CO 2 and weathering fluxes of carbonates and silicates is parameterized in the model using the following equations (see Walker et al., 1981;Berner et al., 1983;Walker and Kasting, 1992): where the superscript "0" refers to the initial (steady-state) value of the weathering flux and pCO 2 , respectively.The parameters n cc and n si control the strength of the weathering feedback (Table 3).For a detailed discussion of the uncertainties associated Figures

Back Close
Full with the weathering parameterization, see Uchikawa and Zeebe (2008).As mentioned above, in steady state, the silicate weathering flux balances the CO 2 degassing flux from volcanism: Thus, the long-term steady-state pCO 2 of the model is set by picking a value for pCO 0 2 , which drives the system towards equilibrium via the silicate weathering equation (Eq.12).Only when the actual model pCO 2 equals pCO 0 2 , will the fluxes be balanced

Atmosphere
The model variable tracking the inventory of atmospheric carbon dioxide, C atm , is related to the partial pressure of CO 2 in the atmosphere by: where q 0 = (2.2×10 15/12) mol ppmv −1 converts from ppmv to mol.Note that for numerical scaling purposes (see Sect. 7.4), C atm is normalized to order 1 in the program by multiplying by (A oc × 100) −1 (arbitrary factor).The differential equation for C atm may be written in the general form (an analogous equation holds for 13 C atm ): where F 's are fluxes due to air-sea gas exchange, volcanic input and weathering (see Sect. 4), and possible carbon input sources.The air-sea gas exchange term for the atmosphere reads: Introduction

Conclusions References
Tables Figures

Back Close
Full where κ as is the air-sea gas exchange coefficient for CO 2 and A k is the area of surface box k; pCO 2 a and P CO 2 k is the atmospheric pCO 2 and the pCO 2 in equilibrium with dissolved CO 2 in surface box k, respectively.The sum runs over all surface boxes.In case of carbon input to the atmosphere from fossil fuel burning or from other carbon sources, for instance, during the PETM, a term of the form: is added where C in is in units of Pg C y −1 .

Sediments
The sediment model calculates %CaCO 3 (dry weight) in the seafloor-bioturbated (mixed) sediment layer of thickness h s as a function of sediment rain, dissolution, burial, and chemical erosion (for more details see Zeebe and Zachos, 2007).The model is particularly useful for long-term integrations and has been constructed similar to other models of this class (e.g.Keir, 1982;Sundquist, 1986;Sigman et al., 1998).However, the current model also includes variable porosity -a feature critical to simulating strong dissolution events that lead to sediment erosion, such as expected for the future or during the PETM (Zeebe and Zachos, 2007;Zeebe et al., 2008Zeebe et al., , 2009)).

Chemical erosion
When dissolution of CaCO 3 exceeds the rain of CaCO 3 plus refractory material such as clay, the sediment column shrinks and previously deposited, underlying sediment Introduction

Conclusions References
Tables Figures

Back Close
Full is reintroduced into the top layer and exposed to dissolution.This is referred to as chemical erosion.As a result, significantly more CaCO 3 is available for dissolution during erosion than originally contained in the top sediment layer.Once the top layer is entirely filled with clay, the sediment column is "sealed" and dissolution ceases.In order to fill the sediment top layer with clay, the sediment volume that was initially filled with CaCO 3 + pore water must be replaced by clay + pore water.Thus, if the sediment porosity φ is constant, the ratio of total CaCO 3 available during erosion to the mass contained in the original surface layer is given by: (Broecker and Takahashi, 1977) where f ci and (1 − f ci ) are the initial CaCO 3 and clay dry weight fraction of the sediment, respectively.However, if porosity varies with %CaCO 3 (as observations show, see below), the ratio of total dissolved to initial CaCO 3 is given by: where φ 0 and φ 1 are the porosities of a pure clay and calcite layer, respectively.The factor (1−φ 0 )/(1−φ 1 ) is of the order 0.3-0.5 and therefore significant as it reduces the erodible CaCO 3 from below the bioturbated layer by 50-70 % compared to the constant φ estimate (Archer, 1996).

Variable porosity
In many locations, it has been observed that porosity decreases with greater CaCO 3 fraction f c (e.g.Mayer, 1991;Herbert and Mayer, 1991;deMenocal et al., 1993).That is, sediment with high CaCO 3 content has a higher concentration of total solids per Introduction

Conclusions References
Tables Figures

Back Close
Full unit volume than low carbonate sediment.The relationship between φ and f c for a sediment layer composed of CaCO 3 , clay, and pore water is given by: where ).The sediment model uses variable porosity as given by Eq. ( 20) and values for φ 0 and φ 1 as given in Table 3.Note that using the non-linear Eq. ( 20) in the model leads to the correct ratio of initial to erodible CaCO 3 (cf.Eq. 19, which was independently derived based on the geometry of the problem), while a linear relationship, for instance, would not.

Sediment model equations (single sediment box)
At every time step, calcite and clay rain of solid density ρ s is added to the top sediment layer of thickness h s (see Table 3 for values).Dissolution of calcite reduces the calcite content and net accumulation is hence rain minus dissolution.At the bottom of the sediment mixed layer, an amount equal to net accumulation is removed via burial.If dissolution of CaCO 3 exceeds the rain of CaCO 3 plus clay, chemical erosion occurs.
The sediment model thus has to provide equations to calculate rain, dissolution, burial, and erosion.At variable porosity, the top layer can be separated into pure calcite plus pore water at porosity φ 1 (volume = A h 1 ) and pure clay plus pore water at porosity φ 0 (volume = A h 2 ).For variable porosity, the model equations can be conveniently written in terms of dh 1 /dt.Conversion to df c /dt merely requires multiplication by a factor (see below).
In case rain exceeds dissolution, no erosion needs to be considered and we can write for a single sediment box: Introduction

Conclusions References
Tables Figures

Back Close
Full where r cs is the calcite rain rate, r d is the calcite dissolution rate, and w c is the calcite burial rate.All rates refer to volume of calcite plus pore water per unit area and time (unit m y −1 ) at porosity φ 1 .Total rates of calcite + clay + pore water are denoted by r s and w.Burial equals rain minus dissolution, i.e. w = r s − r d , and the condition for no erosion is w > 0. The rain rate of calcite, r cs , depends on the carbon export, the rain ratio, and the fraction of water column dissolution.In the low latitudes, for instance, r cs is given by: where F epl is the low-latitude carbon export (in units of mol m −2 yr −1 ), r rain is the rain ratio (C org : CaCO 3 ), ν wc is the CaCO 3 fraction dissolved in the water column (Table 2), converts from mol CaCO 3 to kg CaCO 3 .The rain rate of refractory material, r rs , is calculated correspondingly based on F rrf (Table 3) and the total rain r s is given by r s = r cs + r rs .
The dissolution rate, r d , is calculated as: where R d is given by the following expression at modern seawater Mg/Ca ratio (Sigman et al., 1998): where that the effective rate parameters K sd and n sd relate bottom water undersaturation to dissolution rate (Keir, 1982;Sundquist, 1986;Sigman et al., 1998;Zeebe and Zachos, 2007, see Table 3 for values).They are not be confused with reaction parameters relating porewater undersaturation to dissolution rate such as the calcite reaction order n (typically n = 4.5).
Finally, an expression is needed for the calcite burial, w c , as a function of total burial w.The thickness of the pure calcite layer within ∆z(= w ∆t) can be expressed as f c ∆z (1 − φ) but also as 1 • ∆h 1 (1 − φ 1 ) (calcite fraction = 1), which gives: or expressed per unit time as a rate: As a result, all rates have now been expressed by model-predicted quantities and thus by inserting Eqs. ( 22), (23), and ( 27) into (21), the change in calcite content per time step can be computed.Because we took care of all individual porosities, the relationship between φ and f c , Eq. ( 20), is obeyed automatically.
In case of erosion (w < 0), it can be shown that: where f ci and φ i is the initial calcite fraction and porosity, respectively, and r rs is the clay rain rate (see above).The total dissolution of pure calcite can be derived as: In other words, all calcite in ∆z and calcite rain is dissolved.In addition, calcite is being replaced by the clay in ∆z and by the clay rain (equivalent calcite is also dissolved).Finally, the sediment model can also be formulated in terms of f c by simply multiplying by a factor: where and with F φ = (φ 1 − φ 0 )/(1 − φ 1 ).

Sediment model equations (all sediment boxes)
Let y n be a subset of y (Eq.1), representing the CaCO 3 dry fraction (f c ) in sediment boxes at different depth levels in the different ocean basins.If the total number of depth levels is NSD and the total number of ocean basins is NOC (Table 1), then the total number of equations for all sediment boxes is NSD × NOC.Based on Eq. ( 30), the differential equation for the CaCO 3 dry fraction in sediment box j is: where each sum runs over all sediment boxes j located within the area and depth range of ocean box k.The surface area of sediment box j is denoted by A sed j .

Ocean carbonate chemistry
Carbonate chemistry parameters for modern seawater composition are calculated based on equilibrium constants on the total pH scale (Lueker et al., 2000;Zeebe and Wolf-Gladrow, 2001).The C program uses a simplified and fast numerical routine to compute CO 2 parameters from TCO 2 and TA (Follows et al., 2006).If applied properly, the method yields accurate results that are essentially identical to those obtained with standard routines (Zeebe and Wolf-Gladrow, 2001).The method was originally devised to compute modern carbonate chemistry parameters in biogeochemical models where conditions change little between consecutive time steps (Follows et al., 2006).This is not necessarily always the case in LOSCAR and can lead to failure in rare cases.
For instance, if the model is initiated with a very high TA/TCO 2 ratio, the calculated H + concentration may become negative.The user is warned in such instances and is advised to change the initial conditions.Again, such cases are probably rare.In fairness, it should also be noted that non-standard chemistry conditions (which can occur in LOSCAR), are beyond the original intend of the method (Follows et al., 2006).Apart from the limitation mentioned above, the method is easy to implement, sufficiently accurate, and computationally efficient.Introduction

Conclusions References
Tables Figures

Back Close
Full

Paleocene/Eocene ocean chemistry
Paleocene/Eocene seawater conditions were different from modern conditions owing to factors such as temperature and major ion composition of seawater, including the seawater Mg/Ca ratio (e.g Tyrrell and Zeebe, 2004).These factors can significantly affect thermodynamic quantities such as equilibrium constants and solubility products, which in turn have a major impact on the predicted ocean carbonate chemistry and atmospheric CO 2 .The chemistry routines implemented in LOSCAR allow for variations in, for instance, temperature, salinity, and the concentrations of Mg 2+ and Ca 2+ in seawater.For example, due to warmer surface and bottom water temperatures in the late Paleocene and Eocene, the calcite saturation concentration at a bottom water temperature of 14-17 • C during the PETM is quite different from the modern at 2 • C (see Fig. 3 of Zeebe and Zachos, 2007).This effect is included in LOSCAR by using temperaturedependent equations for the solubility product of carbonate minerals (Mucci, 1983).Pressure corrections for solubility products and equilibrium constants are based on Millero (1995) (Tyrrell and Zeebe, 2004;Zeebe et al., 2009).The effect of seawater Mg 2+ and Ca 2+ on the first and second dissociation constant of carbonic acid is estimated using sensitivity coefficients (Ben-Yaakov and Goldhaber, 1973): and finally: Sensitivity parameters for the effect of Mg 2+ and Ca 2+ on K * are (Ben-Yaakov and Goldhaber, 1973): With these sensitivity parameters, and the modern and paleo-concentrations of Mg and Ca 2+ (see above), the correction to equilibrium constants (Eq.38) can be applied.
Seawater Mg 2+ and Ca 2+ also affect the calcite solubility product, K * sp , and thus the steady-state deep-sea [CO 2− 3 ].Following Mucci and Morse (1984), the stoichiometric solubility product drops with decreasing seawater Mg/Ca ratio.In other words, Eocene K * sp would have been smaller and, given roughly constant deep-sea saturation state, [CO 2− 3 ] would also have been smaller than modern.The data of Mucci and Morse (1984) may be fitted to an equation of the form: where m = modern, α = 0.0833, and x = Mg/Ca.Using modern and P/E-values for [Mg 2+ ] and [Ca 2+ ] as given above, the stoichiometric solubility product of calcite would have been reduced by about 30 %, compared to modern.Another important consequence of changes in oceanic Ca 2+ , for instance, is its effect on the ocean carbon inventory.The long-term carbon inventory and carbonate chemistry of the ocean-atmosphere system is controlled by atmospheric CO 2 and the Introduction

Conclusions References
Tables Figures

Back Close
Full balance between riverine flux and carbonate burial (Zeebe and Caldeira, 2008) ] was 20 mmol kg −1 .In the model, this leads to a pre-PETM ocean carbon inventory that is similar to the modern value, despite a higher baseline atmospheric CO 2 at the time.

Temperature sensitivity
The initial temperature of each individual ocean box is set at the start of the run.
Throughout the run, temperature can be held constant, be manipulated based on user input, or be computed based on a simple expression for the sensitivity of temperature to changes in atmospheric CO 2 as calculated by the model (cf.Archer, 2005).In order to provide a flexible and numerically stable option, the C version of the program includes temperature as an ocean tracer variable.The temperature of ocean box k (T C,k in • C) is assumed to respond to a change in pCO 2 with a certain time lag and relax towards equilibrium temperature.The equilibrium temperature of box k is given by: where the superscript "0" refers to the initial (steady-state) temperature and pCO 2 , respectively, and s is the prescribed temperature increase per doubling of atmospheric CO 2 .The parameter s as used here is conceptually similar to what is generally referred to as "climate sensitivity".However, the precise meaning of s will have to be defined properly for each specific application in the context of the time scales and feedbacks involved (see Zeebe, 2011).Figures

Back Close
Full The differential equation for the temperature of ocean box k then reads: where τ n is the relaxation time, which can take on three different values depending on whether k refers to a surface, intermediate, or deep box (Table 2).

Numerics
As mentioned above, the equations solved in LOSCAR are typically stiff and require an appropriate solver for the problem.The LOSCAR C-version uses a fourth-order Rosenbrock method with automatic stepsize adjustment (Press et al., 1992).For these kind of solvers, it is critical to scale the variables properly.Thus, variables have been scaled to order 1, if necessary, by multiplying by arbitrary factors before passing to the solver.This includes, for instance, atmospheric carbon and temperature (see Sects. 5 and 7.3).The carbonate dissolution flux, R d is proportional to the square root of the CaCO 3 fraction f c (Eq. 24).It turned out that during strong dissolution, f c occasionally became negative when the CaCO 3 fraction approached zero.This issue has been eliminated (in most cases) by using a linear relationship between f c and R d when f c drops below a certain threshold value f c sml .The threshold value can be changed by the user and should be increased if f c still becomes negative during a run.Another option is to increase the solver accuracy by reducing the value of ε slv (the default value is usually not very accurate).
LOSCAR is quick.Running the fossil fuel scenario over 1250 yr (Fig. 5) using the LOSCAR C code compiled under Linux with gcc 4.4.3,without optimization and default ε slv , takes less than 2 s wall clock time on a current standard desktop machine with Intel Core2 Duo E8500 @3.16GHz (no other CPU-demanding processes running).The computational efficiency makes LOSCAR an ideal tool for multi-parameter variations that require a large number of model runs (e.g.Zeebe et al., 2008Zeebe et al., , 2009))

Tuning
In order for LOSCAR to provide model output that resembles observations, several model parameters require tuning.This includes mixing coefficients, biological export fluxes, remineralization fraction (intermediate vs. deep box), rain ratio, and water column dissolution (see Table 2).The tuning is based on comparison between modelpredicted variables and modern observations.For example, parameters were tuned by requiring the ocean tracer variables TCO 2 and TA in the various model boxes to match GLODAP data, averaged over the area and depth range of the corresponding boxes (Key et al., 2004).Note that TCO 2 data were corrected for anthropogenic carbon by subtracting 45 and 25 µmol kg −1 from the surface and intermediate values, respectively (see below for δ 13 C-corrections).The agreement between model and data is satisfactory (see Fig. 4).As a result, the global preindustrial TCO 2 inventory in LOSCAR is 35 920 Pg C vs. 35 760 Pg C based on GLODAP data (Key et al., 2004).Similarly, model PO 4 and oxygen were compared to data summarized in the World Ocean Atlas (WOA05, 2005).Again, the agreement between model and data is adequate, except perhaps for the oxygen content in intermediate boxes, which appears to be underestimated by the model.This could be improved.However, it would come at the expense of a larger mismatch in the deep boxes.This was avoided because for our LOSCAR applications so far, the properties of the deep boxes were more important than those of the intermediate boxes.
Another variable used for parameter tuning is the stable carbon isotope composition of TCO 2 (δ 13 C TCO 2 ), which was matched to the data of Kroopnick (1985).Note that due to the ocean's uptake of fossil fuel carbon (which is isotopically light, i.e. depleted in 13 C), the ocean's δ 13 C TCO 2 is continuously dropping (so-called Suess effect).
Thus, for preindustrial tuning, the early δ 13 C-data sets are more useful than the most recent ones, which are increasingly contaminated with anthropogenic carbon.Nevertheless, Kroopnick (1985) estimated that surface ocean δ 13 C TCO 2 had already dropped by ∼0.5 ‰ and that the average δ 13 C TCO 2 of the preindustrial surface ocean was about Figures

Back Close
Full 2.5 ‰.This surface value was used for model parameter tuning (Fig. 4).As a result, the preindustrial δ 13 C of atmospheric CO 2 is −6.38 ‰ in LOSCAR vs. −6.30to −6.40 ‰ based on ice core and firn data (e.g.Francey et al., 1999).Adequate model values for the steady-state carbonate ion concentration in the deep boxes are important for both the ocean and the sediment model component.After parameter tuning, the preindustrial deep-sea [CO 2− 3 ] as predicted by LOSCAR and calculated based on GLODAP data (Key et al., 2004) are in good agreement (Fig. 4).The preindustrial inventory of CaCO 3 in the seafloor-bioturbated sediment layer (in units of carbon) is about 800 Pg C, close to the value of more complex models (e.g.Archer et al., 1998).
In summary, after model-data comparison including all variables shown in Fig. 4, the values for the parameters labeled "tuned" in Table 2 were obtained.The preindustrial (steady-state) pCO 2 in the model was set to 280 ppmv by assigning this value to pCO 0 2 , which drives the system towards the desired steady-state pCO 2 via the silicate weathering equation (Eq.12).Regarding the Paleocene/Eocene model setup, several key parameters such as deep-sea [CO 2− 3 ] and the calcite compensation depth (CCD) before and during the PETM have been discussed elsewhere and are not repeated here (see Zeebe et al., 2009, Supplementary Information).The pre-PETM inventory of CaCO 3 in the seafloor-bioturbated sediment layer (in units of carbon) is about 620 Pg C. The initial (steady-state) partial pressure of atmospheric CO 2 was set to 1000 ppmv in our P/E-simulations.Although this value falls within the (large) range of available proxy estimates, it is somewhat arbitrary.The user is welcome to change the initial pCO 2 value in the P/E-setup.Figures

Back Close
Full LOSCAR can read in files that supply a time series of fossil fuel emissions in order to project future changes in atmospheric CO 2 , surface ocean pH, calcite and aragonite saturation, and other variables (cf.Zeebe et al., 2008, Supporting Online Material).For example, Fig. 5 shows results obtained with LOSCAR for a fossil fuel emission scenario with a total carbon release of 1000 Pg C over 500 yr.Note that the results differ slightly from those in Zeebe et al. (2008) because ocean temperature was held constant here for simplicity.The initial conditions from which the scenario was started are the preindustrial steady-state conditions shown in Fig. 4. No changes in the biological pump were assumed (PO 4 is constant).The temperature of the low-and high-latitude box is 20 and 2 • C, respectively (Table 2).This temperature difference is mostly responsible for the difference in carbonate ion concentration ([CO 2− 3 ]) and saturation state (Ω) between low-and high-latitude surface boxes.Note that while TCO 2 in the surface boxes responds immediately to the fossil fuel carbon release, there is a delay in TA, which only starts rising once sediment dissolution commences and the calcite compensation depth (CCD) starts shallowing (cf.Ilyina et al., 2009).

Paleocene-Eocene Thermal Maximum
Using appropriate boundary conditions, LOSCAR can also be used to simulate time intervals or events of the past such as the PETM.During the PETM, a large mass of carbon was released into Earth's surface reservoirs (e.g.Dickens et al., 1995;Zachos et al., 2005;Dickens, 2011), while surface temperatures rose by 5-9 • C within a few thousand years.Figure 6 shows results for a PETM scenario with an initial carbon input of 3000 Pg C over a few thousand years, which yielded close agreement with observations (for more details, see Zeebe et al., 2009).Note that the time interval of the integration now covers 200 ky (t = 0 refers to the P/E boundary), rather than a few Introduction

Conclusions References
Tables Figures

Back Close
Full millennia as in the previous example.Changes in boundary conditions compared to the modern setup include a Paleocene/Eocene bathymetry (Bice and Marotzke, 2002), addition of the Tethys ocean, different seawater chemistry (see Sect. 7.2), and different circulation patterns (see Fig. 2).Furthermore, the PETM simulations use different initial conditions for e.g.temperature, steady-state pCO 0 2 , weathering fluxes etc. (Tables 2     and 3).At steady-state pCO 0 2 of 1000 ppmv, but similar carbonate mineral saturation state as in the modern ocean, the steady-state pH of the Paleocene/Eocene ocean would have been lower than modern (Fig. 6).Because of higher seawater Ca 2+ and the effect of Mg/Ca on the solubility product of calcite, the initial carbonate ion concentration in the P/E-simulations is substantially lower than in the modern ocean (cf.Sect.7.2).As a result, steady-state TCO 2 and TA are similar to modern values despite higher pCO 2 (Fig. 6).Note that the Atlantic CCD shoals dramatically during the event, while there is little response in the Pacific CCD, consistent with observations (Zachos et al., 2005;Zeebe et al., 2009;Leon-Rodriguez and Dickens, 2010).The "overshoot" of the CCD, i.e. the fact that its position is deeper at t > 80 ky than its initial position, is a direct consequence of the weathering feedback (see Sect. 4) and is also in agreement with observations (e.g. Kelly et al., 2005).At t > 80 ky, atmospheric pCO 2 is still elevated over the initial pCO 0 Introduction

Conclusions References
Tables Figures

Back Close
Full LOSCAR is a useful tool to tackle carbon cycle problems on various time scales as demonstrated in earlier applications that dealt with future projections of ocean chemistry and weathering, pCO 2 sensitivity to carbon cycle perturbations throughout the Cenozoic, and carbon/calcium cycling during the PETM (Zeebe et al., 2008(Zeebe et al., , 2009;;Uchikawa and Zeebe, 2008;Stuecker and Zeebe, 2010;Uchikawa and Zeebe, 2010;Komar and Zeebe, 2011;Zeebe and Ridgwell, 2011).The present contribution has provided a coherent description of the LOSCAR model.The description will hopefully be beneficial to the readership of the journal, as well as users of the model.I anticipate that future applications will reveal the full spectrum of problems suitable to be studied with LOSCAR.The LOSCAR source code in C can be obtained from the author by sending a request to loscar.model@gmail.com.Full       (Keir, 1982;Sundquist, 1986;Sigman et al., 1998;Zeebe and Zachos, 2007); n sd is not to be confused with the calcite reaction order n, relating porewater undersaturation to dissolution rate (typically n = 4.5).Introduction

Conclusions References
Tables Figures

Back Close
Full   Walker and Kasting, 1992;Toggweiler, 1999).The P/E steady-state setup is inspired by observations and modeling studies of Paleocene/Eocene circulation patterns with significant deep water formation in the Southern Ocean (e.g.Bice and Marotzke, 2002;Thomas et al., 2003;Lunt et al., 2010).Note that a transient contribution of North Pacific deep water (not shown) was included in our PETM simulations (Zeebe et al., 2009).All ocean boxes (except H-box) in the modern and P/E-setup are coupled to sediment boxes (schematically indicated only in the bottom panel for the Pacific by the gray shaded area).Introduction

Conclusions References
Tables Figures

Back Close
Full Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Then k = 1,2,...,NB for TCO 2 , k = NB + 1,NB + Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | The fraction f rim of the export flux is remineralized in the intermediate box, whereas the fraction (1 − f rim ) is remineralized in the deep box.The high-latitude PO 4 export flux can be set directly by assigning a value to the flux parameter F pph .If the value chosen is too large to be supported by the total PO 4 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | dh 1 dt = r cs − r d − w c (21) sd and n sd are "effective" rate parameters (see below), [CO 2− 3 ] sat and [CO 2− 3 ] is the carbonate ion concentration at calcite saturation and in the bottom water, respectively, and c 0 = 1 mol kg −1 so that R d is in units of mol m −2 yr −1 .It is important to note Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 33)where j = 1,2,...,NSD for the first ocean basin (Atlantic), j = NSD + 1,NSD + 2,...,2NSD for the second ocean basin (Indian), and so on.In case of dissolution, TCO 2 and TA are returned to the ocean, giving rise to the sediment source term in the ocean tracer equation (cf.Eq. 2): Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | ) where ∆K * is the change in K * due to the relative change in concentration, ∆c i /c i , of component i .Using ∆c/c = (c − c m )/c m , where m = modern, it follows: Discussion Paper | Discussion Paper | Discussion Paper | ∆K

Fig. 1 .
Fig. 1.Schematic representation of the LOSCAR model (Paleocene/Eocene configuration).A = Atlantic, I = Indian, P = Pacific, T = Tethys ocean, H = High-latitude surface, L = Lowlatitude surface, M = interMediate, D = Deep box.Weathering fluxes and gas exchange with the atmosphere (Atm) are indicated by "w" and "g", respectively.Steps on the faces of ocean boxes indicate sediments (Sed).

Fig. 4 .
Fig. 4. Computed model tracers and observations used for LOSCAR parameter tuning for modern (preindustrial) configuration, see text for details.

Fig. 5 .
Fig. 5. Example of a fossil fuel emission scenario simulated in LOSCAR: total release of 1000 Pg C over 500 years (see Zeebe et al., 2008).Results shown slightly differ from those in Zeebe et al. (2008) because ocean temperature was held constant here for simplicity.L = Low-latitude, M = interMediate, D = Deep, H = High-latitude.A = Atlantic, I = Indian, P = Pacific.Note that the step in the Atlantic calcite compensation depth (CCD, panel h) is due to the spacing of sediment-box depth levels in the model (adding more sediment boxes would make the curve smoother).

Fig. 6 .
Fig. 6.Example of a PETM carbon release scenario simulated in LOSCAR: initial release of 3000 Pg C over a few thousand years (see Zeebe et al., 2009).L = Low-latitude, M = interMediate, D = Deep, H = High-latitude.A = Atlantic, I = Indian, P = Pacific, T = Tethys.See text for details.

Table 1 .
Model-architecture and ocean geometry parameters.

Table 3 .
Weathering and sediment model parameters.