## Abstract

To investigate the weathering of sedimentary organic matter and its role in regulating atmospheric oxygen, a theoretical modeling study is presented that addresses the fundamental controls on atmospheric oxygen uptake: erosion rate, organic matter content, and reaction rate. We compare model results with the previous part of this study that analyzed a drill core of black shale from the New Albany formation (Upper Devonian, Clay City, KY) for total and organic carbon, pyrite sulfur, porosity, permeability and specific surface area. As was observed in the field study, the model predicts that the loss of organic matter by oxidative weathering takes place across a reaction “front” where organic carbon content decreases sharply toward the land surface along with pyrite loss.

The model is based on kinetic control of reaction of organic matter and pyrite with O_{2} dissolved in soil water. The downward diffusion of gaseous O_{2} partitions with dissolved O_{2} in water films on sediment grains via Henry’s law. Once a weathering profile is developed, the downward migrating O_{2} reacts with shale organic matter and pyrite. Pyrite reacts faster with O_{2} than does organic matter (for a given local concentration of oxygen) making the pyrite front generally deeper than the organic matter front. We explore the influence of differing erosion rates, atmospheric O_{2} concentrations, organic matter contents, porosities, tortuosities, and rates of reaction (that could include possible acceleration due to microbes) on the oxygen consumption.

We conclude, based on our modeling, that the erosion rate and the concentration of buried reduced matter, as opposed to the level of atmospheric O_{2}, normally limits the rate of drawdown of atmospheric oxygen. For the vast majority of erosion rates and Phanerozoic oxygen levels, essentially all ancient reduced material is oxidized before reaching the surface. Only in regions of unusually rapid erosion or during very low atmospheric oxygen levels can rates of diffusion of O_{2} in soils and rates of reaction control O_{2} drawdown, leading to weathering that is O_{2}-dependent. In this case erosion and rapid reburial of unoxidized organic matter would occur.

## INTRODUCTION

Major control of atmospheric oxygen levels over geologic time (millions of years) is exerted by the burial of organic matter in sediments and the weathering of ancient organic matter (**OM**) (Ebelmen, 1845; Garrels and Perry, 1974; Holland, 1978, 1984; Berner and others, 2003; Berner, 2004). These processes can be represented by the two overall reactions (with CH_{2}O representing a simplified stoichiometric formula for OM): The first reaction represents burial and the second represents weathering. Compared to reactions involving OM, the second most important control on long-term oxygen evolution in the atmosphere is the formation and destruction of sedimentary pyrite. One form of the pyrite oxidation reaction is: Pyrite forms in anoxic sediments accompanying organic matter burial and is oxidized along with organic matter when such sediments become exposed to the atmosphere.

This study is motivated by the desire to understand the weathering in more detail. Of central importance is whether or not in eroding environments the ancient OM is completely oxidized before it becomes exposed by erosion at the surface, or whether a finite quantity remains unoxidized at the surface where it can be stripped off by erosion, transported downstream and ultimately reburied as relic ancient OM. The response of the rate of O_{2} uptake as a function of atmospheric O_{2} level that would result from these two scenarios is quite different. Consider the case of constant uplift and physical erosion of sediments containing ancient OM. If oxygen levels were low enough, old OM might not be fully or completely oxidized before reaching the surface. Then a fraction of unweathered OM would enter the short time-scale riverine system and soon after become reburied. Under this scenario, there would be less drawdown of oxygen from the atmosphere than the case where all the OM were oxidized. This would imply that low oxygen levels lead to less oxygen drawdown, implying increasing oxygen, and eventual stabilization of oxygen in the atmosphere by negative feedback. Is the oxygen drawdown from the atmosphere dependent on the amount of buried OM being exposed by erosion *and* the oxygen level of the atmosphere (as concluded by Lasaga and Ohmoto, 2002), or is essentially all of the ancient OM oxidized before reaching the surface independent of the oxygen level in the atmosphere? We will conclude that for the vast majority of erosion rates and Phanerozoic oxygen levels that the drawdown of oxygen from the atmosphere is independent of the oxygen level. The broad conclusions of Lasaga and Ohmoto (2002) about stability of oxygen in the atmosphere are still valid, but the oxygen level feedback on weathering-related oxygen drawdown rates is less important than they have claimed.

This problem has been the focus of considerable discussion (Holland, 1978, 1984; Lasaga and Ohmoto, 2002; Holland, 2003; Ohmoto, 2003). The recent work of Lasaga and Ohmoto (2002) addresses this issue using a box model, where a box of a fixed size (see fig. 1A) represents the reaction zone. They also examined the influence of the box size using L_{s} = 1, 3, or 10 m. As stated in Lasaga and Ohmoto (2002):

“Normally chemical reactions lead to gradients in weathering with depth of the soils. In our treatment we will mix the soil so that chemical changes occur uniformly throughout the depth, L

_{s}, of the soil. This depth, L_{s}, marks the areas of high permeability which is most accessible to reaction.”

It is possible that this box model: 1) overestimates the level of oxygen in the reaction zone, 2) effectively brings fresh unreacted OM all the way to the surface due to the “well mixed” nature of the box, and 3) overestimates the evolving OM surface area (by keeping it fixed rather than letting it evolve during oxidation). Compared to our model, their approach apparently favors the erosion, transport, and reburial of OM resulting in the dependence of O_{2} uptake on the level of atmospheric O_{2}. A result from our boundary-layer resolving model, to be completely described below, is shown in figure 1B. A one-dimensional finite difference model was used, with enough grid points to resolve the boundary layer. The model includes oxygen diffusion from the surface downward to a reaction zone, and oxidation of ancient organic matter. The reference frame is fixed to the eroding surface where material is removed at the imposed erosion rate, and uplift brings fresh unoxidized OM into the bottom of the box (well below the reaction zone).

In modeling organic matter weathering we find the most serious issue in need of attention is the nature of the near-surface boundary layer. We find that the conclusions of box models can be erroneous, as they lack the required spatial resolution. Few studies (Petsch, ms, 2000) have attempted to resolve the boundary layer near Earth’s eroding surface that naturally results from the competing rates of: 1) oxygen diffusion into the soil, 2) erosion and re-supply of “fresh” unoxidized material from below, and 3) kinetically controlled oxidation of ancient OM and potentially pyrite. In the present paper we model oxygen diffusion, surface erosion, and kinetically controlled oxidation reactions of OM and pyrite to resolve such boundary layers. Our model assumes that the oxygen transport to the zone of reaction is by gas diffusion through the interconnected pore spaces, followed by Henry’s law partitioning into water films adsorbed on the solid particles and reaction of the dissolved O_{2} with organic matter and pyrite. Also, the influence of microbial acceleration of the kinetic rates can be easily addressed within the context of our model (for example, Daniels and others, 1996; Petsch and others 2001a).

In the sections below, we introduce the mathematical formulation of the model in detail, along with a discussion of the major assumptions. In part I of this study (Wildman and others, 2004), we analyzed a core sample of weathered black shale, and this will serve to guide some of the choices of parameters, although most of our introductory model runs are focused on shales with lower (and more typical) OM content (1% OM). We start the results and discussion section with a model run (the default run) to which many of the other simulations are compared. The typical characteristics of the time evolution and the steady state are described. This is followed by sensitivity studies that reveal how the model results depend on the most important parameters. The surface layer typically becomes devoid of OM due to oxidation. This leads to a buried steady-state transition zone (a “front”), with an OM-free zone near the surface transitioning over short distances to an OM rich zone at depth. The dependence of the front characteristics upon the model parameters is analyzed. In particular, we examine how the front depth and width depend on the rate of reaction, the air fraction in the pores, and the quantity of OM at depth that is being “brought to the surface” by erosion. We also examine the influence of the “grain size” of the OM, the tortuosity formulation, the erosion rate, and the atmospheric oxygen level. Specific comparisons are also made to the box model parameters of Lasaga and Ohmoto (2002). We examine the characteristic “time-to-steady-state” and the balance of fluxes of oxygen (downward flux) and OM (upward flux) at steady state for buried fronts, which could even aid in estimates of local erosion rates. A series of simulations are focused on cases that could bring finite amounts of OM to the surface, where low oxygen levels or rapid erosion rates would be required. Finally, we compare model runs with both OM and pyrite present to the results taken from a core sample of a weathered black shale from part I of this study (Wildman and others, 2004). This comparison is followed by a summary and conclusions. In Appendix A we describe various conversion factors useful to the model. In Appendix B we describe the numerical methods used.

This modeling study was motivated by the results of previous field studies (see Petsch and others, 2000; and Wildman and others, 2004 for a summary of earlier work). In these field studies, distinct patterns of OM loss and alteration are observed associated with weathering of ancient sedimentary rocks. During weathering, the majority of total OM is lost from the rock prior to erosion; radiocarbon analysis of residual OM reveals that a small fraction of ancient kerogen may remain at the surface of some weathering profiles, especially those characterized by high initial TOC contents and rapid rates of erosion. This view of limited export of relict OM is further supported by a growing number of studies that suggest transport of aged (kerogen-like) particulate OM in rivers and export to modern sediments (Kao and Lui, 1996; Leithold and Blair, 2001: Masiello and Druffel, 2001; Dickens and others, 2004; Komada and others, 2004, 2005; Leithold and others, 2005). However, these systems may reflect high erosion rates and rapid transport times typical of small mountainous watersheds, or drastically increased sediment yields associated with historical changes in land use. Thus, this study seeks to address the fate of weathered sedimentary OM in environments with more typical pre- agricultural, global average erosion rates. Additionally, this study improves upon the box model of Lasaga and Ohmoto (2002) that examined feedback between atmospheric oxygen levels and the oxidation of ancient OM, by generating a finite difference model that can resolve the near surface boundary layer of the weathering profile. This approach may then be employed to improve models of atmospheric oxygen (for example Bergman and others, 2004) that are based on the Lasaga and Ohmoto (2002) box model.

Our model is not valid for extremely wet or extremely dry environments. Arid environments may lack the water film that seems to be necessary for OM oxidation to occur. The model domain is assumed to be above the water table. Below the water table oxygen levels would generally be very low if OM were present. We have not included compaction in the model in a rigorous way, but the first-order effects of compaction to limit excessive increases in porosity (that could accompany oxidation of large quantities of OM and pyrite) are approximated by use of an upper limit “lid” on the air porosity. We have also focused on a hypothetical state, where oxygen diffuses through stagnant interconnected pores. We comment on the influence of heterogeneities and cycles of precipitation and drying that should be addressed in future more sophisticated models.

