ocean SOEST

FRED T. MACKENZIE: Research Models: MAGic
[To view the example paper abstract, click here.]

Model Methodology:

In this model, we have broadly partitioned material between the atmosphere, seawater, oceanic crust, and upper mantle. The continents (cratons) are subdivided into surface (including continental shelf and shallow subsurface) and deeper burial reservoirs that cycle at different rates. The following reactions govern transfers and transformations of material within each reservoir (sediments, dissolved components, gases, and crystalline rock): (1) weathering, (2) precipitation, diagenetic, and reduction-oxidation reactions, and (3) hydrothermal-magmatic reactions. Weathering reactions (table 1) consume primary crystalline Na-, K-, Ca-, and Mg-silicates and Ca- and Mg-carbonates through attack by both carbonic and sulfuric acids and produce clay minerals and dissolved cations, bicarbonate, sulfate, and dissolved silica. Simple dissolution of evaporites is also included in this group. Precipitation reactions (table 2) include formation of primary and secondary carbonate (calcite and dolomite) and evaporite (anhydrite, epsomite, halite and sylvite) phases. Low temperature diagenetic and redox reactions include the dissolution of these carbonates, the formation of aluminosilicate phases through so-called "reverse weathering" reactions, formation of organic matter and its subsequent oxidation by oxygen and sulfate, and the oxidation of sulfide (represented as pyrite) by oxygen. Hydrothermal-metamorphic-magmatic reactions (table 3) include reconstitution of primary silicates by reaction of carbonate minerals with clays, decarbonation of calcium and magnesium carbonates, dolomite formation by reaction of calcium carbonate with Mg-aluminosilicates (the ultimate source of Mg for dolomitization), seawater-basalt high and low temperature alteration reactions, and thermogenic sulfate reduction to form pyrite, which ultimately returns the oxygen originally consumed in pyrite oxidation. The latter process is problematic and discussed in more detail below. These reactions are written using normative mineral phases or components to represent overall reaction stoichiometry. These simplifications ignore many intermediate weathering and alteration products, but preserve an overall mass balance during model operation.

Table 2: Low Temperature Precipitation, Diagenetic and Redox Reactionstable2

Table 3: Hydrothermal-Magmatic Reactionstable3


Steady State Fluxes

Because the number of coupled reservoirs and components is large, the model can be best understood in terms of principal subcycles that are themselves linked through transfers of carbon dioxide and oxygen: (1) Calcium-Magnesium-Silicate-Carbonate-CO2, (2) Sodium-Potassium-Silicate-Carbonate-Chloride-CO2, (3) Organic Carbon, (4) Sulfur, and (5) Iron and Phosphorous. These subcycles are illustrated in figures 1through 5, respectively. Quaternary-averaged steady-state fluxes for these subcycles are tabulated in Appendix table A1; net atmospheric carbon dioxide and oxygen fluxes are provided in table A2.

To provide reasonable constraint on initial conditions for MAGic, we first calculated average Quaternary fluxes by assuming a balanced steady-state condition for all reservoirs. This assumption does not limit the overall operation of the model but merely provides a starting point and allows estimation of rate coefficients when appropriate. Values for reservoir masses and fluxes, reported in units of 1018 mols and 1018 mols/Ma, respectively, were obtained from searches of a large and diverse literature (for example, Garrels and Mackenzie, 1971; Garrels and Perry, 1974; Meybeck, 1979; Ronov, 1982; Berner and others, 1983; Holland, 1984; Lasaga and others, 1985; Gregor and others, 1988; Mackenzie, 1992; Compton and others, 2000). The initial model values are based on pre-anthropogenic estimates. These estimates are consistent with the composition of the dissolved load and the mineral and rock sources of individual dissolved components of rivers, the composition of materials cycled from sea to land through the atmosphere, the average and individual lithic compositions of sedimentary rocks, modern marine sediment compositions, and the composition of the geologically more recent ocean and atmosphere. Values consistent with these estimates were also adopted for the existing mass of sedimentary rocks, of oceanic salts, and of the total annual flux of elements through the oceanic environment. In order to close the cycles, best estimates were made of submarine basalt-seawater exchange fluxes and of fluxes involving diagenesis and alteration (so-called "reverse weathering" reactions, denoted RW in tables 2 and A2) of buried sediments (Martin and Sayles, 1994). Many of the assumptions made in the original quantitative model of the sedimentary rock cycle by Garrels and Mackenzie (1972), as well as those made in BLAG (Berner and others, 1983; Lasaga and others, 1985) and GEOCARB (Berner 1991, 1994, 1998, 1999; Berner and Kothavala, 2001), are valid for our model of the sedimentary rock cycle as linked to cycling processes involving submarine basalt alteration and subduction.

Following the approach given in BLAG, GEOCARB, and other related models, many of the fluxes are estimated by first computing a first order coefficient using average Quaternary (t = 0) values for the flux f and mass M:

equ4 with the resulting flux computed as a function of time by

mag5 These coefficients may themselves be held constant (k = k0) or varied with time (k = k(t)). In the latter case, this variation is accommodated by application of one or more dimensionless parameters, for example,


