Rapid timescale for an oxic transition during the Great Oxidation Event and the instability of low atmospheric O2.

The Great Oxidation Event (GOE), arguably the most important event to occur on Earth since the origin of life, marks the time when an oxygen-rich atmosphere first appeared. However, it is not known whether the change was abrupt and permanent or fitful and drawn out over tens or hundreds of millions of years. Here, we developed a one-dimensional time-dependent photochemical model to resolve time-dependent behavior of the chemically unstable transitional atmosphere as it responded to changes in biogenic forcing. When forced with step-wise changes in biogenic fluxes, transitions between anoxic and oxic atmospheres take between only 102 and 105 y. Results also suggest that O2 between [Formula: see text] and [Formula: see text] mixing ratio is unstable to plausible atmospheric perturbations. For example, when atmospheres with these O2 concentrations experience fractional variations in the surface CH4 flux comparable to those caused by modern Milankovich cycling, oxygen fluctuates between anoxic ([Formula: see text]) and oxic ([Formula: see text]) mixing ratios. Overall, our simulations are consistent with possible geologic evidence of unstable atmospheric O2, after initial oxygenation, which could occasionally collapse from changes in biospheric or volcanic fluxes. Additionally, modeling favors mid-Proterozoic O2 exceeding [Formula: see text] to [Formula: see text] mixing ratio; otherwise, O2 would periodically fall below [Formula: see text] mixing ratio, which would be inconsistent with post-GOE absence of sulfur isotope mass-independent fractionation.

Abundant atmospheric O2 at 21% by volume is the most distinctive and consequential feature of Earth’s atmosphere. Produced by cyanobacteria, algae, and plants, O2 is a clear sign of our biosphere that is detectable across interstellar space by telescopic spectroscopy (1). Oxygen permits aerobic respiration, the only known metabolism with sufficient energy yield to sustain complex animal life (2). However, for about the first half of Earth’s 4.5-billion-year-old history, the atmosphere had negligible O2 (e.g., ref. 3). This changed ~2.4 billion years ago.

The timing of the Great Oxidation Event (GOE) and the magnitude of atmospheric O2 concentrations before and after the GOE can be constrained by the geologic record of sulfur isotopes in combination with photochemical models. Archean and earliest Proterozoic sedimentary minerals contain sulfur isotopes with characteristic mass-independent fractionation (MIF) which abruptly disappears 2.4 billion years ago (4). Sulfur MIF in marine sediments likely requires that atmospheric photochemistry produce elemental sulfur, S8 (for explanation, see the introduction in ref. 5) (6, 7). Zahnle et al. (5) used a one-dimensional (1D) photochemical model to show that atmospheric S8 production only occurs when atmospheric O2 is below mixing ratio. An often cited threshold of was from an earlier photochemical model that did not simulate atmospheres with surface O2 mixing ratios between and (6). Therefore, the disappearance of the sulfur isotope MIF signal at 2.4 Ga is strong evidence that O2 first rose above mixing ratio then.

Geologic evidence may suggest that the GOE was not a single monotonic rise of oxygen but characterized by oscillations. Using U-Pb dating, Gumsley et al. (8) updated the chronology of sulfur isotope MIF in the stratigraphic record, finding evidence for two oxic-to-anoxic transitions between and Ga. More recently, Poulton et al. (9) report 2.3 Ga to 2.2 Ga marine sediments with sulfur isotopes consistent with approximately five oxic-to-anoxic transitions. Fluctuating O2 levels coincide with three to four widespread glaciations, indicating extreme climate instability (10). Overall, geochemical evidence tentatively suggests that O2 concentrations and climate were unstable for 200 million years until 2.2 Ga, which marks the most recent estimated timing of the permanent oxygenation of the atmosphere (9). However, interpretations of oscillating O2 have been questioned (11). While the geologic evidence for the O2 oscillations remains equivocal, the data have raised significant questions regarding the feasibility and timescales for Earth’s great oxidation. Some have argued that the oxygen-rich atmosphere is more stable than an oxygen-poor atmosphere (12), which favors a single rise of O2 instead of O2 oscillations.