## THE MODEL

### Mathematical Formulation

We present a model of weathering of ancient organic matter (OM) and pyrite via oxidation. During oxidation, the grains of OM and pyrite reduce in size as a function of space and time, affecting their surface areas and the porosity in the “soil” column. Oxygen is assumed to diffuse downward through the interconnected pore spaces. We assume that the oxygen in the pore space air exchanges with oxygen in adsorbed water films on the grain surfaces via Henry’s law. From the water films, the dissolved oxygen can react with the OM or pyrite via kinetic control. In this simple model we ignore CO_{2} (that would primarily be diffusing upward) and assume the air in the pore spaces is essentially stagnant. Equations describing oxygen diffusion and reaction, as well as the uplift and reaction of the reduced solids are: (1) (2) for each *i*. The subscript or superscript *i* (ranging up to *i*max) indicates the grain type being oxidized. In most of our simulations we used *i* =1 for OM and *i* =2 for pyrite, but the model includes the possibility of multiple grains size bins, each with their own characteristic kinetics and initial sizes. We also define: (3) for Michaelis-Menten or Monod (Boudreau, 1997) kinetics, and (4) for power-law kinetics. Note that: (5) and that (6) for our grain model. The latter conversion factor *b _{i}* is discussed later in Appendix A: the kinetics are surface area controlled and this factor properly accounts for that fact.

Also in the above:

*a*- the oxygen concentration in the air-filled pores [moles O
_{2}/cm^{3}air], - g
_{i} - the concentration of solid undergoing oxidation [(moles C in OM or moles pyrite)/(cm
^{3}bulk)], *t*- time [yr],
*x*- downward distance [cm] with
*x*= 0 fixed at the eroding land surface, *D*_{s}- the effective diffusion parameter (in cm
^{2}/yr as written here) that includes a porosity dependent tortuosity (compare to Berner, 1980) with*D*=_{s}*T̄***D*. Here,_{0}*T̄** is a nondimensional tortuosity factor (as defined in Bear, 1972, and equivalent to θ^{−2}of Berner, 1980), and the diffusivity of oxygen in air is*D*(cm_{0}^{2}/s converted to cm^{2}/yr). We allowed the general form of Archie’s law, which using*n*= m*−1 (using*m*for Berner’s*n*of p. 36: Berner, 1980) becomes*T̄** =*b*(_{T}*φ*)_{a}. Here,^{n*}*b*an empirical constant, normally 1 as used here(except as noted). The most typical Archie’s law uses_{T}*m*= 2, although we find higher values are more consistent with tight shales. Additional issues related to the tortuosity are discussed below, in the section where parameters are chosen. *φ*_{a}- the volume fraction of air in the bulk (air-filled porosity), with
*φ*=_{a}*s*_{air}*φ*, where*s*is the air saturation (_{air}*s*= 1-_{air}*s*) and_{water}*φ*is the total porosity (volume of air and water filled pores, assumed to be interconnected, per bulk volume), *ω*- erosion rate (cm/yr),
*ω*<0, so that material is brought up from below (see also below), *ℜ*_{i}- reaction rate for reduced component
*i*(in units of mole pyrite or C in OM consumed per cm^{3}bulk per year), - Α
_{i} - is a unit conversion factor for power-law kinetics (see Appendix A),
*R*_{max}^{i}- maximum rate of reaction at high oxygen levels, a Michaelis-Menten parameter (mole C in OM or pyrite consumed/m
^{2}/yr), (Boudreau, 1997), *K*_{m}- the oxygen level in the air (mole O
_{2}/cm^{3}air) for half the maximum rate (another Michaelis-Menten parameter). This parameter is adjusted via a temperature dependent Henry’s law constant to relate experimental aqueous phase kinetic measurements as a function of oxygen to local pore-air concentrations of oxygen. This factor is described more below and in Appendix A. *υ*_{i}- a stoichiometric coefficient expressing the number of moles of oxygen consumed per mole of solid material being consumed (pyrite, or C in OM),
*c*_{i}- a conversion factor to adjust units (see Appendix A), and
*d*_{i}- the characteristic grain thickness (m) for OM or pyrite (more on shape factors below, and unit conversions are described in Appendix A).

We have used an effective diffusion coefficient based on air in interconnected pore spaces. Effective transport of oxygen depends critically on the interconnectivity of the air in the pore spaces. If the soil is completely water saturated, oxygen diffusion into the soil is severely limited, due to the low diffusivity of oxygen in water. Even moving water is not very effective in transporting oxygen downward, due to the low oxygen content of water in equilibrium with air. As our comparison core sample (from part I of this paper, Wildman and others, 2004) was extracted well above the water table, and as the water saturation in the shales is likely to be often lower than that which would prohibit the passage of oxygen through the gaseous phase, we limit our modeling to diffusion through the gas. We evaluated the importance of delays associated with diffusion through thin water films adsorbed to the solid particles and found this component of a multiphase diffusion model to be unimportant. From our examination of the literature, most correlations relating effective diffusion coefficients in porous solids to diffusion in pure fluids were measured in materials with porosities greater than 0.1, significantly larger than the porosity in the deeper section of the core analyzed in part I of this paper (Wildman and others, 2004). Even so, we have used extrapolations of tortuosity factors consistent with literature values.

### Choice of Parameters

Many of the parameters are discussed above. We have generally used 17°C for the modeling runs (relevant for the temperature dependence of the oxidation kinetics, the Henry law partitioning, and the diffusion constant), which is a mid-latitude value appropriate for comparison to the field site (Wildman and others, 2004) at Clay City, Kentucky (near 35N latitude, compare fig. 5 of Ekart and others 1999; Hu and Feng, 2003). The global mean effective temperature for weathering is likely to be slightly higher than our choice. Given the increase in reaction rate with temperature for both OM and pyrite reactions, our choice of model temperature is conservative, in the sense that higher temperatures favor less OM present at the eroding surface.

Shape factors were assumed for different grain types. In the model, porosity increased as OM and pyrite were consumed (as the volume fraction, surface areas, and grain thickness *d _{i}* of OM and pyrite are calculated during oxidation, the implied change in porosity is also calculated). In some models we limited the maximum allowed air porosity at a prescribed “lid” value, as a simple way of incorporating the most important effects of compaction that can serve to limit large increases in porosity. Increasing the porosity has the effect of increasing the effective diffusion coefficient because the porosity-dependent tortuosity factor increases (implying a less tortuous path). The pore sizes observed in the SEM images (part I of this study, Wildman and others, 2004) were up to 2 μm in width, but were generally less than that. Lacking a direct measurement of the tortuosity, we have generally used Archie’s law as discussed above. Our model, however, allows for a more general dependence of the effective diffusion coefficient upon the air-filled porosity (compare to Archie, 1942; Troeh and others, 1982; Shimamuri, 1992; Scanlon and others, 2002; Mbonimpa and others, 2003). Note that various authors differ in whether or not a porosity factor is included in the definition of the effective diffusion coefficient. Knudsen diffusion may play a role in increasing molecular diffusivity of oxygen in the smallest pores, but in our core samples the structure suggests such transport is less important than bulk molecular diffusion through the pores.

For OM, we assumed plate-like organic matter particles with aspect ratio *α _{i}* = 5 (width and length)/(thickness

*d*). This assumption is consistent with our SEM photomicrographs of flat organic matter platelets from the Clay City field site described in Wildman and others (2004). This was also the same assumption made by Lasaga and Ohmoto (2002) for their grains of OM, although, unlike their model, our code allows for dynamic changes in grain size during oxidation, implying that the surface area (a key factor in kinetic control) exposed to oxidation decreases while the OM volume decreases, even though the surface area to grain-volume ratio increases. The pyrite present in the field study was mainly in the form of framboidal spheroids. Although such “particles” have very high surface area, we argue that the outer edges of these will be the first to react with any available oxygen, thereby depriving the interior mass from oxidation. This would effectively make the reactive surface area confined to the outer perimeter of the framboids. With this in mind, our simple model used

_{i}*α*= 1 for pyrite (that is the particles of pyrite are modeled as cubes, which have similar size, area, and volume scaling as for spheres, within a small constant factor). Rectangular platelets with width and length

_{i}*L*, thickness

_{i}*d*, and

_{i}*α*=

_{i}*L*/

_{i}*d*, have surface areas of 2

_{i}*α*

_{i}d_{i}^{2}(2+

*α*) and volumes

_{i}*α*

_{i}^{2}

*d*

_{i}^{3}per particle, yielding an area to volume ratio of 2(2+

*α*)/(

_{i}*α*) per particle. In terms of the volume fraction of OM or pyrite, the number of particles of type

_{i}d_{i}*i*in a 1 m

^{3}reference volume, per volume, is

*N̄*=

_{i}*φ*/(

_{i,solid}*α*), with units of number per m

_{i}^{2}d_{i}^{3}^{3}.

*φ*is the volume of

_{i,solid}*i*type grains per bulk volume.

We used *υ _{i}* =1 for OM, and

*υ*=15/4 for pyrite (appropriate for final end-products of oxidation). We used 15.171 g OM per mole of C in OM, as calculated from the assumed OM composition of C

_{i}_{40}H

_{48}O

_{2}N

_{1}S

_{1}(Petsch and others, 2000; Petsch, ms, 2000). Holland’s (1978) summary of average shale compositions range from 0.69 to 0.94 mass per cent of organic carbon (see also Raiswell and Berner, 1986). We have used 1 percent as the choice for many of the introductory model runs.

Values for the diffusion coefficient of oxygen and air were adopted from Washburn (2003, v. 5, p. 62) of the form D_{0} = 0.178 **·**(cm^{2}/s) **·**(T(K)/273.15)^{1.75} at a pressure of 1 atm, then converted to cm^{2}/yr.

There is no uniformly consistent notation for the tortuosity factor that relates diffusion rates in air to diffusion rates in porous media, as evident in Cornell and Katz (1953), Troeh and others (1982), Sallam and others, (1984), Pape and others (1987), Boudreau (1997), Pape and Clauser (1998), Grathwohl (1998), Moldrup and others (2000, 2003a, 2003b), Scanlon and others (2002), Bartelt-Hunt and Smith (2002), and Boudreau and Meysman (2006). With our general form of Archie’s law, which using *φ _{a}D_{s}/D_{o}* =