These time-dependent parameters are used to amplify or attenuate various processes related to land mass exposure, CO2 weathering feedbacks, seafloor spreading rate, sedimentary basin partitioning, et cetera, and are described where necessary in the sections below and summarized in table 4.

TABLE 4 Time-dependent dimensionless coefficients used in weathering and related equations (cBU computed from Budyko and others, 1987; others adapted from Berner, Lasaga, and Garrels, 1983; Berner and Kothavala, 2001)mag10


Fig. 1. Calcium-Magnesium-Silicate-Carbonate-CO2 Subcycle. Numbered labels refer to either reservoir index (identified by larger bold italic font) or flux indices (denoted by f [i], see table A1). The sign of the flux (for example, –f [i]) indicates the direction of the flow with respect to the arrow, and is positive unless so indicated; where fluxes are not labeled due to space constraints, a bold line indicates a bidirectional flow. Weathering of primary continental silicates, represented as normative anorthite (6) and Mg-pyroxene (7), and shelf carbonates calcite (8) and dolomite (9), occurs through attack by carbonic and sulfuric acids (see Sulfur Subcycle in fig. 4), and results in riverine fluxes of dissolved calcium (f [6], f [10], f [12] –f [15]), magnesium (f [7], f [11], f [13], f [15]), and associated bicarbonate. Dolomite deposition in shallow marine sediments is represented by a saturation state-driven flux (f [75]) and by a flux proportional to microbial sulfate reduction (f [71]). Calcite depositional fluxes are computed to be consistent with fluid inclusion data (Mg/Ca), and are associated with either shallow marine (f [74] + f [72], the latter driven by sulfate reduction) and pelagic (f [76], 35) reservoirs. Partitioning between shelf and pelagic calcite is governed by a time-dependent ß parameter. Shelf calcite and dolomite are then either re-weathered or enter a deeper burial regime (f [39], f [23]). Basalt-seawater reactions consist of magnesium uptake (f [81]) and net calcium release (f [96]); some calcium is also consumed in the formation of basalt calcite (f [94], 29); this carbon is eventually consumed during subduction (f [98]), returning CO2 to the atmosphere. Diagenetic alteration of shelf (f [42]) and pelagic (f [111]) sediments results in calcium and bicarbonate fluxes to seawater, a "reverse weathering" magnesium flux to chlorite (f [79]), and CO2 release (f [80]). Dolomite formed in the cratonic burial regime (19) consumes burial calcite (18), but is limited by magnesium supplied by the chlorite reservoir (f [45], 40), and is eventually returned by uplift to the continental weathering regime (f [44]). Burial calcite not consumed by dolomitization is either uplifted (f [47]), or eventually returns CO2 to the atmosphere and calcium to the silicate reservoir by metamorphic decarbonation (f [46]). Pelagic calcite is consumed through subduction (f [103]) and metamorphism (f [110]), with the partitioning between these targets controlled by a (?) parameter. Steady state fluxes are also listed in table A1. Not all fluxes (for example, CO2 exchanges associated precipitation and weathering) are shown.mag7


Fig. 2. Sodium-Potassium-Silicate-Carbonate-Chloride-CO2 Subcycle. See figure 1 caption for description of index notation of reservoirs and fluxes. Bold lines indicate bidirectional flows. Sodium and potassium are derived from weathering of Na-silicates (f [4], f [8], 4), K-silicates (f [5], f [9], 5), NaCl (f [26], 10) and KCl (f [27], 11), and balanced by bicarbonate, pyrite-derived sulfate, and chloride (f [35] = f [4] + f [8] + f [26], f [36] = f [5] + f [9] + f [27]). In seawater, dissolved Na+ (20) and K+ (22) are either re-deposited as evaporites via simple precipitation reactions (f [66], f [67]), or abstracted by seafloor basalt to form altered silicates (f [83], f [84]). Seawater K+ is also taken up by detrital silicates in reverse weathering reactions (f [77]), where it is exchanged for Na+ (f [78]) as a flux of Na-montmorillonite (f [123], 38) to form K-illite (39). Alteration products formed from reverse weathering and basalt uptake are ultimately returned to the craton (f [117], f [118], f [124]), where they are reincorporated in primary silicates to start a new weathering cycle.mag9


Continental Weathering Fluxes and CO2

Weathering fluxes are computed using equations similar to those in GEOCARB III (Berner and Kothavala, 2001), that is, by evaluating a feedback function sensitive to CO2 and additional parameters reflecting the role of vascular land plant fertilization, runoff, paleogeography, and differential susceptibility of carbonate and silicate rock types to weathering. In addition, weathering rate constants were also calibrated (eq 3) by additional time-dependent terms for land area, runoff, relief, and land plants. These parameters are adapted directly from GEOCARB III, which in turn were derived from estimates of these terms in the geological literature. Discussion of these terms is available in the various GEOCARB papers (for example, Berner and others, 1983; Berner and Kothavala, 2001) and is not repeated here.

Diagenetic and Reverse Weathering Fluxes