Evidence for O2 instability and the time-dependent behavior of O2 concentrations has not been reconciled with atmospheric photochemical models. All previous models treated the GOE as successive photochemical steady states (5, 6, 13–19). A photochemical steady state occurs when no atmospheric species changes concentration over time, because their production and loss from reactions and surface sources (e.g., volcanoes or biology) are balanced. Such steady-state calculations have been crucial for understanding the GOE by contextualizing sulfur isotope MIF observations (5, 6), and establishing the relationship between atmospheric O2 concentrations and the degree to which O3 blocks UV photons from Earth’s surface (i.e., O3 shielding) (13, 15, 16), but they do not evaluate time-dependent changes and transient imbalances, or characteristic timescales.

Several theories for the rise of O2 suggest that it relied on a global redox titration over 108 y to 109 y involving oxidation of the upper mantle and/or crust, plausibly driven by hydrogen escape, which led to a tipping point where the source flux of O2 exceeded a kinetically rapid O2 sink from volcanic and metamorphic reductants (20–24). Beyond the tipping point, O2 flooded the atmosphere, reaching a new, long-term balance limited by oxidative weathering.

Here, we developed a time-dependent 1D photochemical model capable of investigating changes of O2 at the tipping point itself over timescales of 102 y to 105 y rather than the longer-term planetary changes which initiated the GOE. We simulate changing O2 as a time-dependent evolution, in contrast to the steady-state approach used in previous studies (e.g., ref. 13), because O2 can change on relatively rapid timescales that are not well characterized by steady states. With our model, we compute the time required for an anoxic-to-oxic atmospheric transition, and the time required for deoxygenation. Additionally, we investigate the stability of O2 concentrations against perturbations to surface gas fluxes produced by biology. Finally, we use our model results to better constrain O2 levels and stability during the GOE (starting at ~2.4 Ga), and during the mid-Proterozoic eon (1.8 Ga to 0.8 Ga).

Results

To investigate the time-dependent behavior of O2 during the GOE, we first computed grids of photochemical steady-state atmospheres. These steady states establish the context for time-dependent photochemical modeling described in subsequent sections. Fig. 1 shows the predicted steady-state surface O2 mixing ratio (Fig. 1), the surface CH4 mixing ratio (Fig. 1), and the precipitation of atmospheric S8 (gray shading) as a function of surface O2 flux between and 1013 molecules per cm2 ⋅s−1, and CH4 flux/O2 flux ratios between 0.27 and 0.49 where gas fluxes are those entering the atmosphere. The surface O2 fluxes reported here are net emissions into the atmosphere which exclude recycling within the biosphere. For reference, a comparable model of the modern Earth requires a surface O2 flux of 1012 molecules per cm2 ⋅s−1, and a CH4 flux/O2 flux of 0.09 (CH4 flux molecules per cm2 ⋅s−1) (16). We consider a CH4 flux/O2 flux ratio close to 0.5 to be more realistic for the Late Archean, prior to the GOE, because this ratio is expected if oxygenic photosynthesis is balanced by methanogenesis. In net (21), where “CH2O” represents organic matter,

Fig. 1.

Colored contours show photochemical steady states of (A) surface O2 mixing ratio and (B) surface CH4 mixing ratio as a function of O2 surface flux and CH4 flux/O2 flux. Gray shading indicates the magnitude of elemental S8 production in the atmosphere, which is considered essential for the preservation of sulfur isotope MIF in marine sediments. Peak S8 production is molecules per cm2 ⋅s−1. Gray shading fades to white for S8 production less than molecules per cm2 ⋅s−1, a negligibly small value. Arrows labeled “Figure 2a” and “Figure 2c” indicate start and end points for time-dependent photochemical models of the oxic transition shown in Fig. 2 . Red, blue, and green stars are the initial conditions used in the simulations shown in Fig. 4 , respectively.

Fig. 2.

Three models of anoxic-to-oxic transitions. (A) Atmospheric oxygenation caused by a step-wise decrease in the methane flux from to molecules per cm2 ⋅s−1 (orange arrow in Fig. 1). (i) Surface O2 and CH4 mixing ratios, and O3 and S8 column abundance over time; (ii) OH surface mixing ratio and tropospheric H2O photolysis rate. (B) Transition caused by step-wise increase in the O2 flux from 1012 to molecules per cm2 ⋅s−1 and a stepwise increase in the CH4 flux to maintain constant CH4 flux/O2 flux = 0.45 (). Transition in C results from a step-wise decrease in the CH4 flux from to molecules per cm2 ⋅s−1 (black arrow in Fig. 1).

Fig. 4.