*b*

_{T}*φ*corresponds to

_{a}^{n*+1}*D/D*of Troeh and others (1982); our

_{o}*T̄** is

*T*; and our

_{TROEH}^{2}*φ*is

_{a}*S*. If the matrix is more resistant to diffusion for low porosities, our use of

_{TROEH}*b*=1,

_{T}*n*=*2 corresponds to

*K*=1,

*m*=3 or

*u*=0,

*v*=3 of Troeh and others (1982), and is more resistant to diffusion than

*n*=*1 by a factor of

*φ*. We view this as a reasonable choice for fine-grained, compact shales.

_{a}There are partially compensating effects of faster reaction at higher temperatures, versus lower Henry’s law fractionation at higher temperatures. However, we performed numerous model runs that indicated that the increase of reaction rates with temperature by far overpowers the decrease in rates arising from lower dissolved oxygen that accompany higher temperatures.

The rate of OM oxidation was based on the experimental study (see also Chang, ms, 1997; Chang and Berner, 1998, 1999) of the kinetics of coal oxidation in aqueous solutions containing various measured quantities of dissolved oxygen, see figure 2A. Coal was used because of the inability to obtain sulfur-free marine kerogen. The relative oxidation rate of coal versus kerogen is not known, but the rates are probably similar. We are forced to use coal because we know of no experimental study of the kinetics of kerogen oxidation in water. In order to simulate postulated increase in kinetic rates of oxidation due to microbial activity, we ran some models with rates up to two orders of magnitude faster than the experimentally derived abiotic rates. Temperature dependence of the rate is accounted for using an apparent activation energy of 42 kJ/mole O_{2} (Chang and Berner, 1999). is the maximum rate of the OM oxidation reaction, which can be related to the *R*_{max}^{i} parameter described above by: where *T _{sim}* is the temperature of the simulation (in Kelvin) and

*T*is the temperature (K) of the experimentally derived kinetics, and

_{exp}*E*is the activation energy appropriate to capture the temperature dependence of the kinetic rate. From Petsch (2000) we used (mole C or O

_{a}_{2}consumed/m

^{2}/yr) and

*K*= 1.787×10

_{m}^{−6}(mole O

_{2}/cm

^{3}air). The conversion between

*K*(if estimated from micromole O

_{wat}_{2}/Liter of water) and

*K*is described in Appendix A, and Henry’s law partitioning is exploited to account for the difference between oxygen concentrations in the water and in the air. Figure 2A shows the OM oxidation rate for coal as a function of O

_{m}_{2}dissolved in water at 24°C from the data of Chang and Berner (1999). We fit the data with a Michaelis-Menten rate law fit as used by Petsch (ms, 2000) and a power law fit. The square-root dependence fit used by Lasaga and Ohmoto (2002) underestimates the rate by about a factor of two for low oxygen levels compared to the other two fits. A vertical bar is placed where oxygen in water is in Henry’s law equilibrium partitioning with current levels of atmospheric oxygen at 25°C (labeled PAL).

Pyrite oxidation kinetics have been studied by Smith and Shumate (1970), McKibben and Barnes (1986), Nicholson and others (1988), Williamson and Rimstidt (1994), Jerz and Rimstidt (2004), and Gleisner and others (2004). The influence of oxygen concentration on the rate of pyrite oxidation was examined by Smith and Shumate (1970), whose abiotic rate data is shown (closed circles) on a log scale in figure 2B, along with the power law fit we used for pyrite oxidation *A _{pow:Texp}* (rate of pyrite oxidation (in mole pyrite consumed/m

^{2}/s) =2.013×10

^{−12}× [O

_{2}in water in micromole/ L]

^{0.6696}, where we ignored data point #1 of their Table 3 that had an obvious irregularity in the ratio of column 1 (O

_{2}pressure in cm Hg) to column 3 (liquid O

_{2}concentration) compared to the other data points). Although Smith and Shumate (1970) do not report the pH for the runs in their Table 3, and Williamson and Rimstidt (1994, their Table 1) assume the Smith and Shumate (1970) Table 3 is for pH=2, our own calculations suggest the pH must have been in the range of 3 to 5, and we expect reasonably low pH for the sediment waters in contact with pyrite. Oddly enough, this nearly 2/3 power law dependence on oxygen concentration was not noted in any of the papers cited above. The influence of microbial activity was found to increase rates of oxidation by up to almost an order of magnitude. Both biotic (circled G’s) and an abiotic (circled A) oxidation data points from Gleisner and others (2004) are also shown in the figure, whose abiotic kinetic rates nearly match the abiotic rate found by Smith and Shumate. Smith and Shumate also measured some biotically controlled rates, with similar kinetic acceleration. Conversion of the power law form above to the factor

*A*is described in Appendix A.

_{i}Other parameters were varied over reasonable ranges. Initial porosity was varied from 2 percent to 10 percent, tortuosity from *n*=*0.5 to *n*=*2, OM content from 0.44 percent to 10 percent by mass of TOC, and grain thickness from 5 to 20 microns.

The mean global physical erosion rate today is 8.5 cm/kyr with mean values for continents ranging within a factor of ten (Berner and Berner, 1996). This value should also apply to shales, which undergo mainly physical erosion due to their content of relatively chemically-resistant components (except organic matter and pyrite). The present global erosion rate is about two times higher than the past, due to the effect of human activity, mainly agriculture, on erosion (Berner and Berner, 1996). About 4.2 cm/kyr is then appropriate for physical denudation of shales in the pre-agricultural earth. This rate reflects only eroded material that is delivered to the sea. Much eroded material is deposited in valleys, hollows, lowland soils, *et cetera*, and is subjected to additional oxidative weathering of any transported organic carbon during storage (see Blair and others, 2004). For carbon cycle modeling, we can ignore this kind of erosion, as only the eroded organic matter carried to the sea and buried without oxidation is relevant to our discussion; such terrestrial carbon reservoirs are transient over geologic time, and the mass of carbon stored in these reservoirs is dwarfed by the mass of crustal carbon. The 4.2 cm/kyr value above is close to the 5 cm/kyr used by Lasaga and Ohmoto (2002) for mean continental erosion rates, who estimated this from Holland (1978). Also, because of relatively recent continental glaciation, present global erosion rates are higher than they have been since about 400 million years ago (Berner, 2004, Chapter 2, p. 19). This means for geologic time, 4 cm/kyr is a likely maximum global average erosion rate. We varied erosion rates from 1 to 50 cm/kyr in our modeling.

Atmospheric oxygen levels were typically set at the present level (20.9476% by volume, or by moles), but were varied from 1 percent to 40 percent. The model presented is relevant to the Phanerozoic time, during which 12 percent oxygen is the likely lower limit (Berner, 2004, 2006). The values below 10 percent are considered only for comparison to the Lasaga and Ohmoto (2002) model where their figures 4 and 5 display results for values as low as 10^{−8} PAL (with 1 PAL corresponding to the present atmospheric level of near 21%) and their figures 6 and 7 show results for levels as low as 0.001 PAL. Note that 10 percent oxygen corresponds to 0.5 PAL.

## RESULTS AND DISCUSSION

### Models of Organic Matter Oxidation

We first present model results on the oxidation of OM without pyrite present. The oxygen reducing capacity in typical shales is dominated by OM compared to pyrite. In addition, typical shales have concentrations of OM on the order of 0.7 to 1 percent (Holland, 1978), much lower than the black shale of our field study (Wildman and others, 2004) that will be discussed later. The first case discussed will be used as a default model run to which most of the sensitivity studies will be compared. Parameters for the default model run are summarized in table 1. Most of these parameters are chosen in a conservative way, in the sense that they favor finite amounts of OM left at the eroding surface. The one exception to that is the parameter *s _{air}* (as described below).

### Basic dynamics and front characteristics.

In figure 3 we have plotted the evolution of the **TOC** percent (total organic carbon mass percent, that is 100 × (mass of C in OM)/( bulk dry mass)) and concentration of O_{2} in the air of the pores as a function of depth. The oxygen level of the atmosphere at the surface was fixed at current-day levels (20.9476 volume %), but was initialized at 0 percent for all depths below the surface node. The erosion rate was fixed at 5 cm/kyr. This simulations used *n** = 2, *d*=10 μm, and *s _{air}* = 0.9 (an air saturation of 90%). The Michaelis-Menten kinetics were used for OM oxidation (compare to fig. 2A), but the rates were adjusted to 17°C as described above. The time step was 0.8 yr and there were 220 grid points. The initial condition was 1 percent TOC at all depths for OM, as shown by the vertical line at 1 percent TOC. The lower boundary condition, well below the reaction zone, was maintained at 1 percent TOC and 0 percent oxygen. The porosity (pore/bulk volume) fraction was initiated at 0.05, and was allowed to vary as the OM was oxidized (taking into account the volume removed by oxidation that became pore space). As most soils generally have considerably larger porosity, these calculations are conservative (they favor OM close to the surface), in that low porosity is associated with inhibited diffusion of O

_{2}downward. The initial conditions are representative for typical shales freshly exposed to the atmosphere.

The succeeding profiles show a deepening “front” where OM makes a transition, from total absence near the surface, up to its value at depth, over a characteristic depth range. The time interval between profiles is 96,000 years (for both fig. 3A and 3B). For the intervals shown, the oxygen profiles exhibit a linear trend (constant gradient) near the surface, indicating diffusive flux down toward the reaction zone where oxygen is consumed and the profile develops some curvature. The system evolves toward a steady state where the OM supply from below is balanced by oxygen supply from above. The reaction zone has a finite thickness that depends on the rate of the oxidation reaction (dependent on the reaction kinetics and evolving OM “grain” size).

Steady-state profiles for several quantities in the default model run are shown in figure 4. In figure 4A, the final profiles for OM and oxygen are shown. The form of the OM front is not symmetric: the transition from OM absence toward increasing OM exhibits higher curvature than the deeper transition toward unoxidized OM. As noted above, the oxygen profile is linear above the OM front, but exhibits curvature in the reacting zone. Figure 4B shows the oxygen consumption reaction rate profile for the steady state, clearly showing that at steady-state the reactions proceed only in the vicinity of the front. Figure 4C shows the steady profile for air-filled porosity. The zones of non-constant porosity give rise to the additional spatially varying diffusive terms for *φ _{a}* and the tortuosity