Reverse weathering as a sink for oceanic elements was proposed more than 40 years ago (Garrels, 1965; Mackenzie and Garrels, 1966a, 1996b). According to this hypothesis, the formation of clay minerals in marine sediments leads to the removal of soluble cations such as Mg2+ and K+, and the discharge of acidity (CO2), the opposite pathway of the weathering reaction. Conclusive field tests of this hypothesis have been difficult because the quantity of authigenic minerals that form is small and the neoformed phases are generally found in the very fine-grained sizes of the sediment. Thus, it is extremely difficult to detect these phases in marine sediments without time-consuming and tedious sediment size separations and detailed determination of their mineralogy and chemistry. Early studies of nearshore muddy sediments did not confirm that reverse weathering reactions occur (Russell, 1970; Drever, 1971), possibly because the studies were insufficiently detailed to permit detection of the small quantities of neoformed minerals present in the fine-grained fractions of the sediment. More recently, work has demonstrated the rapid formation of aluminosilicates in marine sediments and the possibility that the formation of these phases is a potentially important sink of oceanic K and Mg (Mackenzie and others, 1981; Michalopoulos and Aller, 1995). In addition, early diagenetic reactions involving dissolved aluminum in sediments (for example, Stoffyn-Egli, 1982) suggest the involvement of aluminum, silicon, and cations in the formation of authigenic aluminosilicate phases in modern marine sediments.

Another diagenetic reaction of interest is the uptake of Mg2+ in clay minerals during halmyrolysis and diagenesis in marine sediments that is balanced by the release of Ca2+ (Martin and Sayles, 1994), subsequently precipitated as CaCO3. This exchange reaction does not affect the CO2 balance of the atmosphere-ocean system, as do reverse weathering reactions.

Wollast and Mackenzie (1983) found that the balance between dissolved constituents entering the ocean via rivers and submarine basalt-seawater reactions and oceanic sinks, as well as that for particulate materials, required uptake of Mg2+ and K+ in exchange and reverse weathering reactions involving clay minerals. They suggested that these reactions might be responsible for the removal of as much as 50 percent of the weathering inputs of Mg2+, and 50 percent of the river input of K+ from sediment-porewater systems. These findings lead to the conclusion that exchange and reverse weathering reactions involving clay minerals are viable processes in terms of cation removal, at least for Mg2+ and K+, and should be a part of any model that attempts to depict the behavior of the sedimentary cycle over geologic time.

Fig. 3. Organic Carbon Subcycle. See figure 1 caption for description of index notation of reservoirs and fluxes. Bold lines indicate bidirectional flows. Red and blue arrows indicate flows of O2 and CO2, respectively. Marine organic carbon production is governed by relationships similar to those employed in the model of Van Cappellen and Ingall (1996), in which primary productivity (f [1]) is controlled by the availability of organic phosphorous (see Iron and Phosphorous Subcycle, fig. 5). A fraction (f [54]) of the marine productivity flux is then available for organic carbon sedimentation; the remainder is respired and returned as CO2 (f [3]). Organic carbon sedimentation is then partitioned into shallow marine (f [56], 3) and pelagic environments (f [57], 34). Shallow marine organic matter is then respired by two routes. Anaerobic oxidation via microbial sulfate reduction and pyritization (f [73] = (beta) f [50], where f [50] is sulfate reduction flux, see fig. 4) results in carbonate mineral formation (carbon flux to dolomite = 2 x f [71] and calcite = f [72]) and CO2 release (f [70]). Aerobic oxidation returns the remainder of shelf organic carbon as CO2 during organic matter weathering (f [20]). Pelagic organic matter is subducted and also returned as CO2 (f [108]). As in Berner and Canfield (1989), a terrestrial organic matter burial flux (f [21]) is computed from a record of sedimentary organic carbon (constraint curve cBU). Subsequent uplift and oxidation flux of this material (f [48]) is first order with respect to reservoir mass (17). A gas hydrate shunt (f [22], f [49]) has not been implemented in the current version of the model.mag11


Fig. 4. Sulfur Subcycle. Labeling scheme is similar to figure 1. Grayed portions of the figure are examined in more detail in figures 5 and 15. Sulfur is delivered to the ocean as sulfate derived from the weathering of pyrite (f [16], 14) and dissolution of evaporite minerals CaSO4 (f [24], 12) and MgSO4 (f [25], 13). Seawater sulfate is consumed by precipitation of evaporites (f [68], f [69],), and during sulfate reduction and pyritization, both within shallow marine sediment (f [50]) and during hydrothermal reaction of seawater and basalt (f [82], 30). The rate of microbial sulfate reduction is tied to the availability of reactive (ferric) iron (see fig. 5). Sulfide formed by sulfate reduction is then returned and made available for oxidation and weathering. Oxidation of this sulfide (f [16]) consumes the oxygen initially liberated by organic matter production, and yields sulfuric acid for weathering. The reduction of sulfate to basalt pyrite also consumes ferrous iron (silicate) and produces ferric iron; subduction of this ferric iron yields a return oxygen flux (f [119]), explained in more detail in the text (see also figs. 5 and 15).mag12