The photochemical stability of O2. Shading in A shows the steady-state inertial timescale of redox gases (Eq. ), and colored contours are the steady-state surface O2 mixing ratio (same as Fig. 1). (B–D) Time-dependent photochemical simulations with oscillating CH4 surface fluxes, each beginning with steady-state atmospheres indicated in A. O2 stability is directly proportional to the column abundance of redox gases in the atmosphere.

The CH4 flux/O2 flux ratio is smaller than 0.5 on modern Earth largely because of the microbial anerobic oxidation of methane via in ocean sediments, a process that was unimportant in the anoxic mid-Archean ocean with little sulfate (25–27). We include O2 fluxes several orders of magnitude smaller than the modern value ( molecules per cm2 ⋅s−1) because of evidence for smaller primary productivity during the Proterozoic eon (28, 29).

Recall that atmospheric S8 deposition is considered necessary to preserve sulfur isotope MIF in ocean sediments (6). We find that S8 production is not possible above a O2 mixing ratio (the gray-to-white shading boundary in Fig. 1), consistent with previous results (5).

Fig. 1 uses Archean outgassing surface fluxes for CO, H2, H2S, and SO2 listed in Table 1, with the CO2 surface mixing ratio fixed to 1% for all model runs—a reasonable value for the Late Archean according to carbon cycle modeling (30). Additionally, over the same span of surface O2 fluxes and H2 flux/O2 flux ratios, we compute photochemical steady states for the modern fluxes for CO, CH4, H2S, and SO2 listed in Table 1, again fixing CO2 to 1%. The results are shown in .

Table 1.

Fixed surface flux boundary conditions for SO2, H2S, H2, and CO used in this study

Model	SO2	H2S	H2	CO
Archean outgassing*	1010	109	3×10^10	3×10^9
Modern values†	3.5×10^9	3.5×10^8	1.22×10^8	2.65×10^11

All fluxes have units of molecules per square centimeter per second.

*The same fluxes as the “Archean High” values from table 1 in Zahnle et al. (5).

†Surface flux values required to reproduce the concentration of each gas in modern Earth’s atmosphere. These values are also the “Case 1” fluxes described in Gregory et al. (16).

In the following sections, we calculate the time required to transition between different steady-state atmospheres shown in Fig. 1 and .

The Timescale of Atmospheric Oxygenation.

The orange arrow labeled “Figure 2a” in Fig. 1 corresponds to the approximate start and end states of the time-dependent photochemical model run shown in Fig. 2. The model starts with an atmosphere at a steady state, then, at y, we impose a stepwise decrease in the surface methane flux from to molecules per cm2 ⋅s−1 (we keep the surface O2 flux constant at 1012 molecules per cm2 ⋅s−1). This perturbation causes O2 to rise from to mixing ratio over 3,500 y, eliminating photochemical S8 production.

The O2 transition in Fig. 2 , i is modulated by O3 shielding of tropospheric H2O (13). When a stratospheric O3 layer begins to develop, OH production from H2O decreases (Fig. 2 , ii). Decreasing OH concentrations prevent the mutual annihilation of O2 and CH4 (by followed by ), so O2 levels increase. The mixing ratio of CH4 also rebounds. O3 shielding (protecting life on the surface from harmful solar UV radiation) is just barely beginning to operate in this example compared to modern Earth. After the atmosphere reaches a new steady state, the atmospheric column has O3 molecules per cm2, some 26 times smaller than the modern value of molecules per cm2. Note that the extent to which O3 shields tropospheric H2O can be strongly modulated by 3D dynamical effects (31), which we do not account for.

Like Fig. 2, Fig. 2 also shows a transition between a and O2 mixing ratio, but this model uses the modern outgassing fluxes for CO, H2, H2S, and SO2 listed in Table 1 instead of presumptive Archean outgassing values. Also, at t = 0 y, we impose a stepwise increase of the O2 flux from 1012 to molecules per cm2 ⋅s−1 while keeping the CH4 flux/O2 flux ratio at 0.45 (see for context). While the anoxic-to-oxic transition itself still occurs rapidly, the atmosphere simulated in Fig. 2 takes 60,000 y to reach the tipping point, which is much longer than the comparable O2 transition shown in Fig. 2.