*T̄**. The magnitude of the diffusion flux terms (compare to eq 5) have the final profiles as shown in figure 4D, where: and Δt for this simulation was 0.8 yr. Other simulations, especially those with larger OM content, showed that the magnitudes of dφ and dD could exceed df, indicating that it is important to retain each of these terms.

Additional aspects of the time dependence of this simulation are shown in figure 5, where in part 5A we show the depths of various OM concentrations. The upper curve shows the depth evolution where the OM concentration is 95 percent of its initial value (or 95% of the value at depth). (The injection rate from below of bulk volume/bulk area is equivalent to the erosion rate at the surface). The middle curve corresponds to 50 percent of the initial TOC value, while the lower curve corresponds to 5 percent of the initial value. We use the middle curve (50% value) to define the “front depth”. We define the “front width” as the depth difference between the upper and lower curves (the 95% and 5% values compared to the initial value or value at depth). The evolution of this characteristic front width is displayed in figure 5B (its initial value is undefined as the system starts with 100% of the initial value everywhere). Both the front depth and front width show a trend toward a steady-state value. We define a characteristic “time to steady state” as the time required for the front depth to reach 95 percent of its final value, which for this case was 437 kyr. This characteristic time is not exactly unique, in that estimates based on other properties would be slightly different. This is evident in figures 5C and 5D. Figure 5C shows the evolution of oxygen at several fixed depths. At any locale with finite oxygen values, the oxygen concentration is increasing toward the steady state value. Figure 5D shows the evolution of TOC at several depths. It is important to recall that these depths are in the frame of reference fixed relative to the eroding surface. As material is supplied from below, and eroded off the top, some depths experience rapid OM removal by oxidation (for example at the surface) while other depths retain their supplied value (those well below the front).

#### The influence of reaction rate on front width.—

We can study the influence of the rate parameter on characteristics of the front by artificially manipulating the rate constant *R*_{max} of equation 3 (as described by our discussion of fig. 2B, some studies have found that reaction rates are accelerated by biotic effects). The front depth, as defined above, changed less than 5 percent for an increase of the kinetic rate constant by a factor of 100. On the other hand, the changes of the front thickness were significant, as displayed in figure 6. In figure 6A we show the dependence of the front width on the ratio of *R*_{max} to its default value from figure 2A. The front width (around 300 cm for the default model run) decreases to about 30 cm when the rate is increased by 100 times. Results are shown for the default *R*_{max}, as well as 2, 10, and 100 times the default value (all other parameters are as those used for the default model run illustrated in figs. 3–5). Results are also shown for various grid resolutions, where the larger open circles are for the lowest resolution (220 grid points in the domain), the next smaller circles are for 420 grid points, and the smallest open circles are for a resolution of 820 grid points. Figure 6B shows the same results on a log-log scale with a line for reference that has a −0.5 slope. The front width approximately decreases by a factor of 10 for a 100 fold increase in the kinetic rate parameter (an inverse square-root dependence). As the front sharpens, higher resolution is needed, to the point that the 100x results are not quite resolved, even with 820 grid points. The 220 grid results of the default model run (as used for figs. 3–5) at the default kinetic rate are accurate to within about 2 percent for the front width, and are significantly more accurate for the front depth (error less than 0.5%).

#### The influence of air fraction in pores on front depth.—

The air saturation within the pores, *s _{air}*, has a significant influence on the front depth. This is one parameter that we have chosen in a way that does not minimize the front depth for the simulations above. As we restrict our model to the vadose zone (above the water table) and as we presume the air fraction in the pores would be well above 0.5 for much of the year, our simulations up to now have used

*s*=0.9. More sophisticated models would follow the annual cycle with precipitation and evaporation, wetting and drying, as well as water infiltration and air being sucked into the soil layer as a water-saturated pulse drained down toward the water table. The latter would draw fresh oxygenated air down into the soil layers by advection, enhancing the exchange over what our diffusion model can follow. This advective exchange would counterbalance the diffusion inhibition that low air saturations would have during the wet part of the year. For complete water saturation (zero air saturation) oxygen diffusion rates into the soil would be small, allowing little OM oxidation. We illustrate the influence of

_{air}*s*on the front depth in figure 7. For the lower values of

_{air}*s*(0.4 and 0.35), not all the OM is oxidized before reaching the surface (as indicated by the large circles). One caveat regarding the influence of the air/water ratio: the oxidation of OM in arid conditions is expected to be low; a water film seems to be necessary, otherwise the oxidation occurs much more slowly than described by Chang and Berner’s (1999) kinetics. OM oxidation in arid regions may be thus inhibited (compare Leythaeuser, 1973; Clayton and Swetland, 1978).

_{air}A series of simulations (not shown) were done in the presence of larger initial concentrations of both OM and pyrite (TOC of 8.3 mass %, pyrite of 13.1 mass %), where all parameters were held fixed except for *s _{air}*. Compared to the results displayed in figure 7, the high OM and pyrite runs exhibit a reversal of the curvature, that is for air saturation values above 0.5 the front depths remain relatively constant, but drop precipitously for

*s*below 0.45. The front depth varied about 10 percent upon a decrease of

_{air}*s*from 0.9 to 0.5. Further decrease of

_{air}*s*from 0.5 to 0.35 caused the front depth to decrease from around 3.4 m to 1.2 m.

_{air}There are many parameters that must be set for a model run. Changing one parameter can have a similar effect as changing another in certain cases. For instance, setting the porosity to 0.1 and the air saturation to 0.5 will produce the same results as the case with a porosity of 0.05 and an air saturation of 1.0, as long as the volume of material oxidized is negligible. On the other hand, even if large quantities of OM are oxidized, the imposition of an upper limit for the air porosity (as an approximation for the effects of compaction) can be equivalent to lower air saturations, as *φ _{a}* =

*s*

_{air}*φ*, where

*φ*is the total porosity. It is

*φ*that is used in the calculation of the tortuosity (and the effective diffusivity). Our choice of initial total porosity is lower than some shales would have. Model results with higher porosity, but lower air saturation (higher water saturation) could yield the same air porosity near the surface. Some of these issues become more important when larger amounts of OM can be oxidized, and are discussed more in the next subsection.

_{a}#### The influence of the injected magnitude of TOC on front depth and width.—

The OM content of the sediment injected from below at depth (the lower OM boundary condition) has a complex influence on the front depth. If little OM exists at depth, the front tends to be deep. High OM content is associated with shallower front depths. One additional aspect of the model becomes apparent for the cases with large quantities of OM at depth. Although the porosity generally increases as OM is oxidized, we imposed a “lid” value of 0.1 for the air-porosity fraction for many of the simulations, that is, we did not allow the air-porosity to exceed this lid value. The use of a lid porosity is a simple approach to capture the effect of finite, but limited, compaction of the medium. The trend of front depth as a function of TOC (or OM) is not necessarily monotonic. Upon varying the OM concentration at depth, but otherwise using the default parameters as outlined in the first simulation (but using an 0.1 lid air-porosity fraction), a local maximum of front depth occurs around 4.5 percent TOC, as displayed in figure 8A. A local minimum front depth occurs around 2 percent TOC. The broad trend of shallower front depths as TOC increases is consistent with increased oxygen gradients required to oxidize larger quantities of OM. However, from 2 percent TOC to 4.5 percent TOC, the front becomes deeper, which we interpret to be due to the important effect of increased pore spaces opened up for oxygen diffusion as more OM is oxidized near the surface. For TOC above the local maximum of the front depth around 4.5 percent TOC, the imposed lid air-porosity fraction of 0.1 is reached (as discussed above, as a simple approach to deal with compaction), limiting further increase in the air porosity (which also limits the effective diffusivity). We contrast the above behavior to simulations run with a constant porosity of 0.05 (the same as the initial porosity in all the simulations) in figure 8A, but we do not view such constant porosity simulations as realistic, as they disallow porosity increases that would accompany OM oxidation. The use of a strict lid air porosity is not the ideal way to deal with compaction. Future modeling should use a more realistic compaction formulation, with compaction entering more gradually as the porosity increases, yet asymptoting to some value where the solid matrix would not consolidate further. We expect that improved compaction modeling would not remove the relative maxima of front depth in figure 8A for the variable porosity case, but that the cusp-like nature of the curve would be smoothed out. As discussed in the subsection above, the use of a lid value for the air porosity can be equivalent to a system with lower air saturation (greater water saturations) than we set as a parameter.

Figure 8B shows how the front depth compares with the front width for the same simulations. Generally the front width decreases somewhat as the front depth decreases. The “width = depth” line is also shown. When the width is comparable to the depth, finite amounts of OM typically remain unoxidized at the surface. More than 0.1 percent TOC remained at the surface for the constant porosity runs of 6 percent or greater initial, and deep injected value, TOC (the five “+” points closest to the “width = depth” line).

#### The influence of initial OM grain size on front depth and width.—

The initial OM grain thickness, *d*, was varied from 2 to 50 μm, while keeping all other parameters as the default model run (including the initial TOC of 1%). These results are shown in figure 9. The front depth is rather insensitive to the initial grain thickness, but the front width increases rather dramatically with initial grain thickness. This is primarily a surface area effect, as small grain sizes have large areas presented to the oxidative influence. We assumed that the grains retained their shape during oxidation (*α* =5), but of course the thickness *d* decreases via oxidation during the simulations.

#### The influence of tortuosity on front depth.—

The tortuosity exponent *n** is the focus of the next simulations. The above cases all used the rather conservative value *n**=2. This corresponds to a matrix that is relatively tight for diffusion, in that for a given relatively low porosity, the path for diffusion of oxygen is rather tortuous. For low porosities, increasing *n** implies decreasing *D _{s}*. The effective diffusion coefficient

*D*becomes independent of porosity for

_{s}*n**=0. As diffusion is enhanced for lower values of

*n**, the front depth increases as

*n**decreases, as shown in figure 10. The default value of

*n**=2 is the extreme right point on the plot (with a front depth of 977 cm –note front depth scale in meters). We admit that some of these simulations predict front depths that are quite unrealistic. Our model is only relevant above the water table where oxygenated air can diffuse downward. Water tables as deep as the low

*n**results are probably rare in most shales, other than those in arid environments. Even so, the plot well illustrates the influence of changing the tortuosity formulation. The front thickness also increases with increasing front depth (for this case as the square root of the front depth). We have had difficulty finding published tortuosity estimates for tight shales, and thus have used a conservative value for our default value of

