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

Significance Understanding the rise of atmospheric oxygen on Earth is important for assessing precursors to complex life and for evaluating potential future detections of oxygen on exoplanets as signs of extraterrestrial biospheres. However, it is unclear whether Earth’s initial rise of O2 was monotonic or oscillatory, and geologic evidence poorly constrains O2 afterward, during the mid-Proterozoic (1.8 billion to 0.8 billion years ago). Here, we used a time-dependent photochemical model to simulate oxygen’s rise and the stability of subsequent O2 levels to perturbations in supply and loss. Results show that large oxygen fluctuations are possible during the initial rise of O2 and that Mesoproterozoic O2 had to exceed 0.01% volume concentration for atmospheric stability.

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 onedimensional 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 10 2 and 10 5 y. Results also suggest that O 2 between ∼10 −8 and ∼10 −4 mixing ratio is unstable to plausible atmospheric perturbations. For example, when atmospheres with these O 2 concentrations experience fractional variations in the surface CH 4 flux comparable to those caused by modern Milankovich cycling, oxygen fluctuates between anoxic (∼10 −8 ) and oxic (∼10 −4 ) mixing ratios. Overall, our simulations are consistent with possible geologic evidence of unstable atmospheric O 2 , after initial oxygenation, which could occasionally collapse from changes in biospheric or volcanic fluxes. Additionally, modeling favors mid-Proterozoic O 2 exceeding 10 −4 to 10 −3 mixing ratio; otherwise, O 2 would periodically fall below 10 −7 mixing ratio, which would be inconsistent with post-GOE absence of sulfur isotope mass-independent fractionation.
Great Oxidation Event | photochemistry | oxygen Abundant atmospheric O 2 at 21% by volume is the most distinctive and consequential feature of Earth's atmosphere. Produced by cyanobacteria, algae, and plants, O 2 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 O 2 (e.g., ref. 3). This changed ∼2.4 billion years ago.
The timing of the Great Oxidation Event (GOE) and the magnitude of atmospheric O 2 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, S 8 (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 S 8 production only occurs when atmospheric O 2 is below ∼ 2 × 10 −7 mixing ratio. An often cited threshold of 2 × 10 −6 was from an earlier photochemical model that did not simulate atmospheres with surface O 2 mixing ratios between 2 × 10 −6 and ∼10 −15 (6). Therefore, the disappearance of the sulfur isotope MIF signal at 2.4 Ga is strong evidence that O 2 first rose above 2 × 10 −7 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 ∼2. 4  atmosphere is more stable than an oxygen-poor atmosphere (12), which favors a single rise of O 2 instead of O 2 oscillations.
Evidence for O 2 instability and the time-dependent behavior of O 2 concentrations has not been reconciled with atmospheric photochemical models. All previous models treated the GOE as successive photochemical steady states (5,6,(13)(14)(15)(16)(17)(18)(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 O 2 concentrations and the degree to which O 3 blocks UV photons from Earth's surface (i.e., O 3 shielding) (13,15,16), but they do not evaluate time-dependent changes and transient imbalances, or characteristic timescales.
Several theories for the rise of O 2 suggest that it relied on a global redox titration over 10 8 y to 10 9 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 O 2 exceeded a kinetically rapid O 2 sink from volcanic and metamorphic reductants (20)(21)(22)(23)(24). Beyond the tipping point, O 2 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 O 2 at the tipping point itself over timescales of 10 2 y to 10 5 y rather than the longer-term planetary changes which initiated the GOE. We simulate changing O 2 as a time-dependent evolution, in contrast to the steady-state approach used in previous studies (e.g., ref. 13), because O 2 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 O 2 concentrations against perturbations to surface gas fluxes produced by biology. Finally, we use our model results to better constrain O 2 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 O 2 during the GOE, we first computed grids of photochemical steady-state atmospheres. These steady states establish the context for timedependent photochemical modeling described in subsequent sections. Fig. 1 shows the predicted steady-state surface O 2 mixing ratio (Fig. 1A), the surface CH 4 mixing ratio (Fig. 1B), and the precipitation of atmospheric S 8 (gray shading) as a function of surface O 2 flux between 3 × 10 9 and 10 13 molecules per cm 2 · s −1 , and CH 4 flux/O 2 flux ratios between 0.27 and 0.49 where gas fluxes are those entering the atmosphere. The surface O 2 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 O 2 flux of 10 12 molecules per cm 2 · s −1 , and a CH 4 flux/O 2 flux of 0.09 (CH 4 flux = ∼9 × 10 10 molecules per cm 2 · s −1 ) (16). We consider a CH 4 flux/O 2 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 "CH 2 O" represents organic matter, The CH 4 flux/O 2 flux ratio is smaller than 0.5 on modern Earth largely because of the microbial anerobic oxidation of methane via SO 2− 4 in ocean sediments, a process that was unimportant in the anoxic mid-Archean ocean with little sulfate (25)(26)(27). We include O 2 fluxes several orders of magnitude smaller than the modern value (∼10 12 molecules per cm 2 · s −1 ) because of evidence for smaller primary productivity during the Proterozoic eon (28,29).
Recall that atmospheric S 8 deposition is considered necessary to preserve sulfur isotope MIF in ocean sediments (6). We find that S 8 production is not possible above a ∼10 −7 O 2 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, H 2 , H 2 S, and SO 2 listed in Table 1, with the CO 2 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 O 2 fluxes and H 2 flux/O 2 flux ratios, we compute photochemical steady states for the modern fluxes for CO, CH 4 , H 2 S, and SO 2 listed in Table 1, again fixing CO 2 to 1%. The results are shown in SI Appendix, Fig. S1.  3 × 10 10 3 × 10 9 outgassing * Modern 3.5 × 10 9 3.5 × 10 8  In the following sections, we calculate the time required to transition between different steady-state atmospheres shown in  Fig. 2A. The model starts with an atmosphere at a steady state, then, at t = 0 y, we impose a stepwise decrease in the surface methane flux from 4.9 × 10 11 to 4.7 × 10 11 molecules per cm 2 · s −1 (we keep the surface O 2 flux constant at 10 12 molecules per cm 2 · s −1 ). This perturbation causes O 2 to rise from 3 × 10 −8 to 3 × 10 −5 mixing ratio over 3,500 y, eliminating photochemical S 8 production.

A B
The O 2 transition in Fig. 2 A, i is modulated by O 3 shielding of tropospheric H 2 O (13). When a stratospheric O 3 layer begins to develop, OH production from H 2 O decreases ( Fig. 2 A, ii). Decreasing OH concentrations prevent the mutual annihilation of O 2 and CH 4 (by CH 4 + OH → CH 3 + H 2 O followed by CH 3 + O 2 → products), so O 2 levels increase. The mixing ratio of CH 4 also rebounds. O 3 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 3 × 10 17 O 3 molecules per cm 2 , some 26 times smaller than the modern value of 8 × 10 18 molecules per cm 2 . Note that the extent to which O 3 shields tropospheric H 2 O can be strongly modulated by 3D dynamical effects (31), which we do not account for.
Like Fig. 2A, Fig. 2B also shows a transition between a 5 × 10 −8 and 2 × 10 −5 O 2 mixing ratio, but this model uses the modern outgassing fluxes for CO, H 2 , H 2 S, and SO 2 listed in Table 1 instead of presumptive Archean outgassing values. Also, at t = 0 y, we impose a stepwise increase of the O 2 flux from 10 12 to 1.8 × 10 12 molecules per cm 2 · s −1 while keeping the CH 4 flux/O 2 flux ratio at 0.45 (see SI Appendix, Fig. S1 for context). While the anoxic-to-oxic transition itself still occurs rapidly, the atmosphere simulated in Fig. 2B takes 60,000 y to reach the tipping point, which is much longer than the comparable O 2 transition shown in Fig. 2A.
The time required for O 2 to begin to rise in concentration is controlled by the reservoir of reducing gases, primarily CH 4 , H 2 , and CO, in the preoxygenated atmosphere. Big reservoirs of reducing gases slow the timescale of oxygenation, because reducing gases must be mostly removed before O 2 can increase. O 2 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 O 2 begins to rise in Fig. 2A, and why 60,000 y elapse before O 2 rises in Fig. 2B. Fig. 2B starts with more reducing gases, which take longer to eradicate.
We can roughly estimate the time required for O 2 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, N reducing is the O-equivalent column abundance of reducing gases (O equiv molecules per cm 2 ) which is equal to the sum of all reducing gases in the atmosphere (N j ) multiplied by α j , 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) Fig. S1). Transition in C results from a step-wise decrease in the CH 4 flux from 4.9 × 10 11 to 4.5 × 10 11 molecules per cm 2 · s −1 (black arrow in Fig. 1).
α C = −2, from redox stoichiometry of hydrogen, sulfur, and carbon, respectively. It then becomes straightforward to calculate the α j 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). N reducing is approximately equal to the weighted sum of CH 4 , H 2 , 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 A and B, we initiate the rise of oxygen by changing the surface flux of CH 4 and/or O 2 flux. We can quantify this flux perturbation in units of O equiv molecules per square centimeter per second (ΔF Oequiv ), .