The time required for O2 to begin to rise in concentration is controlled by the reservoir of reducing gases, primarily CH4, H2, and CO, in the preoxygenated atmosphere. Big reservoirs of reducing gases slow the timescale of oxygenation, because reducing gases must be mostly removed before O2 can increase. O2 cannot increase while reducing gases are abundant, because large oxygen sinks from reactions with reducing gases prevent it. That is why 3,500 y elapse before O2 begins to rise in Fig. 2, and why 60,000 y elapse before O2 rises in Fig. 2. Fig. 2 starts with more reducing gases, which take longer to eradicate.

We can roughly estimate the time required for O2 to begin to rise, with a back-of-the-envelope calculation of the rate at which reducing gases are eliminated from the anoxic atmosphere. The total reservoir of reducing gases in the preoxygenated atmosphere in O-equivalent units is

Here, is the O-equivalent column abundance of reducing gases (O molecules per cm2) which is equal to the sum of all reducing gases in the atmosphere (N) multiplied by α, the redox state of each gas. Redox state is a relative quantity that requires defining redox-neutral reference species. Following previous models of early Earth (5), we define H2O, SO2, CO2, and N2 as redox neutral, with the oxygen redox parameter . Therefore, , and , from redox stoichiometry of hydrogen, sulfur, and carbon, respectively. It then becomes straightforward to calculate the α for any molecule. For example, . For a more in-depth explanation of atmospheric redox, see section 3 in Harman et al. (32) or chapter 8 in Catling and Kasting (33). is approximately equal to the weighted sum of CH4, H2, and CO because these are the main reducing gases in an Archean Earth-like atmosphere.

The change in column abundance of reducing gases is the difference between the redox columns at the finial and initial atmospheric states.

In Fig. 2 , we initiate the rise of oxygen by changing the surface flux of CH4 and/or O2 flux. We can quantify this flux perturbation in units of O molecules per square centimeter per second (),

Therefore, the time required to oxidize the reducing gases in the atmosphere and permit oxygen to begin rising is approximately

Plugging in values for the O2 transition in Fig. 2 yields y, a value only slightly smaller than the 3,500 y predicted by the time-dependent photochemical model. For Fig. 2, y, which is about a factor of 2 smaller than the figure from 1D photochemistry. Our estimate is too small in this case because the reducing column and its destruction rate are not constant prior to the rise of oxygen (). This calculation illustrates that the time required for oxygen to begin rising, once a tipping point of fluxes is reached, mostly depends on the quantity of reducing gases in the preoxygenated atmosphere.

Fig. 2 shows a more substantial anoxic-to-oxic transition compared to simulations shown thus far (also see Movie S1). We start with the same steady-state atmosphere as in Fig. 2, except we decrease the methane flux by twice as much at t = 0, from molecules per cm2 ⋅s−1 to molecules per cm2 ⋅s−1 instead of molecules per cm2 ⋅s−1 (we keep the surface O2 flux constant at 1012 molecules per cm2 ⋅s−1). O2 begins to rise and eliminates S8 production after ~1,500 y, but O2 will reach higher levels because of the lower CH4 flux. It takes y for O2 to reach its final steady-state abundance of mixing ratio. While the switch from to O2 mixing ratio remains as rapid as in Fig. 2, the predicted increase in O2 concentrations to requires far longer. This timescale is roughly analogous to the time required to deplete H2 and CH4 reservoirs to allow O2 to initially rise in concentration.

In summary, the timescale for O2 to rise in concentration depends on the reservoir of redox gases in the atmosphere, and the magnitude of the perturbation to redox surface fluxes. For O2 to rise from to , reducing gases must first be removed, which can take thousands to 10 thousands of years (Fig. 2 ). Increasing O2 concentrations beyond to near percentage levels requires filling a large O2 reservoir, which occurs on 105-y timescales in our model run (Fig. 2).

The Time Required for Deoxygenation.

Here, we use our time-dependent photochemical model to address the controversy of the reversibility of the oxic transition (9, 11). Fig. 3 shows the reverse of model runs shown in Fig. 2. For each model run, we start with an atmosphere initially at a steady state at the end of the simulations shown in Fig. 2. Then we impose a stepwise change in the O2 and CH4 flux at t = 0 y to return the atmosphere to anoxia. The reversal of Fig. 2 takes ~700, 100, and 40,000 y, respectively, in comparison to the 3,500, 60,000, and 300,000 y required for oxygenation.

Fig. 3.