Fig. 5. Iron and Phosphorous Subcycle. Labeling scheme is similar to figure 1. Iron fluxes are shown in red, phosphorous in blue; P complexed with ferric iron is shown as magenta. The cycling of iron and phosphorous assimilates the approach of Van Cappellen and Ingall (1996), in which phosphorous availability exerts a fundamental control on long term ecosystem productivity. Phosphorous is delivered to seawater by apatite weathering (f [19]). The rate of organic phosphorous fixation, f [55], determines gross productivity (f [1], see fig. 3). The fixation rate is a first order function of the reactive phosphorous pool (27) and a constant parameter describing oceanic ventilation rate (VMIX). The burial fluxes of organic phosphorous in shelf (f [58]) and pelagic (f [59]) regimes are computed from the respective organic carbon burial fluxes (f [56] and f [57]; see fig. 3) and the C/P ratio. The C/P ratio is calculated using the DOA function (degree of anoxicity), in turn dependent on atmospheric oxygen mass (2) and f [1] (see text for details). The fraction of organic phosphorous that is not buried is returned to the reactive phosphorous pool (f [61]). In addition to the burial flux of phosphorous incorporated in organic matter, reactive phosphorous is also consumed during basalt uptake (f [53]), sedimentary apatite precipitation (shelf f [64], pelagic f [65], both tied to f [61]), and by binding to sedimented ferric iron as Fe-P complexes (shelf f [62], pelagic f [63]). All phosphorous is ultimately recycled to the continents as apatite: P in subducted pelagic sediments and basalt is returned via a mantle flux (f [113]).
The link between ferric iron and phosphorous provides critical buffering of atmospheric oxygen. The ferric iron flux (f [17], which includes iron originally derived from pyrite oxidation on land) is the source of iron for subsequent pyrite formation and that bound as Fe-P complexes. Thus iron consumed in the formation of pyrite during microbial sulfate reduction (f [50]) comes at the expense of ferric iron available for complexation with phosphorous (this remainder = f [92]) in shelf (f [52]) and in pelagic sediments (f [51]). Moreover, this partitioning is in turn a function of DOA: as DOA tends toward unity (complete anoxia), the ferric iron fraction decreases to zero. This feedback thus returns reactive phosphorous to the water column and makes it available for fixation, ultimately increasing O2 production, and restoring oxic conditions. Ferrous iron from iron silicate is used as the electron donor during hydrothermal sulfate reduction (f [82], see fig. 4). The oxide produced, together with pelagic ferric iron and basalt pyrite, is eventually subducted (f [130] + f [105] + f [99]). Pyrite iron is returned to the weathering cycle (f [115], fig. 4). The balance of the iron is returned to the basalt as reduced iron silicate (f [128]). Because iron is returned to the weathering cycle in a reduced state, the oxygen that will eventually be consumed must be supplied by a balancing flux (f [119]), possibly through consumption of a more reduced (and undefined) mantle component. Without this balancing flux, the mantle will act as an oxide sink.mag13


TABLE 4 Time-dependent dimensionless coefficients used in weathering and related equations (cBU computed from Budyko and others, 1987; others adapted from Berner, Lasaga, and Garrels, 1983; Berner and Kothavala, 2001)mag14

Seawater Saturation State, Model Operation, and Notation

Because of the potential for significant variation in composition over geologic time, seawater cannot be regarded as a constant ionic medium, and thus the expressions (apparent constants) commonly used to compute CO2-carbonic acid-carbonate equilibria in seawater (for example, Dickson and Millero, 1987) are not appropriate. We have instead adopted Pitzer’s approach in computing thermodynamic properties of seawater (Pitzer, 1973), extended with the parameters supplied by He and Morse (1993) for the carbonic acid system. At each time step, pH and carbonate mineral saturation state (omega = eG/RT = IAP/K) are computed from temperature, PCO2, the current vector of concentrations of Ca2+, Mg2+, K+, Na+, Cl-, SO42-, and total dissolved C (CO2 (aq)+ HCO3- + CO32-). The ocean is thus considered to be perfectly mixed and in equilibrium with atmospheric CO2, a valid assumption given the rapid mixing time of the ocean relative to the long time scale of our calculations. However, direct calculation of seawater activities introduces considerable complexity into the model, primarily in terms of potential depletion of the atmospheric CO2 reservoir and valid evaluation of hydrogen ion activity by the iterative Pitzer routines at each time step. Because of this complexity, we discuss the calculation approach in more detail. The model operation consists of the following sequence:

  1. Initialization and preliminary spin-up:
    1. Read initial (t = t0) reservoir masses y(t0), iteration control parameters, runtime options.
    2. Compute steady state values for rate coefficients (k0).
    3. Compute Pitzer parameters for dissolved seawater components.
    4. Compute time-dependent parameters and feedback functions.
    5. Initialize rate coefficients (k).
    6. Compute balance on function derivatives y'(t0).
    7. Hold t-dependent functions constant and iterate until (sigma)i(yi' / yi) < 2% for atmosphere and seawater components.

  2. Iteration from time t = t0 to t = 0. At each time step ti+1 = ti + (delta)ti:
    1. Update time-dependent parameters, feedback functions, and rate coefficients.
    2. Compute saturation state from seawater composition and current PCO2.
    3. Compute precipitation rates, weathering and exchange fluxes, fixation rates.
    4. Correct fluxes to be consistent with seawater ratio constraints.
    5. Take tentative step on y and evaluate error: If error is within limits, update y, store results and begin new time step, otherwise revise step size (delta)ti and repeat.

Thus during model operation, the step size (delta)ti) is not fixed but varies according to the difficulty and constraints imposed by the equations described below. All time-dependent parameters (ci) are available as smooth spline fits to the original data.
The general model, written in double precision FORTRAN 90, consists of forty coupled nonlinear ordinary differential equations, as well as over 100 fluxes, time-dependent rate constants, and associated parameters. It is solved as an explicit initial value problem using an Adams-Moulton predictor-corrector algorithm (DIVPAG, IMSL FORTRAN Numerical Library version 5.0, Visual Numerics, Inc.). With full compiler optimization and parallelization (INTEL VISUAL FORTRAN IA32), the model requires ~ 5 minutes for the complete 500 Ma run on a 2.1 GHz personal computer (dual AMD Athlon MP processors). To avoid confusion, model time (t) is negative in the past and has implicit units of Ma, thus t0 = –500 Ma, and the time of Quaternary completion indicated by t = 0. The 40 differential equations and constituent fluxes that compose the overall model, organized by reservoir, are given in Appendix table A3.
The notation scheme in this paper was designed to minimize the confusion in referring to the large number of reservoirs, components, fluxes, and associated parameters. The mass of the ith reservoir component is denoted y[i], for which steady state (Quaternary) masses in absolute units are listed in table A3; in the text these masses are usually expressed relative to these Quaternary values (that is, y[2]|t = 0 = 1). Fluxes f and their associated coefficients k are indexed in square brackets (f [1], f [2], ... k [1], k [2], and so on) according to their approximate order of appearance in the program. Fluxes not described in detail in the text appear instead in the figures and Appendix tables A1 and A2.