[4]
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 O 2 transition in Fig. 2B yields τ oxy = 2, 900 y, a value only slightly smaller than the 3,500 y predicted by the time-dependent photochemical model. For Fig. 2B, τ oxy = 29, 000 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 (SI Appendix). 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. 2C 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. 2A, except we decrease the methane flux by twice as much at t = 0, from 4.9 × 10 11 molecules per cm 2 · s −1 to 4.5 × 10 11 molecules per cm 2 · s −1 instead of 4.7 × 10 11 molecules per cm 2 · s −1 (we keep the surface O 2 flux constant at 10 12 molecules per cm 2 · s −1 ). O 2 begins to rise and eliminates S 8 production after ∼1,500 y, but O 2 will reach higher levels because of the lower CH 4 flux. It takes ∼300,000 y for O 2 to reach its final steady-state abundance of 4 × 10 −3 mixing ratio. While the switch from 10 −8 to 10 −5 O 2 mixing ratio remains as rapid as in Fig. 2A, the predicted increase in O 2 concentrations to 4 × 10 −3 requires far longer. This timescale is roughly analogous to the time required to deplete H 2 and CH 4 reservoirs to allow O 2 to initially rise in concentration.
In summary, the timescale for O 2 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 O 2 to rise from 10 −8 to 10 −5 , reducing gases must first be removed, which can take thousands to 10 thousands of years ( Fig. 2 A and B). Increasing O 2 concentrations beyond 10 −5 to near percentage levels requires filling a large O 2 reservoir, which occurs on 10 5 -y timescales in our model run (Fig. 2C).
The Time Required for Deoxygenation. Here, we use our timedependent 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 O 2 and CH 4 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.
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 O 2 to begin to rise is merely the time required to deplete A B C Fig. 3. (A-C) Simulated reversal of the oxic transitions shown in Fig. 2 A, B, and C, respectively. Each oxic-to-anoxic transition is caused by a stepwise change of the CH 4 flux and O 2 flux at t = 0 y. the reservoirs of CH 4 and other reducing gases. Analogously, the timescale of deoxygenation is determined by the reservoir of O 2 and other oxidizing gases in the oxygenated atmosphere. The reversal shown in Fig. 3B starts with only 2 × 10 −5 O 2 , which can be depleted very quickly, allowing the return of an anoxic atmosphere. In contrast, the reversal shown in Fig. 3C takes 40,000 y because the atmosphere starts with 4 × 10 −3 O 2 , 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 O 2 changes. When reservoirs are big, for similar flux perturbations, the O 2 mixing ratio will change relatively slowly over time; however, when reservoirs are small, photochemistry permits rapid O 2 transitions. Therefore, the abundance of redox gases in an atmosphere is closely linked to the photochemical stability of oxygen. Fig. 4A shows the steady-state inertial timescale of redox gases, τ inertia , over the same axes as Fig. 1, which shows mixing ratios. τ inertia 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 O 2 flux, We choose the characteristic flux perturbation to be 10% of the O 2 flux because it is the same order of magnitude as natural redox variations that occur on modern Earth during Milankovitch cycling (see Discussion). An upper limit for the characteristic flux perturbation would be 100% of the O 2 flux. This would decrease all τ inertia values in Fig. 4A by a factor of 10, which would not change our interpretation. Since CH 4 , CO, H 2 , and O 2 are the most important redox gases, the numerator in Eq. 6 is well approximated by 4N CH4 + N CO + N H2 + 2N O2 . Oxygen is the most prone to change for the smallest τ inertia values, coinciding with O 2 mixing ratios between ∼10 −8 and ∼10 −5 shown in the whitish region of Fig. 4A.
The time-dependent photochemical models shown in Fig. 4 B-D illustrate the relationship between τ inertia and O 2 instability. To produce Fig. 4B, we started with the steady-state atmosphere indicated on Fig. 4A, then imposed 17% amplitude oscillations to the CH 4 flux with a period of 10,000 y. This forcing had no perceptible effect on the 3 × 10 −9 atmospheric O 2 . A similar 20% CH 4 flux oscillation also did not significantly perturb an oxic atmosphere starting with 3 × 10 −4 O 2 (Fig. 4D). However, just 5% CH 4 flux oscillations cause approximately four-ordersof-magnitude oscillations in surface oxygen concentrations for an incipiently oxic atmosphere starting with 3 × 10 −7 O 2 (Fig.  4C). O 2 is most unstable where the abundance of all redox gases is smallest relative to a characteristic redox surface flux (the whitish area of Fig. 4A) between ∼10 −8 and ∼10 −5 O 2 mixing ratio. Stability continually increases outside of this range of O 2 concentrations.
While O 2 was relatively stable in the Fig. 4 B and D simulations, it does not mean these atmospheres and initial oxygen concentrations are stable to all atmospheric perturbations. The stability of any O 2 mixing ratio depends on the atmospheric forcings that are likely in nature. In Discussion, we argue that the CH 4 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 ∼10 y are required to significantly affect O 2 concentrations. Imposing 100% amplitude fluctuations to the CH 4 flux with a period of 1 y, starting with the same atmosphere as Fig. 4C, did not significantly alter the atmosphere over time. Atmospheres with between ∼10 −8 and ∼10 −5 O 2 contain some CH 4 and O 2 , which gives the atmosphere inertia against annual to decadal change. 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 10 −6 O 2 mixing ratio, does not mean that such an atmosphere is stable and realistic over 10 5 -to 10 6 -y timescales. Gas fluxes from Earth's surface can vary during these timescales and significantly change O 2 concentrations (Results).