(A–C) Simulated reversal of the oxic transitions shown in Fig. 2 , B, and C, respectively. Each oxic-to-anoxic transition is caused by a stepwise change of the CH4 flux and O2 flux at t = 0 y.

Like the timescale for oxygenation, the timescale of deoxygenation depends on the column abundance of redox-sensitive gases. In the previous section, we established that the timescale required for O2 to begin to rise is merely the time required to deplete the reservoirs of CH4 and other reducing gases. Analogously, the timescale of deoxygenation is determined by the reservoir of O2 and other oxidizing gases in the oxygenated atmosphere. The reversal shown in Fig. 3 starts with only O2, which can be depleted very quickly, allowing the return of an anoxic atmosphere. In contrast, the reversal shown in Fig. 3 takes 40,000 y because the atmosphere starts with O2, which takes longer to deplete.

The Stability of Post-GOE Atmospheric Oxygen.

In the previous two sections, we show that reservoirs of redox gases, primarily methane and oxygen, give the atmosphere chemical inertia, controlling the timescale of O2 changes. When reservoirs are big, for similar flux perturbations, the O2 mixing ratio will change relatively slowly over time; however, when reservoirs are small, photochemistry permits rapid O2 transitions. Therefore, the abundance of redox gases in an atmosphere is closely linked to the photochemical stability of oxygen.

Fig. 4 shows the steady-state inertial timescale of redox gases, , over the same axes as Fig. 1, which shows mixing ratios. is the sum of all redox gases in the atmospheric column (O-equivalent molecules per square centimeter), divided by a characteristic flux perturbation, which we take to be 10% of the O2 flux,

We choose the characteristic flux perturbation to be 10% of the O2 flux because it is the same order of magnitude as natural redox variations that occur on modern Earth during Milankovitch cycling (see ). An upper limit for the characteristic flux perturbation would be 100% of the O2 flux. This would decrease all values in Fig. 4 by a factor of 10, which would not change our interpretation. Since CH4, CO, H2, and O2 are the most important redox gases, the numerator in Eq. is well approximated by . Oxygen is the most prone to change for the smallest values, coinciding with O2 mixing ratios between and shown in the whitish region of Fig. 4.

The time-dependent photochemical models shown in Fig. 4 illustrate the relationship between and O2 instability. To produce Fig. 4, we started with the steady-state atmosphere indicated on Fig. 4, then imposed 17% amplitude oscillations to the CH4 flux with a period of 10,000 y. This forcing had no perceptible effect on the atmospheric O2. A similar 20% CH4 flux oscillation also did not significantly perturb an oxic atmosphere starting with O2 (Fig. 4). However, just 5% CH4 flux oscillations cause approximately four-orders-of-magnitude oscillations in surface oxygen concentrations for an incipiently oxic atmosphere starting with O2 (Fig. 4). O2 is most unstable where the abundance of all redox gases is smallest relative to a characteristic redox surface flux (the whitish area of Fig. 4) between and O2 mixing ratio. Stability continually increases outside of this range of O2 concentrations.

While O2 was relatively stable in the Fig. 4 simulations, it does not mean these atmospheres and initial oxygen concentrations are stable to all atmospheric perturbations. The stability of any O2 mixing ratio depends on the atmospheric forcings that are likely in nature. In , we argue that the CH4 flux oscillations used in Fig. 4 are realistic because comparable fractional changes in the methane flux have occurred over the past 650,000 y.

Flux oscillations over timescales greater than y are required to significantly affect O2 concentrations. Imposing 100% amplitude fluctuations to the CH4 flux with a period of 1 y, starting with the same atmosphere as Fig. 4, did not significantly alter the atmosphere over time. Atmospheres with between and O2 contain some CH4 and O2, which gives the atmosphere inertia against annual to decadal change.

Discussion

Recently, Gregory et al. (16) computed photochemical steady-state atmospheres for a wide range of surface O2 and CH4 fluxes and found bistable O2 concentrations. Their model allows steady-state atmospheres for O2 concentrations below and above mixing ratio but admits few steady-state solutions in between. They hypothesized that feedbacks between O2 and O3 shielding eliminate most solutions with these intermediate O2 concentrations. In contrast, our model can yield a steady-state solution with intermediate O2 concentrations given the right constant surface flux boundary conditions [e.g., Fig. 2 see also Gregory et al.’s (16) figure 8]. The difference might be caused by different steady-state convergence criteria, chemical reaction networks, boundary conditions, or a combination of these factors.