Calcium-Magnesium-Silicate-Carbonate-CO2 Subcycle

The geochemical cycle of calcium and magnesium (fig. 1, table A1) is the most important in terms of its interaction with carbon and the long-term maintenance of CO2. It is also the most complex. Calcium and magnesium are initially weathered from continental silicates and carbonates by carbonic and sulfuric acids (fluxes f [6], f [7], f [10], f [11], f [12], f [13], f [14], and f [15], fig. 1, table A1; Ca-Mg evaporite weathering is described in the sulfur subcycle). These fluxes are computed as a function of time following the template of equations 1 –3:

Ca silicate:

Formula 4 (4)
Formula 5 (5)

Mg Silicate:
Formula 6 (6)
Formula 7 (7)

Formula 8 (8)
Formula 9 (9)

Formula 10 (10)
Formula 11 (11)
Formula 12 (12)
Formula 13 (13)
Formula 14 (14)
Formula 15 (15)

and (gamma) [i] is the fraction of the total pyrite weathering flux (f [16]) for each rock type, and y [6], y [7], y [8], y [9] are the masses of Ca- and Mg-silicate, calcite, and dolomite, respectively (see table A3). As an aside, the actions of GEOCARB dimensionless parameters in equations 12 –15, reflecting the time-dependent forcing imposed by land area/runoff, relief, CO2 feedbacks, land plants, and palegeography, exist for most of the coefficients (k) related to weathering presented in the remaining descriptions. In these cases, weathering parameters are used without change from Berner and Kothavala (2001), and thus will not appear explicitly except as warranted for clarity (see table 4). The trend of key forcing functions is shown in figure 6.

Once present in seawater, these components are redeposited in shelf environments as primary or secondary calcite (f [74], f [72]), dolomite (f [75], f [71]) and CaSO4 and MgSO4 evaporites (f [68], f [69]). Subsequent to deposition, they either reenter the weathering cycle through uplift and exposure (f [12], f [13]), (f [14], and f [15]) or are removed to deeper burial regimes (f [39], f [23]). Pelagic environments receive only calcium carbonate (f [74]), with the total CaCO3 flux (f [93]) partitioned between shelf and pelagic regimes by a time-dependent (ß) parameter:



We should note here that our pelagic reservoir includes slope and rise sediments as well. The value of ß is varied to produce the observed increase in carbonate pelagic deposition during the Mesozoic (Garrels and Mackenzie, 1971; Wilkinson and Walker, 1989). At each time step, the total precipitation flux of calcite (f [93]) is calculated from the current values of seawater Mg (y[24]) and Ca masses (y[23]), their respective derivatives, and other fluxes involving calcium and magnesium (see table A3):


where the derivative of the cation ratio (R'(t)) is computed from the fluid inclusion data of Lowenstein and others (2001) and Hardie (1996).

Shelf dolomite precipitation rate (f [75]) is controlled by seawater saturation state; dolomite is also formed in the same environment during the oxidation of organic matter ("organogenic" dolomitization, f [71]), tied to the rate of microbial sulfate reduction (f [50]):



In equation (19), (f [71]/f [50])|t = 0 is the ratio of steady state Quaternary (pre-glacial) fluxes. In equation (20), c'Ap(t) is a function describing the shelf area available for deposition of carbonate sediment. This function is estimated from Quaternary deposition fluxes in shallow water environments and the submerged shelf area as a function of time, and is normalized to its Quaternary value (thus c'Ap(0) = 1). Similarly, r(omega) = kT (omega)dol –1)n, where (omega)dol is the saturation state of seawater with respect to dolomite, kT and n are the rate constant and reaction order taken from experimental work (Arvidson and Mackenzie, 1999), and (f [75] / r (omega))|t = 0 is the ratio evaluated at their Quaternary steady state values.