*n**. We state this is a conservative choice in the sense that it favors front depths that are shallow, thus tending to favor oxygen-level feedback on the OM oxidation. As long as all OM is oxidized before reaching the surface, no such feedback exists. The default value of

*n**= 2 seems to us to be appropriate for shales that are quite resistant to easy diffusion of air. Future models should consider the effects of heterogeneities, cracks, and fractures in two dimensions.

#### Comparison to parameters in Lasaga and Ohmoto (2002).—

In figure 11 we show the profiles of oxygen and OM for three tortuosity formulations, using their case of 0.44 percent TOC (Lasaga and Ohmoto, 2002; for a global average of shales, carbonates, and sandstones). Besides the choice of 0.44 percent TOC, these models used the default parameters used for the default model run above (including our assumed temperature of 17°C for which the reaction kinetics would be slightly slower than their 25°C assumption) except as noted. Figure 11A used *n**=2, figure 11B used *T̄**=0.2*φ _{a}* (that is

*b*=0.2,

_{T}*n**=1), and figure 11C used

*n**=1. A more open matrix (higher

*T̄**) leads to deeper fronts as expected. We do not claim that such deep fronts (as part C) are physically relevant, as most environments would encounter a water table well above the front depth shown in figure 11C. None of the model runs displayed in figure 11 retain finite amounts of OM at the surface. This is in contrast to figure 4 of Lasaga and Ohmoto (2002), where finite amounts of OM

*are*retained at the surface for the erosion rate of 5 cm/kyr for both

*L*= 1 m and

_{s}*L*= 10 m. We note that we have not used their factor

_{s}*α*

_{conv}, but that cannot account for this difference.

#### The influence of oxygen level and erosion rate on front depth.—

Figure 12 shows the depth profiles for OM and O_{2} for the case with *n**=1, 15 percent atmospheric oxygen, and an erosion rate of 50 cm/kyr. The front widths are generally broader with *n**=1 than with *n**=2. Both the front depth and width can be influenced by erosion rate and oxygen concentration of the atmosphere. We now examine this systematically for different *n** and initial TOC.

We display in figure 13 the variation of OM front depth for two different deep OM concentrations (lower boundary condition of 1% TOC in part A and C, and 5% in part B) for two different tortuosity estimates (*n**=2 in parts A and B, and *n**=1 in part C) as a function of oxygen level in the atmosphere and erosion rate. The front depths increase essentially linearly with oxygen level of the atmosphere. Doubling oxygen, doubles the front depth, as long as the front is buried. As erosion rates increase by a factor of 10, the front depth is essentially 10 times shallower. The front depths for the two deep TOC levels in parts A and B are very closely the same. This is due to the fact that with an air porosity lid imposed of 0.1 (compare fig. 8A and text discussion above) front depths for 1 percent and 5 percent deep TOC were nearly the same value. Also, by comparing to figure 8, we could presume that 0.5 percent initial TOC values would yield considerably deeper front depths, while 10 percent initial TOC would yield shallower fronts. In figure 13C, we note the significant influence of the tortuosity parameter, where the *n**=1 cases yield deeper fronts than *n**=2 (other parameters being equal). Recall that *n**=1 yields the typical Archie’s law, and that *n**=2 is our conservative choice (a tighter matrix inhibits diffusion). For *n**=1 and rapid erosion (50 cm/kyr) finite amounts of OM reach the surface, as noted by TOC percent on that erosion rate curve of figure 13C.

#### The time elapsed for steady state to be achieved.—

In figure 14 we show the characteristic “time to steady state” as the time required for the front depth to reach 95 percent of its final value. Most of our simulations started with *φ* =0.05 and used *n**=2. For these simulations, the time to steady state data shown in figure 14 could be collapsed to a nearly linear trend as a function of front depth, if the ordinate were replotted as the time to steady state multiplied by the erosion rate ω. Plotted in such a fashion, these simulations had a slope of 2 (nondimensional). Additional simulations, particularly with *n**=1 and various values for the initial porosity, did not follow the linear trend just described and highlight the lack of a simple scaling for the “time to steady state”. The two points with time to steady state in excess of 2500 kyr in figure 14 correspond to simulations with low atmospheric oxygen and low erosion rates (5% atmospheric oxygen and 1 cm/kyr erosion rate, the upper point had 5% initial TOC, and the next highest point had 1% initial TOC). These results for time-to-steady-state could be of value when evaluating field sites that may have been disturbed by glacial events, land slides, *et cetera.* “Times to steady state” as long as several hundred kyr for many of model runs suggests that in many landscapes, more recent surface disturbances such as landslides and slumps may reflect weathering profiles that have not yet achieved steady state. Depending on the recurrence interval between such disturbances, the number of landscapes that genuinely achieve steady state may be limited.

#### The balance of fluxes of oxygen and OM.—

In figure 15 we consider a simple scaling relation for the case of “buried fronts”, that is cases where all OM is oxidized before reaching the surface. At steady state, the governing equations simplify considerably. One can derive an expression for the downward flux of oxygen that is balanced by the upward flux of OM (dependent on the concentration of TOC injected at depth and the erosion rate). Such a balance has been verified numerically, and is a good test of the resolution required for accurate solutions. This approach yields where |* _{deep}* implies that

*g*is the value at depth. A more approximate estimate may be useful, where the

_{OM}*da/dx*term is estimated by the atmospheric oxygen concentration (minus zero at depth) divided by the front depth (with appropriate unit conversions). Such an estimate is supported by the linearity of the oxygen gradient in figure 4A, where the oxygen versus depth slope is nearly constant until the front is reached (at depth shown in fig. 5A). If the front depth by itself were a perfect estimate for the

*dx*part of the gradient in this flux, replotting the data using this estimate would collapse all the data to a single line. Of course fronts have various widths, dependent on particle size and tortuosity in the front vicinity, so this estimate is not exact. In figure 15 we show the results of such a scaling estimate from our simulations. The lower line corresponds to the case where the front depth is the perfect depth scaling for the oxygen gradient. Its slope is 470.8 when the ordinate has units of (O

_{2}volume % )* (cm/yr) and the abscissa has units of (mass fraction TOC) * (cm/kyr). The cases with shallow fronts generally lie on the right in this figure, while deep front cases are on the left. The upper line bounds most of our simulations and has a slope of 750 in the same units. If the TOC at depth, the front depth, and an average tortuosity were known above the oxidizing OM, the simple scaling could be used to estimate the local erosion rate! Alternately, results of field studies may be used to estimate the effective surface layer diffusivity if the erosion rate were known independently.

#### The influence of oxygen level and erosion rate on TOC remaining at the surface for n*=1.—

In figure 16 we have plotted the TOC remaining at the surface as a function of the oxygen level of the atmosphere for five different rapid erosion rates (10, 20, 30, 40 and 50 cm/kyr), starting with 1 percent TOC (by dry mass) at depth. These simulations were initiated with 5% porosity, *n** =1, *d*=10 microns, and *s _{air}* =0.9. Only the highest erosion rate yields significant OM left at the surface (that could be transported elsewhere and potentially reburied before being oxidized, for example Blair and others, 2003, 2004) for oxygen levels above 10 percent. For this set of parameters, the model does not predict any residual TOC at the surface, at any of the indicated erosion rates, for atmospheric oxygen levels above 15 percent. The figure clearly indicates that mean continental erosion rates (∼ 4 cm/kyr) are not high enough to produce any oxygen dependent feedback for oxygen levels above 2 percent (by volume), at least for the Chang and Berner (1999) kinetics and the other model assumptions (such as

*s*). The value of

_{air}*n**=1 is typical for Archie’s law in most soils. The fact that the porosity never exceeded 6 percent in the simulations shown in this figure (the initial 5% plus some space opened up by the oxidized OM) is conservative, in the sense that larger surficial porosity would favor deeper fronts, and less OM surviving at the surface for steady state.

#### The influence of oxygen level and erosion rate on TOC remaining at the surface for n*=2.—

In figure 17 we have plotted the TOC remaining at the surface as a function of the oxygen level of the atmosphere for four different erosion rates (5, 10, 20, and 50 cm/kyr). Results for erosion rates of 1 cm/kyr results (not shown) had no TOC remaining at the surface for oxygen levels of 5 volume percent or higher. These simulations used *n** =2, *s _{air}* =0.9, as well as lower boundary conditions of

*d*=10 microns, with 1 percent TOC (by dry mass) for part A, and 5 percent TOC for part B. The fractions of TOC remaining at the surface compared to the deep injected value are quite similar for parts A and B of this figure. This is somewhat dependent on our use of a lid air porosity as a simple means of dealing with compaction. Recall that from figure 8A, the front depths for these two deep TOC values were similar. In general (also from fig. 8A) one would expect lower values of deep injected TOC (<1%) to yield deeper fronts and less TOC remaining at the surface. Results for TOC of 0.44 percent (deep injected value) and a 20 cm/kyr erosion rate are also shown in figure 17A. Only the highest erosion rate yields significant OM left at the surface (that could be transported elsewhere and potentially reburied before being oxidized). Again, within the context of this model, this clearly indicates that mean continental erosion rates ( ∼ 4 or 5 cm/kyr) are not high enough to produce any oxygen dependent feedback for oxygen levels above 12 percent and are thus insufficient to leave finite amounts of OM at the surface for oxygen levels that have been estimated for the Phanerozoic (Berner, 2004). Recall that our choice of the tortuosity parameter

*n**=2 was very conservative and that more typical values (e.g,

*n**=1) yield much deeper fronts (figs. 10 and 13), which would further discourage finite amounts of OM remaining at the surface. In addition, for shales with smaller TOC than 1 percent, front depths would also tend to be deeper (compare to fig. 8A). We conclude that the flux of oxygen consumed by weathering of OM in shales is generally set by the OM content at depth and the erosion rate (as concluded by Holland, 1978). This also implies that an oxygen level dependent feedback for OM weathering is generally not appropriate for mean continental erosion rates and Phanerozoic oxygen levels.

### Models of Organic Matter and Pyrite Oxidation and Field Site Comparison