However, a photochemical steady state, for example, at O2 mixing ratio, does not mean that such an atmosphere is stable and realistic over 105- to 106-y timescales. Gas fluxes from Earth’s surface can vary during these timescales and significantly change O2 concentrations (Results).

For example, in the past 650,000 y, the biogenic methane flux has oscillated with an amplitude of 25% ( molecules per cm2 ⋅s−1) and a 100,000-y period (34, 35). On modern Earth, methanogens in wetlands are a major source of atmospheric methane (36). Every 100,000 y, ice sheets have advanced and retreated, covering and uncovering wetlands, changing the CH4 flux to the atmosphere. These ice ages and methane flux variations are in response to Milankovich cycles with characteristic periods between 20,000 and 100,000 y. This exact same methane oscillation would not have occurred in the Late Archean or Early Proterozoic because modern wetlands did not exist then, but a similar process involving microbial mats is conceivable.

Zhao et al. (37) modeled cyanobacterial mats on Proterozoic land, finding that they could have been a substantial CH4 source to the atmosphere. Ice sheets covering and uncovering microbial mats could have affected global CH4 fluxes. Fig. 4 illustrates the effect of 5% methane flux variations over Milankovitch timescales on an atmosphere starting with O2. O2 oscillates nearly four orders of magnitude between (anoxic) and (oxic) (Fig. 4).

An oscillating methane flux is only one of many possible atmospheric perturbations. The Early Proterozoic geologic record preserves evidence of large igneous provinces (LIPs), or massive volcanic eruptions (8). In , we show that the H2 and CO outgassed from a significant LIP eruption could cause the O2 surface mixing ratio to drop from to in y, causing a return to sulfur isotope MIF. In this simulation, we use the maximum LIP eruption rates reported in the literature (38). In addition to LIPs, a Snowball Earth event concurrent with the GOE would have presumably affected gases produced by the biosphere (10).

Constant surface gas fluxes from biology and volcanism for millions of years in the aftermath of the initial rise of O2 are unlikely. Additionally, our photochemical modeling shows that, for atmospheres with transitional O2 concentrations, relatively small atmospheric perturbations (e.g., a CH4 flux change of 5%) over timescales as short as hundreds of years can cause O2 to change by orders of magnitude (e.g., Fig. 3). Therefore, substantial variability of O2 during the GOE appears possible.

For these reasons, our photochemical modeling results are compatible with recently published evidence of fluctuating sulfur isotope MIF (8, 9) indicating that O2 was unstable between 2.4 and 2.2 Ga. We find that shutoff of S8 aerosol production, which is required to produce sulfur isotope MIF, occurs at O2 mixing ratio, a region of the parameter space where O2 is prone to rapid change (Fig. 4). But, oxygen surface levels between and mixing ratio were likely unstable. Short period changes to the biosphere, or volcanic outgassing rates, could have caused order of magnitude O2 changes over 100- to 100,000-y timescales. Occasionally, big perturbations to the atmosphere, such as an LIP, might have lowered O2 concentrations enough for sulfur isotope MIF to reoccur. Note that the above explanation for the cause of O2 oscillations prior to 2.2 Ga is complicated by S-MIF data presented in Izon et al. (11), which do not suggest the same O2 variability found by Poulton et al. (9).

After 2.2 Ga, and during the mid-Proterozoic, sulfur isotope MIF never returned. Therefore, this time must have had O2 concentrations large enough to prevent O2 collapse. Our modeling shows that larger O2 concentrations give the atmosphere chemical inertia, slowing atmospheric deoxygenation (Fig. 3). It is therefore challenging to reconcile our modeling results with the interpretation of Planavsky et al. (39), who used Proterozoic chromium isotopes to argue that O2 could not have been larger than mixing ratio. Such a small O2 reservoir would have been unstable to LIP eruptions, or variations in the CH4 flux from Milankovitch cycles (Fig. 4), which both have evidence of occurring in the mid-Proterozoic stratigraphic record (40–42). We conclude that, for stability, mid-Proterozoic O2 levels should have exceeded . This conclusion is compatible with mid-Proterozoic Fe isotopes in ironstones, which suggest O2 levels between approximately and mixing ratio (43).