After deposition, pelagic calcite is either subducted to the mantle (f [103]) or removed to a metamorphic regime (f [110]); partitioning of the total pelagic flux (f [112]) between these two paths is controlled by another time-dependent (sigma) parameter (Volk, 1989):



Both of these pathways lead to decarbonation reactions that serve to return material to the weathering cycle. We should note here that in the model the mantle does not constitute a reservoir per se, but simply a collection point for return fluxes to either the silicate (basaltic or continental crust) or atmospheric reservoirs. This simplification, although significant, was made to avoid the necessity of estimating mantle abundances and time-dependent fluxes that are not well known and would introduce additional uncertainty in the model. Basalt-seawater reactions result in magnesium uptake (f [81]) and calcium release (f [96]). Magnesium uptake is assumed to be first order with respect to the Mg2+ concentration in seawater, y[16]; the rate coefficient for this uptake is modified by the time-dependent parameter reflecting seafloor spreading rate, cG, that is,


(see table 4). The cG parameter is identical to the fG parameter in GEOCARB III (Berner and Kothavala, 2001), using the curve concatenated from the data of Engebretson and others (1992) and Gaffin (1987).

The complementary release of Ca2+ from basalt during seawater-basalt interaction (f [96]), although identical to the magnesium uptake flux (f [81]) in the steady state (Quaternary) condition, is allowed to otherwise vary to balance with the uptake fluxes, a description of which we will postpone to our later discussion of the standard run. In addition, diagenetic alteration of shelf (f [42]) and pelagic (f [111]) sediments also results in release of calcium to seawater that is first order with respect to sediment masses (y[3], y[4]):



Magnesium uptake (f [79]) by the sediment to form chlorite, a mineral composition representative of neoformed magnesium-rich silicates, is allowed to vary according to seawater Mg2+ mass (y[16]):


The chlorite thus formed reacts with buried calcite to form burial dolomite. The rate of burial dolomite formation (f(45)) is thus limited by the size of the chlorite reservoir (y[35]):


This dolomite is eventually returned by uplift to the continental weathering regime, and burial calcite not consumed by dolomitization eventually returns CO2 to the atmosphere and calcium to the silicate reservoir by metamorphic decarbonation (f [46]).


Fig. 6. Time variation of forcing parameters (see also table 4). GEOG = cGEOG, paleogeographic effect; Bs = cBs, CO2 weathering feedback (silicates); Bc = cBc, CO2 weathering feedback (carbonates); AD = cAD, area x runoff (carbonates); LA = cLA, carbonate land area; R = cR, uplift/relief; E = cE, terrestrial plant weathering; G = cG, degassing rate; Ap = cAp, shelf depositional area; D = cD, runoff; Aw = cAw, weathering area. Data are from Berner and Kothavala (2001) and Berner and others (1983). ß is shelf-pelagic partitioning factor (ß = 0 constrains all precipitation to occur on the shelf). (Beta) is partitioning the fate of pelagic carbonate sediment to subduction versus metamorphic regimes. Bs and Bc trends are taken from standard run (see fig. 7).


Sodium-Potassium-Silicate-Carbonate-Chloride-CO2 Subcycle

Sodium and potassium river fluxes are derived from the weathering of Na- and K-silicate masses (y[4], y[5]) by both CO2 (f [4], f [5]) and H2SO4 derived from oxidation of FeS2 (f [8], f [9]), and NaCl and KCl evaporites (f [26], f [27]), and are thus balanced by bicarbonate, pyrite-derived sulfate, and chloride (fig. 2 and table A1):

Na silicate:

Formula 28 (28)
Formula 29 (29)

K silicate:
Formula 30 (30)
Formula 31 (31)

In basalt-seawater reactions, the basalt uptake fluxes of Na+, K+, and SO42- (f [83], f [84], and f [82]) are first order with respect to their seawater masses (y[13], y[14], and y[26]); A basalt calcite flux (f [94]) accounts for bicarbonate uptake and CO2 release (f [85]), described in detail in the MODEL RESULTS AND DISCUSSION section.

Formula 32 (32)
Formula 33 (33)
Formula 34 (34)
Formula 35

At steady state, these sulfate and bicarbonate fluxes balance the H2SO4-derived pyrite and CO2 involved in the original weathering reactions. The bicarbonate reacts with Ca minerals within the basalt to form calcite (see Alt and Teagle, 1999; some calcium is derived from seawater itself but the amount is debatable) that is ultimately returned to the cratonic weathering cycle (f [98], fig. 2). Na+ and K+ taken up within the basalt are used to form silicate alteration products, and these are ultimately subducted (f [101], f [102]) and returned to the craton (f [117], f [118]). Sulfate is reduced during oxidization of iron silicates in the basalt to form pyrite. We emphasize here that this is not the pyrite formed by low-temperature, diagenetic reactions in marine sediments in the model. The hematite or magnetite, produced by the reaction involving reduction of sulfate in the basalt, must itself be subsequently reduced after subduction of the oceanic crust to liberate the free oxygen necessary to weather pyrite returned to the craton (see table 3 and fig. 4). If this is not done, the mantle will act as an oxygen sink (Lécuyer and Ricard, 1999). The implications of this treatment of cation and anion seawater-basalt fluxes are also discussed below.