We now present results of modeling pyrite oxidation in addition to OM. Figure 18 displays results of the field data from our previous study (part I of this study: Wildman and others, 2004) as well as results of several model runs that include both OM and pyrite. The field data came from a 13 m deep vertical core drilled into flat-lying layers of the New Albany Shale (late Devonian, 365 Ma) from a hillside near Clay City, Kentucky. In figure 18A and B, the upper 10 m of the depth profiles are shown. Part B shows the raw data, and part A shows the averages of the data at each depth where analyses were made above 10 m (below 10 m there is an apparent change in the original OM and pyrite content of the layers due to non-steady state deposition—see Wildman and others, 2004). The data include mass percentages of TOC, total S, and the porosity volume percent (the porosity could not be measured in the upper 2 m of the profile, due to the friable nature of that section of the core). There is a clear trend with increasing depth: of decreasing porosity and of increasing OM and S. The S in the lower section is primarily in pyrite, which increases quite smoothly over the depth range below 2.5 m. The OM front is sharper, as observed in the TOC increase below 2.5 m. Aside from oxidation of OM near the surface, some of the variation in TOC at depth is possibly due to variations in the initial TOC, “fast pathways” for oxygen diffusion, or other causes. Significantly, the pyrite is essentially absent when the OM has increased significantly just below the “front”. This is consistent with kinetics of pyrite oxidation that are more rapid than ancient OM oxidation (compare to fig. 2). Except for the data nearest the surface, where the OM is certainly “modern”, the other near-surface data between 1 and 2 m have small amounts for OM (≤ 0.3% TOC). The data indicates that the deep portion of the core had a porosity of near 5 percent, which we have used in most of our modeling as initial and lower boundary conditions. If the volume of OM and pyrite (proportional to measured TOC and S) at around 9 m represent an initial value, the total removal of this volume by oxidation of OM and pyrite near the surface would open up the porosity to just over 20 percent. Our modeling (fig. 18C) shows general agreement between the model and the porosity data below 2 m. Very near the surface, additional processes can contribute to increased porosity (due to flora and fauna), but our focus is primarily below this zone. Other minerals could weather as well, contributing to porosity increase, but such weathering is likely a minor component of the porosity increase compared to the weathering of pyrite and OM. Weathering of pyrite does produce sulfuric acid that could dissolve silicate minerals, also leading to greater porosity. Compaction could also play some role in the higher porosity section above 3 m, but this is difficult to evaluate from our data.

Figures 18C–18F show model profile results at steady state. These results used MM (Michaelis-Menten) kinetics for OM and power-law kinetics for pyrite oxidation. In figure 18C, the TOC values for OM superimpose the results of two models, with the bold curve arising from MM kinetics, and the thin curve showing nearly identical results from a power-law fit to the data of COAL25-125 of Chang and Berner (1999). The oxygen concentration has a linear gradient near the surface. Diffusive flux down this gradient feeds oxygen to the zone of reaction in the OM and pyrite fronts. Some curvature of the oxygen profile in this zone arises from the varying consumption, as well as the variations in porosity and tortuosity. The oxygen concentration vanishes where the OM and S (in pyrite) recover their deep values (supplied from below). The increase in TOC with depth is more gradual than pyrite. The pyrite front is very sharp, and generally deeper than the OM front, both consistent with the imposed faster kinetics of pyrite oxidation than for OM. This model used *s _{air}*=0.45, present levels of atmospheric oxygen, a 5 cm/kyr erosion rate, initial OM thickness of 10 microns, initial pyrite thickness of 20 microns, 8.3 weight percent TOC, 13.1 weight percent pyrite (about 7.00 wt.% pyrite S) and

*n**=2. The porosity increase from depth toward the surface in the model is simply a result of volume loss of OM and pyrite.

There are some important differences between the observed data and the model results. For the field data, the TOC and the S appear to vanish together at the upper part of the front. The model results for these quantities are different in detail: The transition depth where S increases to finite values is significantly deeper than the corresponding transition depth for TOC. In addition, modeled S and TOC achieve their maximum near the same depth at the bottom of the front. In addition, the pyrite front is sharp in the model results, unlike the data. This may be due to protective coatings on the pyrite in the field site or other factors discussed below. We have not attempted to solve an inverse problem to match the data as closely as possible. The kinetics, tortuosity formulations, erosion rates, and compaction models could be varied in an attempt to fit the data better, but such an exercise is beyond the scope and purpose of this study.

We also present, for comparison purposes, several other simulations with OM and pyrite present. In figure 18D we show results for the same parameters used for figure part C, except we used *s _{air}*=0.9, initial OM thickness of 5 microns, and an imposed “lid”

*φ*of 0.1. The resulting front depths are somewhat deeper, and the porosity “front” is much sharper. The curvature of the oxygen profile in the reaction zone is different from the last case (arising from differences in the air porosity and resulting tortuosity).

_{air}Two additional simulations (figs. 18E and 18F) with 10 percent atmospheric oxygen (about half present values) and 10 cm/kyr erosion rates (about double pre-agricultural values) are presented, otherwise using model parameters as for figure 18C, unless noted. For figure 18E, *s _{air}*=0.9, and in figure 18F,

*s*=0.45 (note different vertical scale). The higher air saturation run of part E exhibits a broader OM front, while the lower air saturation run retains finite quantities of OM at the surface for these conditions of rapid erosion rate and low atmospheric oxygen. Note that the large initial deep values of TOC and pyrite tend to favor shallower fronts than would cases with more typical shale values, as shown in previous plots.

_{air}We have based our OM oxidation modeling on the kinetics measured by Chang and Berner (1999), as was true for the modeling study of Lasaga and Ohmoto (2002). We must admit that there may be some portion of ancient OM that weathers more slowly than was measured by that study. This potentially resistant portion may also be influenced by just how the OM was deposited and preserved. Specifically, there is a strong correlation observed in modern marine sediments between organic carbon content and mineral surface area (Mayer, 1994, 1999; Hedges and Keil, 1995; Mayer and others, 2002; Dickens and others, 2006), with an approximate slope of 0.5 to 1 mg OC m^{−2} mineral surface area. Association with mineral surfaces (adsorption into micropores, interlayers within clay minerals, steric hindrance, *et cetera*), should exhibit slower oxidation kinetics than unassociated material. This relationship appears to extend from unconsolidated sediments to lithified, ancient sedimentary rocks (Kennedy and others, 2002, 2006). Additionally, OM in sedimentary rocks is not uniform in source or composition. Shales reflect OM derived from a variety of sources: intact algal remains, partially biodegraded material preserved during diagenesis, remains of sediment bacteria, and detritus from vascular plants. The myriad different OM sources, extents of degradation prior to lithification, and the compositions these materials exhibit, suggests that assigning a single uniform degradation rate constant is a simplification. In addition, particle aggregation may influence accessibility of oxygen.

Samples from the Clay City weathering profile exhibit surface areas of about 30 m^{2}/g (table 2, Wildman and others 2004). Typical OC-mineral surface area relationships would suggest a mineral-associated OC content of 1.5 to 3 percent, which is well below our observations of unweathered shale at this site. Typical of many OM-rich rocks, the majority of OM in unweathered Clay City shale appears not associated with mineral surfaces. TOC contents as low as 0.2 percent in weathered sections of the profile indicate that any mineral association does not significantly inhibit OM loss during weathering, Also, electron photomicrographs (Wildman and others, 2004) indicate abundant micron-sized “blebs” of organic matter not intimately associated with minerals. A quantitative microprobe scan produced an area percent of discernible black material (organic matter) in rough agreement with chemical analysis for organic carbon.

Furthermore, while mineral associations may cause OM to be more resistant to oxidation or degradation, these associations cannot render OM inert. If globally a large proportion of sedimentary organic matter were entirely resistant to oxidation during weathering, then appreciable OM would be continually eroded, re-deposited and accumulate over geologic time. This did not happen because it would have led to continual change in the δ^{13}C of the ocean and the O_{2} content of the atmosphere, which have not been found (Broecker, 1970; Berner, 2004). A simple calculation illustrates why re-deposition over long periods, of an appreciable proportion of total organic matter as totally O_{2}-resistant OM, is impossible. Over the past 500 million years the global rate of burial of organic carbon has ranged roughly between about 1 and 6 x 10^{18} mole C /Myr (Berner, 2006). Taking an average burial rate of 3 x 10^{18} mole C /Myr, assume that only 10 percent, or 0.3 x 10^{18} mole C /Myr, of this at each time of burial was old OM resistant to oxidation under all conditions. As this situation continued over many millions of years this would represent an increase in the mass of organic carbon and an excess of atmospheric O_{2} that was not removed by oxidation during weathering and that accumulated in the atmosphere. Over 500 million years at a rate of 0.3x 10^{18} mole C /Myr, this amounts to an excess production of 150 x 10^{18} moles of O_{2}. Yet there are only about 38 x 10^{18} moles of O_{2} in the present atmosphere. This argument cannot be refuted by the assumption that the proportion of OM burial that is unreactive increased with time to ultimately account for all OM burial. This could happen only if photosynthesis ceased and deposition included no new organic matter. The unlikelihood of the accumulation over the past 500 million years of a large proportion of totally O_{2}-resistant OM is also supported by the fact that the predicted carbon isotopic record for this situation is far different than that actually measured (Broecker, 1970).

An additional factor is that anaerobic decomposition of sedimentary OM occurs via biological methane formation (Martini and others, 2003) at depths below those that we have considered. This methane upon ultimate evasion to the atmosphere is quickly oxidized by atmospheric O_{2}. The overall sum reaction is the uptake of O_{2} by reaction with OM. Thus, there is an even greater chance that organic matter will not survive uplift and erosion than what is considered by our model, which involves only direct oxidation by O_{2}.

Resistance to oxidation in the weathering profile is not the same thing as complete resistance to oxidation. There may be some fraction of sedimentary OM that resists degradation during weathering, becomes incorporated into riverine OM pools, but is largely oxidized during transport in the river and before reburial. Oxidation of a large proportion of terrestrial organic matter during transport in rivers and after deposition in nearshore marine sediments is demonstrated by the data of Berner (1982). Ultimately, the organic matter becomes oxidized, but the locus of oxidation shifts from the weathering profile to the rivers and ocean for some fraction of sedimentary OM. This fraction is not considered in our weathering model. Some landscapes do transport ancient, particle-associated OM to the oceans for reburial (compare Masiello and Druffel, 2001; Raymond and Bauer, 2001a, 2001b; Blair and others, 2003, 2004; Leithold and others, 2005; Komada, Druffel and Hwang, 2005; Goñi and others, 2005). The results of this model suggest that these landscapes may be characterized by high erosion rates, recent disturbance that removed weathered rocks and regolith, or greater water saturation. However, unoxidized, reburied ancient OM does not accumulate on a global scale in a way that increases the total mass of the crustal carbon inventory. Otherwise oxygen levels and δ^{13}C of the oceans would have increased continually over geologic time, as discussed above.

