Triple oxygen isotope constraints on atmospheric O2 and biological productivity during the mid-Proterozoic.

Reconstructing the history of biological productivity and atmospheric oxygen partial pressure (pO2) is a fundamental goal of geobiology. Recently, the mass-independent fractionation of oxygen isotopes (O-MIF) has been used as a tool for estimating pO2 and productivity during the Proterozoic. O-MIF, reported as Δ'17O, is produced during the formation of ozone and destroyed by isotopic exchange with water by biological and chemical processes. Atmospheric O-MIF can be preserved in the geologic record when pyrite (FeS2) is oxidized during weathering, and the sulfur is redeposited as sulfate. Here, sedimentary sulfates from the ∼1.4-Ga Sibley Formation are reanalyzed using a detailed one-dimensional photochemical model that includes physical constraints on air-sea gas exchange. Previous analyses of these data concluded that pO2 at that time was <1% PAL (times the present atmospheric level). Our model shows that the upper limit on pO2 is essentially unconstrained by these data. Indeed, pO2 levels below 0.8% PAL are possible only if atmospheric methane was more abundant than today (so that pCO2 could have been lower) or if the Sibley O-MIF data were diluted by reprocessing before the sulfates were deposited. Our model also shows that, contrary to previous assertions, marine productivity cannot be reliably constrained by the O-MIF data because the exchange of molecular oxygen (O2) between the atmosphere and surface ocean is controlled more by air-sea gas transfer rates than by biological productivity. Improved estimates of pCO2 and/or improved proxies for Δ'17O of atmospheric O2 would allow tighter constraints to be placed on mid-Proterozoic pO2.

The oxygenation of Earth’s atmosphere–ocean system is widely considered to have occurred through two major step function increases in atmospheric oxygen (O2) (1, 2). In the first of these steps, referred to as the “Great Oxidation Event,” Earth’s atmosphere transitioned from a weakly reducing state in which O2 was poorly mixed and vanishingly low in abundance (3, 4) to a state in which atmospheric O2 was present in sufficient abundance to prevent the production and transfer of mass-independent fractionation (MIF) of sulfur isotopes to marine sediments (5). Photochemical models indicate that the disappearance of these signals from the rock record at ∼2.3 billion years ago (Ga) (6) requires atmospheric O2 levels above ∼10−5 times the present atmospheric level (PAL) (7). A second increase in ocean–atmosphere O2 levels is now considered to have occurred at some time between ∼800 to 400 Ma, during which time atmospheric O2 began to approach roughly modern values, and the ocean interior became periodically and then more or less permanently well oxygenated (2, 8, 9).

However, atmospheric O2 levels in between these two broad steps in Earth system oxygenation remain debated and poorly constrained. Evidence for widespread anoxia in the ocean interior has been taken to suggest atmospheric O2 levels below ∼40% PAL (1, 10). Estimates based on stable chromium (Cr) isotopes in shales and oolitic ironstones (11, 12), cerium (Ce) systematics in marine carbonates (13, 14), and iron (Fe) retention patterns in ancient soil horizons (15) are all consistent with very low atmospheric O2, perhaps below ∼1% PAL. In contrast, trace element abundances and isotopic signatures in shales (16–18) and iodine systematics of marine carbonate concretions (19) have been used to suggest higher atmospheric O2 levels, well above ∼1% PAL and perhaps even >10% PAL. A reanalysis of the Ce anomaly data also suggests atmospheric oxygen partial pressure (pO2) > 1% PAL (20). Constraining atmospheric O2 prior to the emergence and ecological expansion of metazoan (animal) life is essential for understanding the relative roles of intrinsic biological factors and environmental drivers (or, more likely, their interplay) in controlling the expansion of biological complexity and, in particular, the rise of animals (21–23). Understanding how primary productivity has changed through time as atmospheric O2 rose is also a fundamental unanswered question about Earth’s history (e.g., refs. 24 to 26).

Here, we employ a tool for estimating ancient pO2 and primary productivity in deep time: the triple oxygen isotope composition of marine and lacustrine sedimentary rocks. This approach utilizes oxygen isotope signatures that deviate significantly from the usual mass-dependent fractionation (MDF) line. Today, about 30% of this oxygen MIF (O-MIF) signal originates from ozone (O3) photochemistry, and we expect that this percentage was considerably higher in the past (). The O-MIF value can be defined by the relationship (27)

Here, δxO represents the deviation in parts per thousand (‰, or “permil”) from the xO/16O ratio in an isotopic standard, and λ is a reference line slope that we take to be equal to 0.528, following Young et al. (28) and Pack (27). For small fractionations, this reduces to