Steady state chloride balance is preserved through simple precipitation of NaCl and KCl (f [66], f [67])—no volcanic HCl component is currently included. These fluxes are varied during model operation in first order proportion to their respective ion activity products. In addition, K+ is also exchanged for Na+ during diagenetic reactions with Na-montmorillonite and formation of K-illite (f [77]), a reaction that also consumes bicarbonate and releases acid as CO2 (f [78]). Steady state requirements also dictate that these Na- and K-phases ultimately be returned to the craton (f [127], f [124]), where they are reincorporated in primary silicates to start a new weathering cycle.

Organic Carbon Subcycle

Organic matter cycling exerts a critical control on both long-term atmospheric CO2 and O2, and our model reflects our belief that net ecosystem production is ultimately limited by phosphorous availability (Smith and Mackenzie, 1987; Tyrrell, 1999). The controls on organic matter burial are also complex, and are described in the Iron and Phosphorous Subcycle (see below). Although the availability of organic matter in marine sediment also reflects the total sedimentation rate, primarily inputs of clastic detritus from terrestrial weathering (Canfield, 1991; Berner, 2004; Mackenzie and Lerman, 2006), the current model contains neither clastic weathering nor organic carbon inputs from land to the coastal ocean, and the primary influence of terrestrial weathering on the marine organic carbon cycle is through the supply of inorganic phosphorous (the apatite weathering flux, f [19]). The control of the rate of marine organic matter fixation (f [1]) closely follows the model of Van Cappellen and Ingall (1996), in which the rate is proportional (via (gamma)RR, the Redfield ratio) to f [55], the rate of organic phosphorous uptake:


The organic phosphorous flux is in turn a first order function of the mass of reactive inorganic phosphorous in the ocean (y[28]) and a rate function (VMIX) that describes the ocean’s ventilation rate:


In all the runs described below, VMIX is fixed to a value of 3.0. In the steady state model, organic phosphorous fluxes are computed by dividing Quaternary organic carbon fluxes by 250 (table A1). This approach adopts Van Cappellen and Ingall’s (1996) argument that C/P ratios vary from 200 for fully oxic to 4000 for completely anoxic sediments. The flux of organic material available for marine deposition, f [54], is computed as:


where n = 2.5 (Van Cappellen and Ingall, 1996), and the respiration flux (f [3]) in turn is computed as the difference between (f [1]) and (f [54]). The net burial of organic material on the shelf (f [56]) reflects the difference between f [54], the flux of organic carbon oxidized in the formation of pyrite (viz. f [50]), in turn tied to reactive iron availability (described below), and the pelagic organic carbon flux, f [57]. The partitioning of organic matter burial between shelf and pelagic regimes is currently fixed so that the shelf receives 95 percent of the total. Pelagic sedimentary organic carbon is eventually subducted and returned to the atmosphere (f [108], (f [114]).

The mass of organic carbon residing in shelf sediments (y[6]) may be rapidly recycled and returned as CO2 (f [20]):


In addition, we have provided for the development of burial of terrestrial organic matter in non marine basins (f [21]) using an approach similar to that of Berner and Canfield (1989). We have used data derived from the organic carbon content of the rock record (Budyko and others, 1987) to calculate the total derivative (cBU(t)) of this reservoir. This constraint and the assumption of a first order uplift-weathering flux (f [48] = k[23] y[17]) yield the required burial flux:


The cycling of organic carbon, coupled with sulfur and iron, results in large exchange fluxes involving atmospheric reservoirs of CO2 and O2 (table A2; eq A2, table A3). For purposes of clarity, we note that the balance of atmospheric oxygen reflects reactions in which it may not be directly involved as an oxidant. During photosynthesis, electrons lost during ATP and NADPH production are returned to chlorophyll via oxidation of water, yielding free oxygen; the chemical energy stored in ATP and NADPH is then used to convert CO2 and water to sugars and carbohydrates (carbon fixation). A fraction of this organic carbon, after some transit through the food web, may ultimately be deposited as sedimentary organic matter. Subsequent consumption of this organic matter by aerobes will consume free oxygen, thus balancing the cycle. However, in anaerobic marine environments, this organic matter is commonly oxidized by bacterial sulfate reduction to yield dissolved H2S. Over a complex subsequent reaction path(s) this dissolved sulfide will ultimately react to form pyrite, FeS2. While the weathering of pyrite on land will consume oxygen, burial of pyrite actually represents oxygen accumulation in the atmosphere. However, this oxygen was originally produced during photosynthesis. This can be seen by recognizing that the sum of the overall reaction for sulfate reduction (eq 76' is eq 76 (table 2) written on a per mol carbon basis without calcium),


and pyrite oxidation (eq 78, table 2),


is oxidation of organic matter by oxygen:


Thus, there is no explicit term of pyrite burial in the oxygen balance (table A2 and eq A2), but the role of FeS2 weathering in consuming oxygen and its formation representing oxygen accumulation in the atmosphere are an integral part of MAGic.


Sulfur Subcycle

Sulfur is delivered to the ocean as sulfate derived from the weathering of pyrite:


and the first order dissolution of evaporite minerals CaSO4 and MgSO4. Although the sulfuric acid released from pyrite oxidation participates in silicate and carbonate rock weathering reactions (eqs 5, 7, 9, 11, 29, and 31), the fluxes involved are minor compared to CO2, and have no associated feedback functions. Sulfate reduction occurs both biogenically within shallow marine sediment and abiogenically during hydrothermal reaction between seawater and basalt, as mentioned above. The bicarbonate released in biogenic sulfate reduction is used to form calcite and dolomite. Reduced sulfur formed by these pathways is returned through uplift and exposure of marine sediments and through release of volcanic H2S and SO2. Oxidation of this reduced sulfur thus completes the cycle, consuming the oxygen initially liberated by organic matter production, and yielding sulfuric acid for weathering.


Iron and Phosphorous Subcycle

Iron and phosphorous cycling are closely linked to that of organic carbon and sulfur. Sedimentary phosphorous is represented as organic P, P bound with Fe oxyhydroxides, and apatite phases, and partitioned between shelf and pelagic environments. The rate of incorporation of P in organic matter, f [55], is allowed to vary in proportion to the availability of so-called "reactive" phosphorous. The overall approach closely follows the model of Van Cappellen and Ingall (1996), which uses the ocean’s overall oxidation state or "degree of anoxicity" (DOA) to represent the combined influence of the ocean’s oxygen tension, primary production flux, and ventilation rate. DOA effectively parameterizes the coupling of iron and phosphorous governing the efficiency of sedimentary phosphorous release and recycling, the fraction of reactive iron present as reduced (pyrite) versus oxyhydroxide phases, and the C/P ratio of organic matter undergoing burial. This coupling yields a critical negative feedback control on atmospheric oxygen levels. If PO2 becomes too high (as DOA approaches zero), then the phosphorous sorbed to reactive FeOOH undergoing deposition and burial will remain within the sediment, and its removal from seawater will limit productivity and allow PO2 to fall. Conversely, if PO2 becomes too low, for example, due to enhanced demand in the weathering of shelf organic matter, the increase in anoxia (as DOA approaches unity) will lead to a greater fraction of reactive iron being sulfidized. Loss of previously FeOOH-bound P to the water column will in turn increase productivity and PO2.



Example Paper:

MAGic: A Phanerozoic Model for the Geochemical Cycling of Major Rock-Forming Components

Rolf S. Arvidson*, Fred T. Mackenzie** and Michael Guidry**
* Department of Earth Science MS-126, Rice University, P.O. Box 1892, Houston, Texas 77251-1892 USA; Rolf.S.Arvidson@rice.edu
** Department of Oceanography, School of Ocean and Earth Science and Technology, University of Hawaii at Manoa, Honolulu, Hawaii 96822 USA


A dynamical model (MAGic) is presented that describes the elemental cycling of sedimentary materials involving sodium, potassium, calcium, magnesium, chloride, carbon, oxygen, iron, sulfur and phosphorous through much of the Phanerozoic. The model incorporates the basic reactions controlling atmospheric carbon dioxide and oxygen concentrations, continental and seafloor weathering of silicate and carbonate rocks, net ecosystem productivity, basalt-seawater exchange reactions, precipitation and diagenesis of chemical sediments and authigenic silicates, oxidation-reduction reactions involving carbon, sulfur, and iron, and subduction-decarbonation reactions. Although MAGic contains feedback and forcing functions adapted from the GEOCARB models (Berner, 1991, 1994; Berner and Kothavala, 2001), these functions are incorporated in a reservoir-reaction scheme that is considerably more detailed. Coupled reservoirs include shallow and deep cratonic silicate and carbonate rocks and sediments, seawater, atmosphere, oceanic sediments and basalts, and the shallow mantle. Model results are reasonably consistent with recently published constraints provided by fluid inclusion, isotopic, floral, and mineralogical records. We have used these results to evaluate sensitivity to uncertainties in the history of the earth-ocean-atmosphere system over the past 500 Ma: the advent of pelagic carbonate sedimentation, the importance of burial versus early diagenetic dolomite formation, the importance of reverse weathering, and the relationship of these processes to seafloor spreading rates. Results include a general pattern of dolomite abundance during periods of elevated seafloor spreading and alkalinity production, elevated atmospheric CO2 concentrations for most of the Phanerozoic similar to those predicted by GEOCARB, and covariance of seawater sulfate to calcium ratios with magnesium to calcium ratios. These trends are broadly consistent with proxies for seawater composition and the mass-age data of the rock record itself.



[ << Back to Models]

Return to Fred Mackenzie's home page

[ Top of Page ]

SOEST Home | Directory | Oceanography | Research at the School


Last Revised: Feb 2008
by TKC