We originally constructed our weathering model to attempt to explain the organic matter and pyrite distributions of a specific Devonian shale (Wildman and others, 2004). We acknowledge that the shale in the field location is not strictly homogeneous in either physical or chemical properties. We did find, however, that the degree of inhomogeneity of the section from which the core was taken was not significant for the carbon and pyrite sulfur at depth. In addition, the permeabilities were measured (generally on the order of 10^{−17} m^{2}) in both the vertical and horizontal directions, with the most extensive measurements made on a 12.28 m depth sample (compare to discussion in Wildman and others, 2004). The permeability decreased with increasing effective pressure (confining minus fluid pressures), and the permeabilities reported were those extrapolated to the lower effective stresses representative for field conditions. The difference between the horizontal and vertical permeability was within a factor of two at all effective pressures.

The surface slope was about 20° in the region where the vertical core was taken, while the shale layering was essentially horizontal. Our modeling implicitly assumes vertical diffusion downward into a shale that may develop an oxidized zone near the surface with horizontal “fronts”. One could expect that there may be some inhomogeneity in the case of gaseous diffusion due to the layering of the shale layers, and that diffusion along the layers may be easier than across the layers, in which case the relevant distance for diffusive transport would be the horizontal distance from the sloping surface to the location of interest, which is considerably longer than the vertical distance to the surface. Alternately, if the pore space interconnectivity was totally homogeneous, the relevant distance would be the shortest distance to the sloping surface, which for a 20° slope would be less than the vertical distance by a factor of 0.94. This factor translates into a ∼6 percent difference that may be relevant when comparing model results with the observations: we view the “shortest distance to the surface” as more relevant than the purely horizontal distance within the shale layers due to the similarity of horizontal and vertical permeabilities, as well as the rather similar appearance of the gradients in OM and pyrite sulfur between a vertical core and samples taken horizontally along a single layer (compare to Wildman and others, 2004). Also, vertical and horizontal thin sections, taken from the zone of the “OM front” through which oxygen must be passing, show similar textures. All this is to admit that we could not expect model/field comparisons to be precise, given the inherent geometry and heterogeneity of the field site.

It would be of interest to measure oxygen concentration and isotopic profiles from the surface to well below the OM front for the Clay City site. Such measurements have been done for sandy and sandy clay loam soils in Israel (Angert and others, 2001), where characteristic oxygen penetration depths are well more than a meter. Other profiles with very low oxygen present below a few cm have been measured in rice paddies, where water saturation is high. As oxygen consumption at the OM and pyrite front must create an oxygen concentration gradient to drive diffusive transport, we expect the Clay City O_{2} profile to look similar to figure 18C.

## SUMMARY AND CONCLUSIONS

In order to better understand and model the process of the global oxidation of organic matter, we have modeled the oxidation of organic matter and pyrite in near-surface eroding environments. We conclude (compare to our discussion of fig. 17) that the flux of oxygen consumed by weathering of OM in shales is generally set by the OM content at depth and the erosion rate (as concluded by Holland, 1978). This also implies that an oxygen level dependent feedback for OM weathering is generally not appropriate for mean continental erosion rates and Phanerozoic oxygen levels. Long-term oxygen levels are subject to a complex surface boundary condition, and box models (Lasaga and Ohmoto, 2002) cannot resolve complex boundary layer behavior at the Earth’s surface.

Lasaga and Ohmoto (2002) assume that quadrupling the oxygen concentration in the atmosphere would double the flux drawn down to consume ancient OM (the square-root dependence in their eq 32 or 37). From our calculations (compare figs. 16 and 17) based on the same kinetics for OM oxidation (Chang and Berner, 1999), it is clear that for average erosion rates, and for atmospheric oxygen concentrations of more than 5 percent, increasing the oxygen level has no influence on the amount of oxygen consumed by the oxidation of ancient OM. Lasaga and Ohmoto also conclude that “for the erosion rate to be such that oxygen levels make little difference on the oxidation flux” and for their chosen surface area of OM in rocks, erosion rates would need to be less than 0.8 cm/kyr for *L _{s}* =10 m or less than 0.08 cm/kyr for

*L*=1 m, considerably smaller than average erosion rates (where

_{s}*L*is their assigned characteristic soil depth for a well mixed box where OM is oxidized, compare to fig. 1A). Our generally conservative modeling results would estimate much higher erosion rates for which the oxygen level has little feedback on oxygen consumption flux.

_{s}Some other modeling studies of atmospheric O_{2} and biogeochemical cycles (for example Bergman and others, 2004; Arvidson and Mackenzie, 2006) have not taken into account our finding of the lack of oxygen concentration feedback on ancient OM weathering, These studies, like that of Lasaga and Ohmoto, include a negative feedback due to O_{2}-dependent weathering which is incorrect unless unusually global rapid erosion rates were present in the past.

Factors that favor little ancient organic matter at the surface also lead to buried fronts. Deeper fronts are encouraged by slow erosion, small (deep) organic matter contents, high air porosity, high atmospheric oxygen levels, low water saturation, and low *n** (tortuosity parameter), as is clear from the scaling shown in figure 15. Small organic grain size and rapid reaction rates have a small effect that also makes fronts slightly deeper. Front widths are generally smaller for rapid reaction rates, and to a lesser degree when grain sizes are small (for OM or pyrite) or when fronts are shallow. It would be of interest to convolve the spatially dependent erosion rates with the ancient OM at depth taking into account how much reduced material could be eroded and rapidly reburied, in order to estimate the drawdown rates of oxygen due to erosion on a global scale.

The inclusion of pyrite oxidation in the model does not reproduce the broad pyrite front observed in the field study. There are a number of reasons that may contribute to this discrepancy. Our simple model did not include aqueous phase nor solute transport: each of which could influence pyrite oxidation. Pyrite oxidizes to dissolved species such as SO_{4}^{2−} and Fe^{2+}, and even OM weathering produces some dissolved OM that will be subject to oxidation in solution. Such dissolved species will travel with moving fluids and will influence, and be subjected to, additional chemical changes in the sediments. Wetting and drying, along with air purging associated with rainfall events could potentially bring some oxygen down to the deeper levels episodically, causing pyrite to oxidize below its main front. The framboidal nature of the pyrite may also have an influence on its oxidation, not captured by our assumed particle size on the larger end of the observed range. The formation of iron oxide rims on the initial pyrite oxidized may inhibit diffusion of oxygen to the pyrite rich cores. Such a rim could slow the effective kinetics of pyrite oxidation, thus leading to a broader front. It is also possible that the pyrite concentration was simply not constant in the deep sediments, and that some aspects of the pyrite front may reflect initial heterogeneity of the sediments. Fractures in the matrix may also contribute to enhanced oxygen diffusion.

Our simulations, while a marked improvement over simple box models, are still dramatically simplified compared to natural systems. We have attempted to capture what we view to be the major controlling processes governing the oxidation of organic matter in near-surface eroding shales. However, modern organic matter at the surface may use up some oxygen, making the effective oxygen concentration at the upper boundary somewhat lower than the atmospheric value. We have also assumed the oxidizing layer is above the water table. If the unoxidized OM were below the water table, the diffusion of oxygen through the water would be orders of magnitude slower than diffusion through air. We have also assumed that the pore spaces are interconnected. We would encourage direct measurements of the diffusion of oxygen through shale samples in order to determine the tortuosity. Such measurements were beyond the scope of this study. It would also be of interest to measure oxygen profiles in the weathered sections above oxidizing OM.

Future extensions of our simple model should address additional complexities. These include models that include episodic rainfall, soil air purging, and aqueous phase solute transport coupled with kinetically controlled mineral dissolution and precipitation reactions. The influence of both lateral and vertical heterogeneity of physical and chemical properties should also be addressed in more detail (compare with dual porosity models discussed in van Genuchten and Sudicky, 1999). Of course other mineral weathering reactions should also be considered in this effort. Regarding the experimentally derived kinetics, it is critical to re-measure low [O_{2}] aqueous phase kinetics for OM and pyrite. The possible existence of other slower OM kinetics than were found by Chang and Berner (1999) should also be investigated.

## APPENDIX A

### Conversion Factors

We first give conversion factors between mass fractions, and the concentrations *g _{i}* for the cases of

*i*corresponding to pyrite or C in OM. For volume fractions and densities, type

*i*corresponds simply to pyrite or OM (not C in OM).

*φ*_{B,solid}- (1-
*φ*) is the volume fraction of solid in the bulk (cm^{3}solid/cm^{3}bulk). *ρ*_{solid}- 2.5 (g solid/cm
^{3}solid) is taken as the average density of the solid portion of the medium.

The concentration *g _{i}* ((moles C in OM, or moles pyrite)/(cm

^{3}bulk)) is related to

*MF*, the mass fraction of

_{i}*i*in the dry “rock”, by: (A1) Note that the total organic carbon mass fraction (TOC) is just

*MF*for the appropriate

_{i}*i*.

*MW*_{i}- (12.01 g C / mole C) or (119.975 g FeS
_{2}/mole pyrite), where “g” for grams is not to be confused with*g*for concentrations._{i}

To convert between volume fraction and concentration, we use: (A2) where (A3) The 15.171 g OM/(mole OM) comes from the assumption that 1 mole of OM contains 1 mole of C and corresponds to the element ratios (C_{40}H_{48}O_{2}N_{1}S_{1})÷40 →C_{1}H_{1.2}O_{0.05}N_{0.025}S_{0.025}, whose molecular weight is 15.171 g OM/(mole OM) (Petsch and others, 2000; Petsch, ms, 2000).

*φ _{i,solid}* is the volume fraction of particles of type

*i*(OM or pyrite) in the bulk.

To convert between volume fraction and mass fraction, we use: (A4) The factor 0.7916 corresponds to (12.01 g C/(mole OM)) ÷(15.171 g OM/(mole OM)).

We now give the forms of other factors needed to convert units.

*b*_{i}- a conversion factor to adjust units, where