We use the logarithmic form (Eq. in all calculations, as the fractionation in ozone is not small. The positive Δ′17O signal from ozone (∼26‰ at the surface) is transferred to CO2, and O2 acquires a negative Δ′17O signal (∼−0.43‰ at the surface) by mass balance (27). The Δ′17O signal from atmospheric O2 can be transferred to the rock record via weathering of pyrite, FeS2, and the subsequent burial of oxidized sulfur phases at Earth’s surface (25, 29, 30). Some authors have recently suggested that H2O2 may be responsible for oxidizing pyrite in rivers draining the Himalayas (31), so we keep track of Δ′17O of H2O2 as well. Small changes in Δ′17O can also be induced by mass-dependent processes with different values for λ as discussed further below (see the section Dilution of the O-MIF Signal in O). We use the term “O-MIF” only to refer to the large changes induced by ozone formation.

On the present Earth, Δ′17O values are small, roughly −0.05‰ in seawater sulfate (32, 33) and −0.43‰ in tropospheric O2 (27). By comparison, the existing oxygen isotope record from sedimentary sulfate minerals deposited during the Proterozoic shows values as light as −0.9‰ (25), requiring tropospheric O2 values significantly lighter than those of the modern atmosphere. The overall magnitude of these signals is controlled by atmospheric O2 and CO2, which act to control the production and storage of isotopic anomalies in the atmosphere and by the globally integrated productivity of the biosphere, which acts to eliminate MIF (34, 35). As a result, if one knows (or can assume) atmospheric pCO2, one can either use an assumed atmospheric pO2 to estimate global gross primary production, GPP (25, 36), or invert for an estimate of atmospheric pO2 (15, 37). (GPP is the total rate of synthesis of fresh organic matter, which, for an aerobic ecosystem, is just the rate of oxygenic photosynthesis.) However, previous attempts to solve for pO2 (ibid.) used a box modeling approach that mischaracterized the effect of stratosphere–troposphere gas exchange at low pO2, made potentially problematic assumptions about the link between GPP and air–sea O2 transfer, and left out other key aspects of atmospheric chemistry.

Here, we revisit this problem with a detailed one-dimensional photochemical and ocean–atmosphere gas exchange model. This approach includes several advances that allow for a more realistic representation of the oxygen cycle. First, our photochemical model explicitly calculates the lowering of the ozone layer as atmospheric O2 decreases (38) and directly simulates the photochemically catalyzed isotopic exchange between O2 (or O3) and gaseous H2O. Second, our model includes limitations imposed by gas exchange rates across the air–sea interface. At O2 levels lower than today, the gas transfer rate, rather than GPP itself, should control the extent of biospheric recycling of Δ′17O values.

Production of the O-MIF (Δ′17O) Signal during Ozone Formation

O-MIF in Earth’s atmosphere is mainly produced by the ozone formation reaction (39–41)

Here, “O” is the ground state of atomic oxygen, and “M” represents a third molecule needed to carry off the excess energy of the collision. Ozone is a bent molecule with a resonant double bond between three oxygen atoms (). If we let “Q” represent a minor isotope of O (i.e., 17O or 18O), the isotopic analogs of reaction [] can be written asand

Two different forms of ozone are created: asymmetric OOQ (reactions [] and []) and symmetric OQO (reaction []). The rate constants for these reactions are given in . Asymmetric ozone in reaction [] forms ∼16% faster than symmetric ozone in reaction [] (the so-called “ƞ effect” in the chemical physics literature (40, 42–44). This effect depends on both temperature (45) and pressure (46). Reaction [] forms asymmetric ozone as well, but its rate constant is roughly equal to that of reaction [] (for 17O), or even slower (for 18O), for reasons that are related to the rapid isotope exchange reaction

The rate of isotope exchange as well as the rate of ozone formation is influenced by the symmetry number of two for the O2 molecule and by the vibrational zero-point energy difference between lighter OO and heavier OQ. The net result is that O3 becomes enriched in minor isotopes (Q atoms) by about 70‰ (17O) and 75‰ (18O), producing Δ′17O values of ∼30‰ following Eq. . This large enrichment has been confirmed by stratospheric observations made from balloon flights (47–49). The enrichment is then transferred to CO2 via the reaction sequences (50)and

Here, hv represents an ultraviolet photon, and Q(1D) is an excited electronic state of atomic oxygen. Thus, atmospheric CO2 acquires positive Δ′17O inherited from ozone. The excess Q atoms in CO2 are ultimately derived from O2; hence, O2 acquires negative Δ′17O. However, because O2 is ∼500 times more abundant than CO2 in the modern atmosphere, the negative Δ′17O values in O2 are much smaller than the positive Δ′17O values in CO2. Q(1D) can also react with other species (e.g., H2O and CH4), thereby spreading the O-MIF signal throughout the photochemical network.

Dilution of the O-MIF Signal in O2 and CO2 via Biospheric Recycling

Recycling by the Marine Biosphere.

The O-MIF signal in O2 and CO2 can be diluted by processing through the terrestrial and marine biospheres (35). Consider the marine biosphere first. CO2 dissociates in solution to form carbonic acid and other forms of dissolved inorganic carbon (DIC), all of which readily exchange O atoms with water. [The time scale for isotopic exchange is essentially equal to the time scale for CO2 hydration, which is of the order of minutes (51). By comparison, the residence time of CO2 in the surface ocean is ∼20 d, as shown in . (CO2 and O2 behave similarly because their air-sea transfer veocities are nearly the same.) Thus, whenever CO2 dissolves in water, it loses its O-MIF signal. However, the situation is quite different for O2. Dissolved oxygen does not exchange O atoms with water unless it is used for respiration (35), which can be represented by reaction [] running to the right

The reverse of respiration is photosynthesis, which is represented by reaction [] running to the left.

Respiration dilutes the O-MIF signal in O2 by reducing much of it to water and adding some O atoms to other molecules. However, it also creates a small Δ′17O signal by way of a phenomenon called the “Dole effect.” Respiration induces strong MDF because lighter O isotopes are taken up faster than heavier ones (52). This causes atmospheric O2 to be enriched in minor isotopes by ∼23.5‰ for 18O and ∼12‰ for 17O (28). The slope of the MDF line for respiration is ∼0.515 (28), which is slightly lower than the slope of the normal terrestrial MDF line given by Eq. . This difference imparts a negative Δ′17O signal to O2, which is important when we compare our model to observed values in Fits to Modern Earth. This process is less of a complication at low pO2 levels because air–sea gas transfer rates decrease while the O-MIF in O2 signal becomes substantially larger; hence, we have elected to ignore it in our low-O2 calculations. We justify this assumption explicitly in .

The processing of the O-MIF signal by the marine biosphere is complicated because air–sea gas transfer limitations also play a role. The net downward flux of gas X across the atmosphere–ocean interface can be expressed as a “piston velocity,” kX, times a concentration difference over a thin boundary layer at the ocean’s surface:

The piston velocity can be related to the product of the 10-m average wind speed and the Schmidt number of the fluid (53) (). The remaining terms in Eq. are the Henry’s Law coefficient αx and the species atmospheric partial pressure pX. [X] is the dissolved concentration of the gas in the surface ocean. We found that, ≅ 5.0 m ⋅ d−1, and ≅ 4.7 m ⋅ d−1.

Note that Φ(X) in Eq. is a net flux, the difference between a gross downward flux proportional to αX⋅pX and a gross upward flux proportional to [X]. The upward and downward fluxes are independent (i.e., they do not impede each other). Thus, for CO2, which exchanges O atoms upon contact with water, the dilution term is straightforward: it is equal to the gross downward flux ⋅⋅pCO2, which can be used to calculate a maximum deposition velocity in the photochemical model (). Thus, for CO2, the dilution term is entirely independent of GPP. We explain how this is implemented in Methods.

For O2, the situation is more complicated because isotopes are not exchanged unless the O2 is taken up by respiration, so one cannot simply use the gross upward and downward fluxes implied by Eq. . The relevant air–sea exchange flux should include only the O2 that is reprocessed by respiration and photosynthesis. Furthermore, O2 that is produced by photosynthesis within the surface ocean and then consumed in situ by respiration should not be counted as contributing to this exchange flux. An additional consideration is that the time-averaged net O2 flux at the atmosphere–ocean interface must be approximately zero if the system is to remain in steady state. (We can neglect the ∼10 Tmol ⋅ yr−1 of O2 that is consumed by oxidative weathering on the continents [table 10.1 in ref. 54.]) Thus, the surface ocean must remain near Henry’s Law equilibrium with the atmosphere on average according to Eq. . However, in reality, the surface ocean departs substantially from equilibrium at different localities and at different seasons, driving measurable atmosphere–ocean O2 exchange (55).

We can account for these complex O2 exchange processes in our globally averaged model by following an approach pioneered by Bender et al. (52). On the modern Earth, marine GPP produces ∼12,000 Tmol ⋅ O2 ⋅ yr−1 (52), a value that falls nicely in between more recent estimates of 8,000 to 17,000 Tmol ⋅ yr−1 (56). We use the older values here because they have been used by other authors (25, 34, 35, 37), so one can compare with their results. Bender et al. (52) estimate that 12% of marine GPP is recycled internally within the surface ocean, so the relevant sea–air (or air–sea) flux of O2 is 12,000 × 0.88 = 10,600 Tmol ⋅ O2 ⋅ yr−1.

To estimate the recycled O2 fraction at lower pO2 levels, one must understand the logic by which Bender et al. (52) calculated their 12% modern value. They compared the timescale for air–sea exchange of O2 with that for uptake by respiration. The timescale for air–sea exchange is the dissolved O2 reservoir size in the ocean mixed layer divided by the (gross) upward O2 flux from Eq. . Both quantities are proportional to pO2 (at least on a time average), so the exchange timescale is just a constant equal to the depth of the ocean mixed layer divided by the piston velocity. The timescale for uptake by respiration is equal to the dissolved O2 reservoir divided by GPP, so it depends linearly on pO2 if GPP remains constant. Today, the timescale for air–sea exchange is ∼1/14 that for respiration, so most photosynthetically generated O2 escapes to the atmosphere. At lower O2 levels, this ratio decreases, so more photosynthetic O2 is recycled. The calculation is shown explicitly in . The results are shown in Fig. 1. Importantly, GPP is assumed to remain constant all the way down to 0.01 PAL O2, thereby maximizing dilution of the Δ′17O signal.

Fig. 1.

Atmosphere–ocean O2 exchange rates for different O2 mixing ratios. The red curve represents the PV limit described in Recycling by the Marine Biosphere. The orange curve shows the calculated O2 flux from . The solid portions of the red and orange curves show the air–sea O2 exchange rates employed in the low-O2 calculations described in the text. The line labeled “Estimated terrestrial mat GPP” shows the effective O2 exchange rate with the Proterozoic terrestrial biosphere in the Fig. 4 simulations. The dashed grey line is the effective PV for O2 exchange between water-saturated mats and the atmosphere as described in Recycling by the Terrestrial Biosphere.

Fig. 4.

Calculated Δ′17O values for O2 in the marine biosphere–only simulations (A) and the terrestrial biosphere–included simulations (B). The red-contoured interval represents the range of values that are consistent with the 1.4-Ga Sibley sulfate data, assuming 8 to 15% incorporation of O2. The dashed lines represent various suggested constraints on pCO2 discussed in the text.

At pO2 levels of 1% PAL and below, the O2 exchange flux calculated by this method approaches the “piston velocity (PV) limit,” so we simply use the latter value as we do for CO2. Again, this maximizes the rate of air–sea O2 transfer, and thus minimizes the Δ′17O signal in O2. The limiting O2 flux, Φ, is readily calculated by setting the dissolved O2 concentration equal to zero in Eq. : Φ ⋅⋅pO2 ≅ 6.5 × (pO2/1 bar) moles ⋅ m−2 ⋅ d−1, using the O2 solubility from . Converting time units to years and multiplying by the area of the surface ocean, 5.1 ×1014 m2, yields the PV limit values shown in and in Fig. 1. More details about model implementation are shown in Methods.

Critically, our estimated air–sea O2 exchange fluxes at low pO2 levels are only weakly connected to marine GPP. The exchange flux could be lower than calculated here if Proterozoic GPP was substantially lower than today, but it cannot be higher unless Proterozoic GPP was higher than today, and pO2 was also relatively high. This same reasoning implies that the O-MIF data from sedimentary sulfates are largely independent of marine GPP. Thus, previous attempts to estimate Proterozoic GPP from such data (e.g., ref. 25) need to be rethought.

Recycling by the Terrestrial Biosphere.

According to Farquhar et al. (57) and Bender et al. (52), modern terrestrial GPP is ∼14,100 Tmol ⋅ yr−1. We again adopt values from Bender et al. (52), which remain close to more recent estimates: 12,250 Tmol ⋅ yr−1 (58) and 13,900 Tmol ⋅ yr−1 (59). However, photorespiration, which recycles O2 without contributing to GPP, is also important. (In photorespiration, the enzyme RuBisCO oxygenates RuBP instead of carboxylating it, thereby short-circuiting the process of photosynthesis.) Thus, terrestrial uptake of O2 rises to 20,400 Tmol ⋅ yr−1 when this process is considered (52). CO2 is not produced by photorespiration, so the recycling of CO2 by land plants should go at the slower, uncorrected GPP rate. Isotope exchange between CO2 and plant water is complicated by the fact that the residence time of CO2 in leaves is often shorter than the time scale (minutes) for CO2 hydration. However, this problem is alleviated by the ubiquitous presence in plants of the enzyme carbonic anhydrase, which greatly accelerates hydration (51). Therefore, we assume that the Δ′17O signal in CO2 is almost zeroed out by this process as it is by recycling through the marine biosphere.

The processing of atmospheric O2 by the terrestrial biosphere should have been slower during the mid-Proterozoic because vascular plants had not yet evolved. That said, significant portions of the exposed land surface may have been covered by microbial mats composed largely or partly of O2-generating cyanobacteria (37, 60). Estimating the GPP of such mats is problematic because it is difficult to determine what fraction of the land surface should have been covered and because the area of the continents may have changed over time. Hence, we present calculations both with and without terrestrial mats. Based on a detailed study of a modern microbial mat in Indonesia, Finke et al. (60) estimated a terrestrial microbial mat GPP of 20 to 200 Tmol ⋅ O2 ⋅ yr−1 for an Archean Earth with 1/6 the current land area. This translates to 120 to 1,200 Tmol ⋅ O2 ⋅ yr−1 if the continents had grown to their present size by 1.4 Ga. Other authors have estimated Proterozoic microbial mat fluxes as high as 4,000 Tmol ⋅ O2 ⋅ yr−1 (26). We use a terrestrial O2 flux of 660 Tmol ⋅ yr−1 (∼5% of modern terrestrial GPP) in the calculations discussed in the next paragraph and shown in Fig. 4B. This number is close to the midrange terrestrial O2 flux in Planavsky et al. (26).

Terrestrial microbial mats are subject to gas exchange constraints analogous to those imposed at the atmosphere–ocean interface. The surface of a water-saturated mat is much less permeable, though, than the turbulent ocean surface. Based on observations of the dissolved O2 profile within a modern, water-saturated microbial mat by Finke et al. (60), we estimate that the effective PV for O2 exchange between the mat and the atmosphere is about 1/60 that for atmosphere–ocean exchange (). The O2 storage capacity for a mat is also much smaller than that of the surface ocean, so more of its photosynthetically produced O2, ∼70%, should be released to the atmosphere under low-pO2 conditions (60). Because of the low exchange velocity, this O2 could not have easily flowed back into the mat if pO2 was low. Thus, to balance O2 release, productive mats must also have released CH4 (61) or H2/CO (62) in stoichiometrically equivalent amounts; otherwise, that O2 would have accumulated in the atmosphere. Here, we assume that the O2 flux from mats was balanced by CH4. The 660 Tmol(O2) ⋅ yr−1 produced in our simulations that include mats corresponds to a CH4 flux of 330 Tmol ⋅ yr−1, or about 10 times the modern CH4 flux ().

Results

Fits to Modern Earth.

We implemented the chemistry and transport described in Production of the O-MIF (Δ′ and Dilution of the O-MIF Signal in O in a recent version of our one-dimensional photochemical (main) model (63), now linked to an isotopic model which duplicates this same chemistry for species containing one minor O isotope, either 17O or 18O. O2 and CO2 were treated as variables in both the main code and the isotopic code, but their boundary conditions were handled differently (see Methods). For our modern atmosphere calculations, the surface O2 mixing ratio was set to 0.21 in the main code, and the CO2 mixing ratio was set to 300 ppmv. We henceforth define these concentrations as 1 PAL. A relatively low value of pCO2, near the preindustrial level of 280 ppmv, is appropriate for these calculations because the Δ′17O signal in modern atmospheric O2 has been gradually accumulating on a timescale commensurate with the lifetime of O2 against isotopic recycling, ∼1,200 y (35). In both the main code and the isotopic code, O2 and CO2 were given biospheric recycling fluxes using numbers from Dilution of the O-MIF Signal in O. The O2 flux from the combined marine and terrestrial biosphere was set equal to 31,000 Tmol ⋅ yr−1. The CO2 flux was 22,500 Tmol ⋅ yr−1, which represents the sum of a gross marine flux calculated from the downward part of Eq. along with a flux equal to terrestrial GPP.

Model predictions for Δ′17O of O2, O3, and CO2 are compared with measurements in Fig. 2. The stratospheric ozone measurements are from balloon flights carried out over the last ∼40 y (see ref. 49 and references therein). They measure bulk ozone, separated cryogenically from air and converted quantitatively to O2. Isotopic measurements of surface ozone were made indirectly by oxidizing nitrite and measuring the composition of the resulting nitrate (64). This measures only the transferrable (end atom) isotopic composition of O3. Symmetric ozone (OQO) was assumed by those authors to have Δ′17O = 0‰. The fractionation values in the ozone data are affected by local meteorological factors such that the records from different places and times vary significantly. The shape of the modeled Δ′17O profile for ozone reflects temperature- and pressure-dependent correction factors described in . The US 1976 Standard Atmosphere temperature profile was used in all calculations (), including those at lower pO2 levels. Sensitivity studies were performed to determine the effect on our results ().

Fig. 2.

Comparison between model-generated Δ′17O values for our present-day simulation (curves) and atmospheric observations (symbols). Filled blue squares, blue circles, and the blue star are O2 measurements for different altitudes from Pack et al. (67), Thiemens et al. (68), and Pack (27), respectively. The blue curves are model-predicted values with (solid) and without (dotted) corrections for the Dole effect. Measurements for CO2 are from Kawagucci et al. (66), with filled black circles showing measurements from Kiruna, Sweden (68°N) on February 22, 1997, and open black circles showing measurements from Sanriku, Japan (39°N) on August 31, 1994. Measurements for O3 are from Krankowsky et al. (49), with filled green circles, open green circles, and filled green squares showing measurements from Brazil, France, and Sweden, respectively. The green star is the measured surface fractionation of bulk O3 from Vicars and Savarino (64). The light green curve shows the Δ′17O of stratospheric O3 calculated by Liang et al. (65). Note the scale change at Δ′17O = 0‰.

Importantly, to generate Fig. 2, the rate constant for the formation of asymmetric 49O3 by reaction [] has been adjusted downward by a factor β = 0.9772 (). Without this correction factor, our model predicts Δ′17O values of ∼37‰ for O3 at 25 km altitude, which is well above the average of the data. With the correction, Δ′17O is ∼28‰ at that altitude, which is within the range of observations; however, we underpredict Δ′17O of O3 above 30 km. The correction also brings us closer to the tropospheric ozone data (see Fig. 2). We feel justified in making this correction because we are ultimately interested in Δ′17O of O2, not O3. Overpredicting the (positive) Δ′17O of O3 would cause us to also overpredict the magnitude of the (negative) Δ′17O value of O2. Young et al. (28) employed a similar strategy because their model also overpredicted Δ′17O of O3 and because they were also focused on O2. The model of Liang et al. (65) overpredicted Δ′17O of stratospheric O3 as well (Fig. 2). Even with this correction factor, our model still slightly overpredicts Δ′17O of surface ozone: we get ∼27.5‰, whereas the reported value is ∼26‰ (64). The actual value of Δ′17O for surface ozone depends strongly on the Δ′17O of symmetric O3, which has not been measured experimentally. The discrepancy between calculation and measurement will only be resolved when we fully understand the isotopic fractionation in reactions [4-6].

Predicted Δ′17O values for CO2 increase from 0‰ at the surface to ∼10‰ at 40 km, in agreement with the data from Sweden but not as close to the data from Japan (66). CO2 is expected to be unfractionated near the surface because of rapid exchange of O isotopes with seawater. Comparisons of Δ′17O between model and observations for H2O2 and nitrate are shown in .

For O2, the Δ′17O value predicted near the surface by our model, −0.14‰, is well above observed values. The O2 data are scattered, but Pack (27) suggests that −0.43‰ is a reasonable average value. This discrepancy is expected because according to Young et al. (28), 2/3 of the isotopic signal in modern O2 is caused by the Dole effect, described in Dilution of the O-MIF Signal in O. Accordingly, we subtract 0.29‰ from our calculated values to produce the solid curve in Fig. 2, which agrees well with the observations (27, 67, 68). Recall that this correction is expected to be smaller at lower pO2 levels (Dilution of the O-MIF Signal in O).

Results for Lower pO2 Levels (Marine Biosphere Only).

In a second set of calculations, we lowered atmospheric pO2 from 1 to 10−3 PAL, keeping pCO2 constant at 1 PAL. The interpretation of model results becomes complicated at pO2 < 10−3 PAL because terms other than air–sea exchange become important in the O2 budget. In these simulations and in the ones described in the next section, we set the O2 surface exchange flux equal to the values shown in using the PV-limited flux when pO2 is lower than 0.03 PAL. Surface fluxes of other trace gases (CH4, N2O, CO, and H2) were kept equal to their modern values (see Methods). The selected results are shown in Fig. 3. Notably, as O2 decreases, the ozone layer moves downward, and its peak density decreases (Fig. 3), consistent with earlier results (38, 69, 70).

Fig. 3.

Vertical profiles of (A) ozone (O3) number density, (B) Δ′17O of O3, (C) Δ′17O of O2, and (D) Δ′17O of CO2 from photochemical model across a range of assumed ground-level atmospheric pO2 values (in PAL). All results assume an atmospheric pCO2 of 300 ppmv.

Because the CO2 mixing ratio is held constant in these calculations, the O2:CO2 ratio decreases at lower O2 levels, causing O2 to become more negatively fractionated (e.g., the Δ′17O in surface O2 is about −8.5‰ at 10−3 PAL pO2) (Fig. 3). This depletion in Δ′17O in O2 can be transferred to O3, so the fractionation in ozone also decreases with decreasing O2. For example, the surface ozone fractionation decreases from 27.5 to 19.3‰ as O2 decreases from 1 to 10−3 PAL. Tropospheric CO2 stays near 0‰ because isotopic exchange with seawater remains rapid in all these calculations. The large decrease in Δ′17O for tropospheric O2 at low O2 levels is a result of mass balance. In today’s atmosphere, the positive Δ′17O signal in O3 is transferred to CO2 in the stratosphere as discussed in Production of the O-MIF (Δ′. CO2 then flows downward into the troposphere and ultimately transfers that signal to seawater. Residual O2 becomes negative in Δ′17O as a result. However, at low pO2, the ozone layer moves lower in the stratosphere or even into the troposphere (38, 70). CO2 remains neutral or even slightly negative in Δ′17O (Fig. 3), while ozone transfers its positive Δ′17O signal to water vapor—and from there to the ocean—by way of the HOx chemistry described in . This latter process is highly efficient, and the atmospheric O2 reservoir is small, so O2 becomes strongly negative in Δ′17O.

Next, we repeated our calculations over a grid with pO2 ranging from 1 to 10−3 PAL and with pCO2 from 1 to 300 PAL, yielding 42 assumed pO2/pCO2 combinations. Following the methodology described earlier in this section and in Methods, CO2 surface exchange fluxes remain at the PV limit, so they become large at high pCO2 values. The O2 exchange fluxes were taken from . The calculated Δ′17O values for ground-level O2 are shown in Fig. 4. Some authors (e.g., ref. 31) have suggested that the main oxidizing agent for pyrite is H2O2, not O2. The Δ′17O value of H2O2 tracks that of O2 closely (), however, suggesting that it may not matter which species actually causes the oxidation.

When pyrite is oxidized to sulfate by O2, only 8 to 15% of the oxygen in the sulfate comes from O2; the rest comes from water (71). By contrast, Killingsworth et al. (72) suggest that the incorporated O2 fraction is 18 ± 9%. We use the lower numbers from Balci et al. (71) here to remain consistent with previous studies (25, 37). Higher O2 incorporation in sulfate would lead to higher predicted pO2 values. The 1.4-Ga Sibley sulfates were deposited in a shallow evaporative setting (25). The most negative Δ′17O value observed in these sediments is −0.9‰. Indeed, all these data are below −0.3‰, below the range of Phanerozoic sulfates (25). The O2 needed to produce this fractionation must therefore have Δ′17O values in the range of −6 to −11‰. Most of this O2 is taken up in the last oxidation step going from sulfite to sulfate (71), so we assume that a similar uptake fraction would apply to H2O2 (). This range is represented by the red-contoured interval in Fig. 4. The Δ′17O values higher than −6‰ in O2 are excluded by the Sibley data, assuming that the Balci et al. (71) oxidant incorporation factors are correct. The Δ′17O values lower than −11‰ in O2 cannot be excluded, but they would require bacterial reprocessing of dissolved sulfate and consequent dilution of the O-MIF signal within the water column beneath which the sulfates were deposited. We ignore the Dole effect correction to Δ′17O of O2 at pO2 levels below 1 PAL. More than half of this correction comes from the terrestrial biosphere, which should have been greatly diminished at 1.4 Ga if pO2 was low (Recycling by the Terrestrial Biosphere). Any Dole effect caused by the marine biosphere at that time would push our calculated Δ′17O values lower for O2, causing our predicted pO2 values to be higher than those estimated in this section.

Not surprisingly, predicted Δ′17O values for O2 depend on both pO2 and pCO2. Thus, estimating pO2 from Proterozoic data depends on having an accurate estimate for pCO2. Unfortunately, constraints on pCO2 span an order of magnitude. A credible lower limit on pCO2 can be inferred from climate considerations, but the upper limit on pCO2 is unclear. The two dashed horizontal lines in Fig. 4 indicate a range of plausible pCO2 levels. The lower line is at ∼12 PAL (3,550 ppmv), and the upper line is at ∼24 PAL (7,100 ppmv). These estimates come from three-dimensional climate model simulations (73). According to their (zero-CH4) climate model results, 3,550 ppmv CO2 at 1.42 Ga would have produced a glacial climate with a mean surface temperature of 10.7 °C, whereas 7,100 ppmv CO2 would have produced a mean surface temperature of 18.6 °C—sufficient to keep the continents (but not the poles) ice free. These calculations account for the fact that solar luminosity at that time was ∼90% of the present value (74).

If one considers these values to be lower and upper bounds on pCO2 and assumes no dilution of the O-MIF signal, pO2 should be in the range of 0.8 to 20% PAL (Fig. 4). These values are roughly an order of magnitude higher than the 0.1 to 1% PAL values estimated by Planavsky et al. (37) from these same data. Part of the difference is attributable to the more detailed photochemical model used here and to our improved treatment of air–sea gas exchange. However, much of it stems from our more restricted range of plausible pCO2 values. Planavsky et al. (37) assumed a range of pCO2 values randomly distributed between 2 and 500 PAL, along with a substantial terrestrial oxygen flux. If we allow pCO2 values as low as 2 PAL, our predicted pO2 could be 0.1% PAL or even lower. These low CO2 levels would have resulted in global glaciation unless Earth’s greenhouse effect at that time was supplemented by large amounts of CH4 and/or N2O (73, 75). However, N2O photolyzes rapidly below ∼0.1 PAL O2, so its greenhouse contribution should have been modest (63, 76). CH4 is a more likely Proterozoic greenhouse gas; its possible effect is considered in the following section.

A weakness in the above argument is that Liu et al. (73) considered their 24 PAL pCO2 estimate to be a lower limit for mid-Proterozoic pCO2 (assuming no methane), not an upper limit. The climate record suggests that the “Boring Billion” period between 1.8 to 0.8 Ga was ice free, except for a few possible exceptions. Glaciations have been reported in the King Leopold formation in northwestern Australia at 1.8 Ga (77) and in the Vazante Group in Brazil, which was once considered Neoproterozoic in age but which has been redated at 1.1 to 1.3 Ga (78). Atmospheric CO2 levels above ∼30 PAL would permit pO2 values as high as today, according to Fig. 4. Thus, other arguments are needed to place an upper limit on pCO2 if we wish to constrain mid-Proterozoic pO2 levels.

Additional Results for Lower pO2 Levels (Terrestrial Biosphere Included).

Finally, we repeated the calculations shown in Fig. 4 but with a terrestrial (microbial mat) biosphere included. Terrestrial GPP was set at 660 Tmol (O2) ⋅ yr−1 with an accompanying CH4 flux equal to half that value (see Recycling by the Terrestrial Biosphere). We assumed O2 exchange rates with the marine biosphere were the same as in Results for Lower pO . The effect on Δ′17O of O2 is shown in Fig. 4. The red contours at −6 and −11‰ shift upwards, particularly at low pO2 levels where air–sea O2 exchange is slow. Most of this change is caused by higher rates of O2 recycling, which dilutes the O-MIF signal.

However, the interpretation of these results is changed even more if one considers the effect of higher CH4 levels on climate. The high CH4 flux assumed here, 10 times the modern value, causes atmospheric CH4 to increase from 1.6 ppmv in the base, 1 PAL O2 model to as much as 60 ppmv at some O2 levels (). The nonlinear response of CH4 concentration to its surface input rate has been seen in other models [e.g., Pavlov et al. (79)]. Even at low pO2, the calculated CH4 mixing ratio is 15 to 20 ppmv, enough to generate 2 to 3° C of greenhouse warming (76). One CO2 doubling produces roughly this same temperature increase when the climate is relatively warm; consequently, we reduced the lower pCO2 limit in Fig. 4 from 12 to 6 PAL. When this change is made, the permissible range for pO2 extends all the way down to 0.1% PAL or lower. If the CO2 level and CH4 flux assumed here are correct for the mid-Proterozoic, few constraints can be placed on pO2 from the existing O-MIF data. Terrestrial GPP values of 4,000 Tmol ⋅ yr−1 or higher (26) can be ruled out, though, because they produce Δ′17O values in O2 less negative than −6‰ ().

Discussion

Geologic Constraints on Mid-Proterozoic pCO2.

Past atmospheric CO2 concentrations can be inferred from geologic data as well as from climate models (54). One approach is to study the composition of ancient soils preserved as paleosols. The two samples that are closest in time to the period of interest are the Flin Flon paleosol at 1.8 Ga (80) and the Sturgeon Falls paleosol at 1.1 Ga (81). The original authors used the fact that these paleosols retained iron during weathering to derive estimates for pO2. Sheldon (82) reanalyzed them for pCO2 using a similar kinetic approach. Interpolating between his (very high) pCO2 value at 1.8 Ga and his much lower value at 1.1 Ga yielded a pCO2 of 3,600 ppmv (12 PAL in our units) at 1.4 Ga (82). This value should be considered as a lower limit, as Sheldon’s analysis method assumes that all the CO2 that enters the soil reacts with silicate minerals. A reanalysis of earlier Archean paleosols by Kanzaki and Murakami (83) using a different methodology yields pCO2 values that are 10 to 20 times higher than those calculated by Sheldon for those same samples. Thus, the paleosol data support the lower limit on pCO2 shown in Fig. 4 but may not provide a useful upper limit.

Our lower limit on mid-Proterozoic pCO2 is also supported by a study of the carbon isotopic composition of microfossils preserved in the Ruyang Group in Shanxi Province, China (84). Carbon isotope fractionation decreases when CO2 is less available; hence, the highly depleted δ13C values (−32 to −36‰) measured in the microfossils require relatively high pCO2. They estimate a lower limit of 3,600 ppmv (12 PAL) at 1.4 Ga, similar to Sheldon’s estimate from paleosols.

Upper limits on mid-Proterozoic pCO2 have been estimated by looking at fossil cyanobacteria as well as their effect on the nature of the sedimentary record. Riding (85) argued that a decline in pCO2 below 10,000 ppmv (33 PAL) at ∼1,400 to 1,300 Ma resulted in blooms of planktic cyanobacteria that induced “whitings” of carbonate mud in the water column whose sedimentary accumulation began to dominate carbonate platforms at that time. A later study of cyanobacterial sheath calcification suggests that 3,600 ppmv (12 PAL) is an upper limit on pCO2 at 1.2 Ga (86). These authors argue that falling CO2 levels led to the evolution of intracellular carbon-concentrating mechanisms that pumped CO2 into cell carboxysomes, raising the pH of the surrounding cytosol and triggering in vivo carbonate deposition.

We conclude that atmospheric pCO2 was likely falling during the mid-Proterozoic and that concentrations of 12 to 24 PAL at 1.4 Ga are in the right ballpark based on geologic proxy data. If so, then the very low pO2 levels allowed by Fig. 4 may not be supported. Moving the lower limit on pCO2 back to 12 PAL in Fig. 4 would imply a lower limit of ∼0.5% PAL for pO2. However, even this limit would not be firm because of possible dilution of the O-MIF signal during sulfate precipitation in the water column.

Additional Thoughts on Mid-Proterozoic CH4.

Our results depend critically on the concentration of atmospheric CH4. When CH4 is low, that is, when the biogenic methane flux is comparable to (or lower than) today, we can derive a lower limit on mid-Proterozoic pO2, ∼0.8% PAL, from the Sibley data if we assume no recycling of the O-MIF signal in the water column (Results for Lower pO). At 10 times the current methane flux, we get no lower limit on pO2 whatsoever (Additional Results for Lower pO Levels (Terrestrial Biosphere Included)). In the latter calculations, we assumed a substantial methane flux from terrestrial microbial mats. This assumption is speculative for two reasons: 1) The fractional land coverage and productivity of such mats during the mid-Proterozoic is uncertain. 2) Whether mats grown under low pO2 could actually generate this much methane is unclear. The cultivation of mats under such conditions might help to answer the latter question.

Some authors (e.g., ref. 79) have also suggested that methane could have been generated in significant amounts in Proterozoic marine sediments. Their argument was that, if GPP was relatively high and the recycling of organic matter by aerobic decay and sulfate reduction were both low, more of the organic matter should have been recycled by fermentation and methanogenesis. This argument has been criticized by other researchers (87–89) who argue that methane production in sediments remains low even under low-O2, low-sulfate conditions, partly because of the fact that biological productivity must have been lower because of slower recycling of critical nutrients such as phosphorus (88, 89). Unless this objection can be countered, we accept that marine methane production must have been low.

Conclusions

We provide a more mechanistically grounded framework for interpreting the triple oxygen isotope record through Earth’s history. Our model improves on the treatment of atmospheric ozone photochemistry and introduces key constraints on air–sea gas exchange. The O-MIF record provides a measure of the atmospheric ratio of O2:CO2, not of pO2 itself. A reasonable lower limit on pCO2 at 1.4 Ga, ∼3,600 ppmv or 12 PAL, can be estimated from climate models and is supported by geochemical proxy data. CO2 partial pressures lower than this value can be tolerated climatically but only if CH4 concentrations were high (>15 ppmv). Given pCO2 ≥ 12 PAL and low CH4, mid-Proterozoic O2 should have been >0.8% PAL if the Sibley data are taken at face value. The dilution of the Δ′17O signal by the reprocessing of sulfate within the water column above where the sulfates were deposited could decrease this lower limit on pO2, but there is no easy way to quantify this effect. High CH4 and lower pCO2 would also allow pO2 to have been lower.

Upper limits on pO2 during the mid-Proterozoic are difficult to estimate from the O-MIF data, again because of uncertainties in pCO2. Marine GPP cannot be estimated reliably from sedimentary O-MIF data because of complications imposed by constraints on air–sea gas transfer. Terrestrial microbial mat GPP is constrained by similar considerations. If atmospheric pO2 was low, then terrestrial GPP could only have been high if mats emitted large quantities of CH4 or H2/CO to the atmosphere to soak up the emitted O2. Despite all these limitations, the sedimentary O-MIF record remains a useful tool for investigating the composition of the Precambrian atmosphere. If better proxies for pCO2 can be obtained, along with more direct measures of Δ′17O of atmospheric O2, the O-MIF record may eventually lead to a quantitative understanding of the history of atmospheric O2.

Methods

Our photochemical model includes a main model for normal chemical reactions and an isotopic model designed for oxygen isotope reactions. The main model is derived from Stanton et al. (63). It extends from the ground up to 100 km in altitude and contains full atmospheric chemistry for CHONS (carbon, hydrogen, oxygen, nitrogen, and sulfur) species (up through C1 for carbon). Chlorine chemistry has been removed from the model, as it is less important for the preindustrial atmosphere than for today. The main model includes 42 species linked by 175 chemical reactions (). Boundary conditions for the long-lived species are described later in this section. The main model is run to steady state, assuming a fixed solar zenith angle and using a fully implicit time integration scheme.

The boundary conditions in the model are critical and must be chosen carefully. At the upper boundary (100 km), most species are given an effusion velocity of zero (equivalent to zero flux). Major species (O2 and CO2) that are rapidly photodissociated above this level are allowed to flow up, whereas their photolysis products (CO and O) flow down in stoichiometric ratios. The lower boundary condition for each long-lived species can be either fixed mixing ratio, fixed flux, or fixed deposition velocity. O2 and CO2 are always assigned fixed mixing ratios in the main model. They serve as control variables in our calculations. Most other long-lived species are assigned constant deposition velocities following Kasting et al. (90). In the modern atmosphere simulation, some biogenic or partly biogenic trace gases (CH4, N2O, H2, and CO) are given mixing ratios equal to their observed values (). The main model then calculates surface fluxes needed to sustain them. At lower pO2 levels, these fluxes are held constant in both the main and isotopic models. For CH4, the modern mixing ratio is fixed at 1.6 ppmv; this yields a corresponding upward flux of 33 Tmol ⋅ yr−1, which is comparable to estimates for the CH4 flux on modern Earth (91). In the calculations that include terrestrial microbial mats (Additional Results for Lower pO), the CH4 flux was increased to 330 Tmol ⋅ yr−1, balancing the assumed mat O2 flux.

The key to calculating the dilution of the O-MIF signal by biospheric recycling is to give special treatment to O2 and CO2. Unlike shorter-lived gases such as CH4, their surface mixing ratios do not depend in a simple way on their surface fluxes. The surface fluxes calculated by the photochemical model, given fixed mixing ratio boundary conditions, are tiny compared to the recycling fluxes described in Dilution of the O-MIF Signal in O. Moreover, in the (globally averaged) photochemical model, each gas must be flowing either up or down at the lower boundary, whereas, on the real Earth, O2 and CO2 are flowing both upward and downward in different places and at different times of the year. We simulate this exchange by inserting fictitious tropospheric chemical production terms for O2 and CO2 equal to their calculated recycling fluxes () in both main and isotopic codes. These sources are distributed evenly between 0 and 11 km, as should be appropriate for a long-lived gas injected into a rapidly mixed reservoir. The injected O2 and CO2 have nowhere to go except down through the lower boundary. Dividing those downward fluxes by the number density of each gas yields a deposition velocity, which is then saved for use in the isotopic code. For pO2 ≤ 1% PAL, to ensure that we maximize the rate of O2 exchange, we gradually increase the fictitious production term until the calculated O2 deposition velocity is equal to the PV described in Recycling by the Marine Biosphere (also see ).

The structure of the isotopic model is like that of the main model but with significant deviations in the chemical scheme and boundary conditions. The isotopic model contains 358 reactions that duplicate the chemistry in the main model but for isotopically substituted species (). Both minor O isotopes, 18O and 17O, are represented in the isotopic model by a single parameter, Q. Q is assumed to be sufficiently scarce so that species containing a Q atom never react with each other; hence, one never produces doubly substituted isotopic species. This is a good assumption for both 18O (∼0.2% of total oxygen) and 17O (∼0.04%). The isotope routine is run twice for each calculation, once for 17O and once for 18O. The δxO values found in this manner are combined to calculate Δ′17O by Eq. . Species in the isotopic model are initialized with values proportional to their concentrations in the main code. Those species containing one O atom are set equal to their main model counterpart, species containing two O atoms are set equal to twice their main model values, and so on. This ensures that each isotopically substituted species starts out with zero fractionation (because the isotopic composition of the species is related to the number of Q atoms [one] divided by the total number of O atoms in the molecule).

The boundary conditions for the species in the isotopic model are the same as in the main model except for O2 and CO2 (). Instead of assuming fixed mixing ratios, we use fixed deposition velocities for OQ and COQ. Each gas is given the same deposition velocity as that derived in the main code, along with the same distributed tropospheric recycling flux. Thus, their surface concentrations can now change in response to fractionating chemistry that occurs within the atmosphere. Because they start close to the solution, the isotopic equations can be iterated to steady state using Newton’s method. shows the boundary conditions used in the main/isotope codes and a flowchart of the model operation described in this section. Further details on the methodology are provided in .