Our results are not sensitive to the changing solar UV photon flux between the GOE ( Ga) and the mid-Proterozoic. Recalculating Fig. 1 using the solar UV flux at 1.3 Ga (44) results in surface O2 and CH4 surface mixing ratios within a factor of 2 of Fig. 1.

Our work also has implications for the most likely oxygen levels before the GOE. Johnson et al. (45) analyzed molybdenum isotopes in the Archean sedimentary record for signs of continental oxidative weathering. Their work is compatible with two end-member interpretations: 1) If Archean O2 was evenly distributed over the globe, then the surface O2 mixing ratio was and , or 2) if O2 accumulation was geographically restricted, then the O2 surface flux was greater than 0.01 Tmol ⋅y−1 ( molecules per cm2 ⋅s−1). Our modeling suggests interpretation 1 is unlikely because we find that O2 is likely unstable over geologic time for this range of oxygen levels.

In our modeling, we do not explicitly consider redox reservoirs in the oceans, sediments, crust, and mantle, for good reason. These reservoirs are coupled to the atmosphere and can modulate O2 levels. However, the timescale of equilibration between the atmosphere and other redox reservoirs is often relatively long (e.g., 100 My for organic carbon in continental sediments), so we consider them to be approximately constant over the timescale of O2 transitions ( y).

A caveat is that the coupling between redox reservoirs in the atmosphere, crust, or sediments might depend on atmospheric composition. An example is the pyrite oxidation rate, which depends on the partial pressure of oxygen (16). We do not explicitly consider such feedbacks, which could affect the timescale of changing O2 levels.

An additional, related caveat is that our model does not consider biological feedbacks. The rise of O2 would limit habitats for anaerobes, and permit more widespread aerobic respiration, potentially dropping the CH4 flux/O2 flux ratio farther than we have modeled here. Also, a stronger ozone UV shield would make new habitats for cyanobacteria and allow the expansion of life on land, promoting chemical and oxidative weathering. All these changes, which we do not explicitly model, would modulate oxygen levels. In this article, we impose changes in the CH4 and O2 flux that are supposed to be representative of a changing biosphere, but a better model would determine more realistic changes in gas fluxes by directly coupling 1D photochemistry and biology.

Conclusions

Our time-dependent photochemical modeling of the GOE suggests that oxygen can rise and fall over geologically short periods of time. For an anoxic-to-oxic transition, once a tipping point of imbalanced redox fluxes is reached, the reservoir of reducing gases in the atmosphere must be eliminated before O2 can begin to rise. This takes hundreds to 10 thousands of years. O2 accumulation to just hundredths or tenths of percent levels requires filling a large O2 reservoir, which may occur on a 105-y timescale. Atmospheric deoxygenation occurs over similar periods of time, mainly controlled by the magnitude of the initial O2 abundance.

We also find O2 instability, especially for mixing ratios between and . For these O2 concentrations, photochemistry demands that both CH4 and O2 be relatively small in concentration. This small reservoir of redox-sensitive gases permits rapid changes to the atmosphere’s redox state. For example, for an atmosphere starting with O2, 5% amplitude oscillations to the methane flux with a period of 10,000 y cause oxygen to fluctuate four orders of magnitude between anoxic and oxic. Additionally, we show that a LIP eruption could cause the collapse of O2 and the return of sulfur isotope MIF for an atmosphere starting with O2 mixing ratio ().

We emphasize that the short-term (102 to 105 y) variability in O2 levels considered here occurred on the backdrop of the billion-year oxidation of the crust and mantle, and long-term organic burial, which are argued to be the ultimate causes of the rise of oxygen on Earth (e.g., ref. 20).

Overall, our modeling is compatible with, but does not prove, proposed geologic evidence for fluctuating and unstable atmospheric O2 after the initial rise of oxygen 2.4 billion years ago. A single, unidirectional, oxidation event remains plausible, although it would require strong and perhaps biological feedbacks promoting permanent substantial changes in the global CH4 flux/O2 flux ratio. While this is evident between the Archean (CH4 flux/O2 flux 0.5) and modern (CH4 flux/O2 flux 0.1) biospheres, the dynamics of the Proterozoic biosphere remain largely unexplored. Also, our results suggest that a stable, post-GOE, mid-Proterozoic atmosphere would need an O2 mixing ratio exceeding a value in the to range.

Materials and Methods