(A5) *α _{i}* is a shape factor discussed above in the “Choice of Parameters” subsection.

*ρ*_{i}- is the density of OM or pyrite (in g/cm
^{3}). We used the density of graphite for OM (2.26 g/cm^{3}) as did Lasaga and Ohmoto (2002). For pyrite, we used Robie and others, (1978) (5.011 g/cm^{3}). *N̄*_{i}*φ*/(_{i,solid}*α*) is the nucleation density of particle type_{i}^{2}d_{i}^{3}*i*in a 1 m^{3}reference volume, per volume (with units of number per m^{3}) in terms of the volume fraction of OM or pyrite.

The conversion factor in the reaction term is: (A6) where *ρ _{i}* is in cgs units, and the m/cm conversion is required by the fact that

*d*is in units of meters.

_{i}For cases using power-law kinetics we fit the form *A _{pow:Texp}* *([O

_{2}])

*to the experimental data. The parameter*

^{n}*A*is used for the power-law prefactor at the temperature of the experimental measurements, and has the units appropriate to make

_{pow:Texp}*ℜ*with units described in the Mathematical Formulation subsection. Given the power-law form of the rate law, the units are a bit complicated and raised to noninteger powers

_{i}*n*. We used experimental data for oxidation rates that are in units of [(mole pyrite or C in OM consumed)/(m

^{2}exposed pyrite or OM surface area)/s] and experimental measurements of the oxygen concentration in micromole of O

_{2}/Liter of water: ([O

_{2}]). This must be subsequently converted into different units given the form and units of

*A*in our governing equations and the fact that the experimental data used oxygen concentrations in water, but our final power-law form uses oxygen concentrations in the gas phase, so the Henry law partitioning must also be taken into account, through “convptc” (defined below). This leads to the following:

_{i}

*A*_{i}*A**(3.15×10_{pow:Texp}^{7}s/yr) * eactfact_{i}**c**(convptc)_{i}is a unit conversion factor for power-law kinetics used to get the proper units for^{n}*A*. The factor eactfact_{i}_{i}used to convert the rate from the temperature of the experiment to the temperature of the simulation and is defined in the “Choice of Parameters” subsection.

For Michaelis-Menten kinetics, the parameter *K* is the oxygen concentration for a rate with half the maximum value. If the parameter *K _{wat}* were found for the oxygen concentration ([O2]) in [micromole of O

_{2}/Liter of water], a conversion is necessary to find the

*K*for the form of equation (3), where oxygen concentrations

_{m}*a*are in [mole O

_{i}_{2}/cm

^{3}air].

Regarding the conversion of *K _{wat}* to

*K*we have:

_{m}*K*=

_{m}*K*/convptc with (A7) The mole fraction of O

_{wat}_{2}in water is

*X**(partial pressure of O

_{1}_{2}in air, which for a total pressure of 1 atm. is the volume fraction of O

_{2}in air), with

*X*from Gevantman (2006). For a total pressure of 1 atm., the partial pressure of O

_{1}_{2}in air (in atm.) is numerically equal to the volume fraction of O

_{2}in air. M

_{wat}is the molecular weight of water.

## APPENDIX B

### Numerical Methods

We used an implicit finite-difference method to solve the equations presented above. Some special difficulties arise in the equations given above due to the “stiff” nature of the equations (compare Press and others, 1986). During a transient phase of the solution, as the system is evolving toward steady state, the reaction terms may be much larger than the advection and diffusion terms. As the system comes close to steady state, a balance develops between advection (erosion bringing up fresh material) and reaction for the pyrite and OM equations, while a balance of reaction with diffusion terms develops in the oxygen equation. For equivalent surface areas and oxygen concentrations, the oxidation rate of pyrite is at least an order of magnitude larger than that of the OM. The numerical solutions obtained for typical OM and pyrite “grain” sizes and shapes leads to “fronts” that are much sharper for pyrite than for OM (see fig. 18). In addition, each of the components of the system may essentially disappear in some locations (oxygen at depth, pyrite and OM near the surface). This latter factor necessitates the use of thresholds in the code: levels of components below which one must assume the component is absent. Some of the favored implicit methods for ordinary differential equations (for example Gear 1971, p. 271) involve solving nonlinear algebraic equations at each grid point and time step, and are not readily adapted to partial differential equations. When one of the components disappears (for example pyrite), the number of unknowns at that grid point would have to be modified. The logistical complications of dealing with disappearing components and changing matrix array sizes forced us to devise a specialized direct iteration implicit method. Care is always needed in adjusting the grid resolution and time step size to assure numerical accuracy.

The method used to solve the finite difference equations is as follows. Each of the terms in the governing equations (1)-(5) may be written in their finite difference form. For a particular vertical level (*x _{j}*) the trial values (

*N**) for the new time level

*N*was written as (1) (2) where

*o*denotes the old time level,

*N*-1. (B3) In the above expression,

*ϖ*is a weighting factor for the time level of the diffusion term, with

_{d}*ϖ*=1 for totally implicit, and

_{d}*ϖ*=0 for an explicit step. We generally chose

_{d}*ϖ*=0.5 or 0.75. The term inside the first [] brackets is the 1D laplacian of

_{d}*a*in central difference form. The term inside the second [] accounts for spatial variations in the effective diffusivity and porosity and is multiplied by the gradient of the oxygen concentration, all of these are in a one-sided backward difference from level

*j*.

The erosion term, at vertical level *j*, is written as: (B4) for each *i*, with a forward difference on *g* (upstream weighting, in the sense that in the frame fixed to the eroding surface, erosion brings fresh material up from below). Here, *ϖ _{adv}* is a weighting factor for the time level estimates of

*g*, with

_{i}*ϖ*=1 for complete weighting on the guess value *, and

_{adv}*ϖ*=0 for a solely explicit estimate using the old,

_{adv}*N*-1, time level. The estimates

*g**differ from the

*g*values in equation B2 as follows. Prior to the time step a value for the advected

^{N*}*g*value at

*x*is estimated by averaging a time projection value (

_{j}*g*=

_{i,j}^{TP}*ϖ*(1.5

_{TP}*g*-0.5

_{i,j}^{o}*g*)+ (1-

_{i,j}^{N−2}*ϖ*) (

_{TP}*g*)), with

_{i,j}^{o}*N*−2 indicating an older time level, and a spatial projection value (

*g*=

_{i,j}^{SP}*ϖ*(

_{SP}*g*-(

_{i,j}^{o}*ωΔt/Δx*)(

*g*−

_{i,j+1}^{o}*g*)+ (1-

_{i,j}^{o}*ϖ*) (

_{SP}*g*))) to yield

_{i,j}^{o}*g**=

_{i,j}*ϖ*(

_{xot}*g*)

_{i,j}^{SP}*+*(

*1- ϖ*)

_{xot}*g*. This combination of spatial and temporal projections helps accelerate the convergence of the iteration within the time step. Similar estimates were made for

_{i,j}^{TP}*d**for use in the reaction term. Within each time step, iterations are performed. After the first iteration within a time step, the next value used for

*g**is simply

_{i,j}*g*as found by the previous sub-time-step iteration. Generally two such sub-time-step iterations were performed.

_{i,j}^{N*}Special care is taken with the reaction term to avoid underflow and overflow when values of *a, d*, or *g* become very small. Reaction terms are zeroed out if the size of the grains had *d _{i,j}* < 1 nm or if

*a*< 10

_{j}^{−40}mole O

_{2}/cm

^{3}air. In the code, both

*g*and

_{i}*d*are dynamic. The equation governing

_{i}*g*implied changes in

_{i}*d*and thresholds for

_{i}*d*also imply thresholds for

_{i}*g*. The

_{i}*g*

^{2/3}form (eq 6) was generally used for

*g/d*in the reaction terms (eq 3 or eq 4). Some runs were done with a special estimate for denominator values of

*d*of the form

*d̄*

*=d*ϖ*after the first sub-iteration, and similarly for

_{bar}+ d^{o}(1-ϖ_{bar})*a*in the denominator for Michaelis-Menten kinetics. For the first sub-iteration, the value of

*d**described previously was used for the denominator value. For the case of power law kinetics,

*a*in the reaction term was partitioned into old and estimated values as: (B5) This form was derived from the Taylor expansion about the old value of

*a*for small changes in

*a*with a time projection of ϖ

_{bar}(0 for explicit and 1 for full time step projection). The old time step (

*N*-1) values were used for the porosity and the effective diffusivity in the above expressions. Special forms of the finite difference expressions given above are required at the top and bottom boundaries of the 1D computational domain, to accommodate boundary conditions. The equation for

*a*was solved by transferring the

*a*terms to the left-hand side, while values depending on

^{N*}*a**or

*a*were placed on the right-hand side of a system that was solved by a tridiagonal algorithm. After the

^{o}*a*equation was solved, the

*g*equations were solved using the old and new time-step values of

_{i}*a*for the reaction terms in an identical form as used in the

*a*equation (with proper adjustment for the stoichiometry and porosity).

## Acknowledgments

We would like to thank Neil Blair, Richard Wildman, and an anonymous referee for reviewing the manuscript. We acknowledge Richard Wildman, James O. Eckert, Ulrich Mok, and J. Brian Evans for technical help supporting this and the previous paper. We thank Jay J. Ague, Mark Brandon, Ron Smith, and Ken Caldeira for helpful discussions. EWB would like to thank David Rossman (for a clever perl program used to search plot files) and Magdalena Gleisner (for supplying the numerical values for the data of Gleisner and others, 2004). This research was primarily supported by Department of Energy (Chemical Sciences, Geosciences and Biosciences Division, Office of Basic Energy Sciences, Office of Science) Grants and DE-FG02-01ER15216 to D. M. Rye, J. J. Ague, and E. W. Bolton; and DE-FG02-01ER15173 to R. A. Berner; National Science Foundation Grant EAR 0104797 to R. A. Berner; and Grant NSF-EAR 06160707 to S. T. Petsch, T. I. Eglinton and K. J. Edwards. EWB received additional support from Yale University, NASA grant NNG05GQ97G to Rob Rye and E.W. Bolton, and some of this work was performed as part of the NASA Astrobiology Institute’s Virtual Planetary Laboratory Lead Team, supported by the National Aeronautics and Space Administration through the NASA Astrobiology Institute under Cooperative Agreement No. CAN-00-OSS-01.