Discussion
For example, in the past 650,000 y, the biogenic methane flux has oscillated with an amplitude of 25% (6 × 10 9 molecules per cm 2 · 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 CH 4 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 CH 4 source to the atmosphere. Ice sheets covering and uncovering microbial mats could have affected global CH 4 fluxes. Fig. 4C illustrates the effect of 5% methane flux variations over Milankovitch timescales on an atmosphere starting with 3 × 10 −7 O 2 . O 2 oscillates nearly four orders of magnitude between ∼10 −8 (anoxic) and ∼10 −4 (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 SI Appendix, we show that the H 2 and CO outgassed from a significant LIP eruption could cause the O 2 surface mixing ratio to drop from 2 × 10 −5 to 4 × 10 −9 in ∼100 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 O 2 are unlikely. Additionally, our photochemical modeling shows that, for atmospheres with transitional O 2 concentrations, relatively small atmospheric perturbations (e.g., a CH 4 flux change of 5%) over timescales as short as hundreds of years can cause O 2 to change by orders of magnitude (e.g., Fig. 3B). Therefore, substantial variability of O 2 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 O 2 was unstable between 2.4 and 2.2 Ga. We find that shutoff of S 8 aerosol production, which is required to produce sulfur isotope MIF, occurs at ∼10 −7 O 2 mixing ratio, a region of the parameter space where O 2 is prone to rapid change (Fig. 4). But, oxygen surface levels between ∼10 −8 and ∼10 −4 mixing ratio were likely unstable. Short period changes to the biosphere, or volcanic outgassing rates, could have caused order of magnitude O 2 changes over 100-to 100,000-y timescales. Occasionally, big perturbations to the atmosphere, such as an LIP, might have lowered O 2 concentrations enough for sulfur isotope MIF to reoccur. Note that the above explanation for the cause of O 2 oscillations prior to 2.2 Ga is complicated by S-MIF data presented in Izon et al. (11), which do not suggest the same O 2 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 O 2 concentrations large enough to prevent O 2 collapse. Our modeling shows that larger O 2 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 O 2 could not have been larger than 2 × 10 −4 mixing ratio. Such a small O 2 reservoir would have been unstable to LIP eruptions, or variations in the CH 4 flux from Milankovitch cycles (Fig. 4), which both have evidence of occurring in the mid-Proterozoic stratigraphic record (40)(41)(42). We conclude that, for stability, mid-Proterozoic O 2 levels should have exceeded ∼10 −4 . This conclusion is compatible with mid-Proterozoic Fe isotopes in ironstones, which suggest O 2 levels between approximately 2 × 10 −4 and 2 × 10 −3 mixing ratio (43).
Our results are not sensitive to the changing solar UV photon flux between the GOE (∼2.4 Ga) and the mid-Proterozoic. Recalculating Fig. 1 using the solar UV flux at 1.3 Ga (44) results in surface O 2 and CH 4 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 O 2 was evenly distributed over the globe, then the surface O 2 mixing ratio was >3 × 10 −8 and <2 × 10 −7 , or 2) if O 2 accumulation was geographically restricted, then the O 2 surface flux was greater than 0.01 Tmol · y −1 (3 × 10 7 molecules per cm 2 · s −1 ). Our modeling suggests interpretation 1 is unlikely because we find that O 2 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 O 2 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 O 2 transitions (<10 5 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 O 2 levels.
An additional, related caveat is that our model does not consider biological feedbacks. The rise of O 2 would limit habitats for anaerobes, and permit more widespread aerobic respiration, potentially dropping the CH 4 flux/O 2 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 CH 4 and O 2 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 O 2 can begin to rise. This takes hundreds to 10 thousands of years. O 2 accumulation to just hundredths or tenths of percent levels requires filling a large O 2 reservoir, which may occur on a 10 5 -y timescale. Atmospheric deoxygenation occurs over similar periods of time, mainly controlled by the magnitude of the initial O 2 abundance.
We also find O 2 instability, especially for mixing ratios between ∼10 −8 and ∼10 −5 . For these O 2 concentrations, photochemistry demands that both CH 4 and O 2 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 3 × 10 −7 O 2 , 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 O 2 and the return of sulfur isotope MIF for an atmosphere starting with 2 × 10 −5 O 2 mixing ratio (SI Appendix).
We emphasize that the short-term (10 2 to 10 5 y) variability in O 2 levels considered here occurred on the backdrop of the billionyear 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 O 2 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 CH 4 flux/O 2 flux ratio. While this is evident between the Archean (CH 4 flux/O 2 flux ≈ 0.5) and modern (CH 4 flux/O 2 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 O 2 mixing ratio exceeding a value in the 10 −4 to 10 −3 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, The above system of partial differential equations (PDEs) describes how the number density (n i ) of each chemical species i changes over altitude and time.
In our photochemical model, we solve a simplified version of Eq. 7 which assumes that the total number density does not change over time (∂n/∂t ≈ 0). 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 i ) instead of number densities (see Appendix B.1 in ref. 33 for a derivation), [10] To approximate Eq. 8, 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. 8; 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. Thermal diffusion coefficient of species i. We neglect Dimensionless this term (α Ti = 0) T Temperature Kelvins The system of ODEs derived from finite differencing Eq. 8 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 CH 4 has a chemical lifetime of ∼10 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 (∂f i /∂t = 0). To find steady states, we begin with some initial atmospheric composition, then integrate Eq. 8 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 O 2 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. 8. 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 SI Appendix.
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. Data, Materials, and Software Availability. Source code has been deposited in Zenodo (https://zenodo.org/record/6824093#.Ys3KuOzMKCc) (48). PhotochemPy code and the corresponding Python scripts are available at https://zenodo.org/record/6824092. The most up-to-date code is available at https://github.com/Nicholaswogan/PhotochemPy (49).