To investigate the transition between an anoxic and oxygen-rich Earth, we use a photochemical model with one spatial dimension of altitude, approximating a global average vertical profile. One-dimensional photochemical models are typically governed by a simplification of the continuity equation for molecules,

Table 2 defines all the variables and their units. Here, the flux () is given by

Table 2.

Variables in Eq. 7

Variable	Definition	Units
f i	Mixing ratio of species i	Dimensionless
ni	Number density of species i	Molecules per cubic centimeter
z	Altitude	Centimeters
t	Time	Seconds
n	Total number density	Molecules per cubic centimeter
Pi	Total chemical production of species i	Molecules per cubic centimeter per second
Li	Total chemical loss of species i	Molecules per cubic centimeter per second
R_i, rainout	Production and loss of species i from rainout	Molecules per cubic centimeter per second
Q_i, lightning	Production and loss of species i from lightning	Molecules per cubic centimeter per second
Φ_i	Vertical flux of species i	Molecules per square centimeter per second
K	Eddy diffusion coefficient	Square centimeters per second
Di	Molecular diffusion coefficient	Square centimeters per second
Hi	=N_akT/μ_ig, The scale heights of species i	Centimeters
Ha	=N_akT/μ¯g, The average scale height.	Centimeters
Na	Avogadro’s number	Molecules per mole
k	Boltzmann’s constant	Ergs per kelvin
μ	Molar mass. μ¯ is mean molar mass of the atmosphere, and μi is the molar mass of species i	Grams per mole
g	Gravitational acceleration	Centimeters per square second
α_Ti	Thermal diffusion coefficient of species i. We neglect this term (α_Ti=0)	Dimensionless
T	Temperature	Kelvins

The above system of partial differential equations (PDEs) describes how the number density (n) of each chemical species i changes over altitude and time.

In our photochemical model, we solve a simplified version of Eq. which assumes that the total number density does not change over time (). This assumption is valid for atmospheric transitions which maintain approximately constant surface pressure and atmospheric temperature. The continuity equations can then be written in terms of mixing ratios (f) instead of number densities (see Appendix B.1 in ref. 33 for a derivation),

To approximate Eq. , the model replaces the spatial derivatives with finite difference approximations, turning the system of PDEs into a larger system of ordinary differential equations (ODEs). This is the “method of lines” approach to solving a PDE. Catling and Kasting (33), their Appendix B.2, provides a detailed description of how to do this with Eq. ; therefore, we will omit a detailed description here, except to point out a sign error. The first two terms for the equation for B in their equation B.16 should have minus signs instead of plus signs.

The system of ODEs derived from finite differencing Eq. can be evolved forward in time with numerical integration. However, the photochemical ODEs are “stiff,” meaning that some dependent variables (i.e., the mixing ratios) change much more quickly than others. For example, in the modern atmosphere, OH typically has a chemical lifetime of about 1 s, while CH4 has a chemical lifetime of y. Stiff equations require special, high-stability, “implicit” integration methods. For more details on stiff equations and the implicit methods used to solve them, see ref. 46.

Often, we solve for steady states of the photochemical continuity equation (). To find steady states, we begin with some initial atmospheric composition, then integrate Eq. forward in time until the atmosphere ceases to change, that is, a steady state is reached. The assumption of photochemical steady state is approximately valid for most periods of Earth’s history, because the atmosphere changes slowly enough to be in a quasi-steady state.

However, the Paleoproterozoic rise of O2 was a relatively fast atmospheric transition that is not well modeled as a photochemical steady-state process. Therefore, describing it requires accurately solving the continuity equation over time.

To model the photochemistry of the GOE, we modified the photochemical model contained within the Atmos modeling suite (described in Appendix B of ref. 33) so that it can accurately solve the time-dependent behavior of Eq. . We call the modified version of the model PhotochemPy. Instead of using a Backward Euler as in Atmos, we used CVODE BDF ODE solver from Sundials Computing (47). CVODE BDF is an implementation of the backward differential formulas (BDF) and is a gold standard for solving large chemical kinetics problems. For details, see .

PhotochemPy is open source under a Massachusetts Institute of Technology license. The version of the code (v0.2.14) used in this paper and the corresponding Python scripts to reproduce work done in this article are at https://zenodo.org/record/6824092. However, the most up-to-date version of the code can be found at the following GitHub link: https://github.com/Nicholaswogan/PhotochemPy.