Insights from a dynamical system approach into the history of atmospheric oxygenation

Domain: earth_science
Published: 2024-08-09
Editorial outcome: major_revision

Abstract
Atmospheric oxygen levels are traditionally viewed to have been relatively stable throughout Earth’s history with several-step increases. Emerging evidence, however, suggests extremely dynamic atmospheric oxygen levels through large swaths of Earth’s history. Here, we provide a new perspective on atmospheric oxygen evolution using a dynamical analysis to explore the relative importance of previously proposed feedbacks on the global oxygen and carbon cycles. Our results from a stochastic analysis of oxygen mass balance in this framework suggest there are multiple steady states for atmospheric oxygen, but only three stable states. One stable state under anoxic conditions (<10−5 present atmospheric level (PAL)), one at low (~10−3 to 10−2 PAL), and one near modern value atmospheric oxygen levels. Our findings also suggest two unstable states (tipping points) for atmospheric oxygen: one around 10−5 and another one around 10−1 PAL.

Introduction
Oxygenation of the Earth’s atmosphere was traditionally widely viewed to have occurred in two major steps: one in the late Archean and early Proterozoic (between 2.4–2.1 billion years ago (Ga)), known as the Great Oxidation Event (GOE), and another one at the end of the Proterozoic Eon between 0.8–0.6 Ga, referred to Neoproterozoic Oxygenation Event (NOE), coincident with the first major rise of complex, multicellular life1–3.

After the GOE, atmospheric oxygen levels are generally accepted to have passed the threshold (~10−6−10−5 present atmospheric level (PAL)4) for the formation of mass-independent fractionation of sulfur (MIF-S), as a geochemical proxy for an anoxic atmosphere5,6. And almost certainly reached values >1% PAL or higher during the geological interval known as the Lomagundi-Jatuli Event that records the most pronounced excursion in the isotopic composition of carbon in marine carbonate (δ13C)1,3,7,8. A wide range of models have been invoked to explain the GOE—from the emergence of cyanobacteria to changes in the cycling of bio-essential elements (e.g., phosphorous)9–11. Unlike the GOE, increases in atmospheric oxygen levels during the NOE are poorly constrained. Increases in biological complexity around this time have been commonly invoked as the driver or result of this oxygenation event (e.g., 12). Regardless of the causes of the GOE and NOE, an increase in atmospheric oxygen has been traditionally assumed to have been unidirectional with major steps occurring during the GOE and the NOE1. The rise of land plants has also been proposed to drive a step jump in atmospheric oxygen levels (e.g., 13), consistent with the results from both biogeochemical modeling and proxy records.

Recent evidence from geological records challenges the traditional view of stepwise atmospheric oxygenation. More specifically, results from a suite of geochemical proxies that capture the change in the surface ocean and atmospheric oxygen levels, point toward dynamic atmospheric oxygen14–19. For instance, analysis of Paleoproterozoic-age sedimentary rocks suggests fluctuations in the production of MIF-S during and after the GOE—implying swings in atmospheric oxygen levels14. In the Proterozoic Eon, the results from redox-sensitive elements and other geochemical proxies also suggest pulses of surface ocean and atmospheric oxygen levels15,16,20–22. Lower proxy-based estimates of atmospheric oxygen are on the low end or below estimates of stable oxygen states driven by 1-D atmospheric models (e.g., 23). Although our picture of Earth’s redox evolution is continually developing, it is becoming increasingly difficult to explain some of the recent insights from the rock record that suggest variable atmospheric oxygen (e.g., 14) and are inconsistent with the traditional view of largely invariant atmospheric oxygen levels during the Proterozoic (e.g., 1).

To provide a new view on the range of stable states of atmospheric oxygen, we conducted a dynamical analysis of global oxygen mass balance. Dynamical system analysis, in its simplest definition, is the study of the state of a system as its main parameter (initial condition, in this case, pO2) is varied. This tool is mainly useful for studying non-linear systems where a small change in a parameter in the system can lead to a substantial change in the state of the system. This approach has not yet been widely used in the Earth Sciences—and more specifically in studying the evolution and dynamic of Earth’s biogeochemical cycles of major elements. However, it is widely used to study non-linear systems across a wide range of fields such as engineering and economics and is ideally suited to explore systems with strong positive and negative feedbacks—like the global oxygen cycle.

We use a simple oxygen mass balance expression as the basis for our analysis:

d(pO2)/dt = JNPP + JHydrogen − JCH4,O2 − Jvolc − Jocean,sink − Jorg,weath − Jpyr,weath (1)

Here d(pO2)/dt is the change in atmospheric oxygen level in Tmol per year, JNPP, JHydrogen, is the source of oxygen respectively from the burial of organic carbon on land and in the ocean, and the rate of escape of hydrogen atoms to space. The production of oxygen through organic burial is the sum of oxygen production in the ocean and on land. The oceanic oxygen production is calculated by considering the strength of net primary production (NPP (Tmol year−1)), the burial efficiency of organic matter, and a factor that accounts for the effect of anoxia on phosphorous cycling and organic matter burial (see “Methods”). Similar to the previous studies, the terrestrial oxygen production includes the fire feedback which is parameterized as a function of pO2, carbon dioxide level (pCO2), and temperature. The rate of hydrogen escape is a function of atmospheric methane level (pCH4), following previous studies.

The JCH4,O2, Jvolc, Jocean,sink, Jorg,weath, Jpyr,weath, are sinks of oxygen which respectively represent atmospheric methane oxidation, volcanic outgassing (e.g., H2S), oceanic sinks of oxygen (oxidation of sulfide, ferrous iron, and methane), rate of oxidative weathering of organic matter, and pyrite. The rate of atmospheric methane oxidation (JCH4,O2) is similar to previous studies, and it is a function of pO2, pCH4, and an oxidation parameter (Ψ). The rate of oceanic sink of oxygen is parameterized to be a function of NPP and the amount of anoxia in the ocean, which is a function of pO2. Finally, consistent with previous modeling results, the oxidative weathering of pyrite and organic matter was assumed to change, albeit non-linearly, with pO2.

Tied to each of the fluxes presented in Eq. (1) are previously proposed feedbacks in the coupled carbon, oxygen, and sulfur cycles. For instance, oxygen fluxes from the biosphere are increased when surface oxygen levels are decreased—given that anoxia favors organic carbon burial. In contrast, fires can suppress atmospheric oxygen levels—higher atmospheric oxygen levels favor fires. All of the utilized feedbacks are tied to oxygen fluxes (full parameterization schemes used for the oxygen sinks and sources are presented in the “Method” section). Since the flux of oxygen corresponds to hydrogen escape (JHydrogen), and atmospheric methane oxidation (JCH4,O2) depends on the concentration of atmospheric methane, we coupled the oxygen mass balance to a simple mass balance model for atmospheric methane that includes the flux of methane from the ocean as a source of methane and atmospheric methane oxidation and hydrogen escape flux as the sinks of atmospheric methane (pCH4). The details of the parameterization used for the methane mass balance model are also explained in the “Methods” section.

The feedbacks we are exploring have been extensively discussed and there is general consensus that these factors play a critical role in modulating Earth’s atmospheric oxygen levels4,10,24–26. There is, however, uncertainty in the strength of each of these feedbacks. To account for these uncertainties, we employed a stochastic approach, exploring the range of previously proposed values for key terms in our oxygen mass balance4,10,23,24,27–29. Specifically, in our stochastic approach, instead of using a fixed value for input parameters, the values for the parameters will be sampled randomly within the range suggested in the literature (Table S1), assuming a uniform distribution, and the most probable range of dpO2/dt will be obtained for the full range of pO2 values. We did 1000 simulations for each of the stochastic runs. There was no noticeable difference in the model results with more than 1000 runs, which suggests that our selected number (1000) is sufficiently large to capture the variability in the range of model input parameters. The stochastic approach used in our analysis allows us to be relatively agnostic about the debate over specific parametrizations and explore possible Earth system states given the current consensus on the long-term oxygen cycle. Finally, while our model accounts for the major feedbacks on the global oxygen mass balance, there are other environmental factors (e.g., tectonic activity, continent configuration, etc.) that likely impact atmospheric oxygen levels and have varied throughout Earth’s history and also likely impacted some of the major feedbacks. However, our stochastic approach would likely capture, at least in part, the effect of these factors in altering the feedback strengths.

Results
Our results indicate five possible steady states for the atmospheric oxygen levels (dpO2/dt=0) (Fig. 1). Generally, stable states (states I, III, V; Figs. 1, 2) are defined as the steady-state points (dpO2/dt=0) where the sign of dpO2/dt changes from positive to negative. In these states, the system tends to remain in equilibrium because the positive-to-negative transition indicates a decrease in pO2 values when starting from values larger than the steady state, eventually reaching stability. Conversely, when starting from values smaller than the steady state, the positive-to-negative transition signifies an increase in pO2 values, pulling the system back toward equilibrium (states I, III, V; Figs. 1, 2). Unstable states are also steady-state points (dpO2/dt=0) but differ in that the sign of dpO2/dt changes from negative to positive (states II, IV; Figs. 1, 2). In these states, the system tends to move away from equilibrium, as the negative-to-positive transition implies an increase in pO2 values when starting from values larger than the steady state, leading the system away from stability. Conversely, when starting from values smaller than the steady state, the negative-to-positive transition indicates a decrease in pO2 values, pushing the system further from equilibrium.

Our results from the stochastic simulation were filtered to be consistent with the geochemical proxies. Specifically, we only considered the results from the stochastic simulation that include a steady state point under high oxygen modern-like conditions, a steady state point under anoxic conditions (pO2 < 10−5 present atmospheric level (PAL)) which is consistent with the presence of anoxic atmosphere during Earth’s early stage (e.g., mass-independent fractionation of sulfur), and a steady-state point between the two end members of high and anoxic conditions that would generally correspond to the geologic times where there is no evidence for anoxic nor high oxygen conditions.

The results from the stochastic analysis of the global oxygen mass balance suggest that atmospheric oxygen can be in balance under essentially anoxic conditions (state I; ~8×10−9 (25th percentile) to 2×10−7 (75th percentile) with median ~3×10−8 and minimum and maximum values of ~1×10−13 and 4×10−7 PAL), very low oxygen conditions (state II; 1×10−5 to 1×10−4 with median ~3×10−5 and minimum and maximum values of 3×10−7 and 2×10−4 PAL), low oxygen condition (state III: ~1×10−3 to 1×10−2 with median ~3×10−3 and minimum and maximum values of ~8×10−5 and 3×10−2 PAL), relatively high oxygen condition (state IV: ~1×10−1 to 4×10−1 with median ~2×10−1 and minimum and maximum values of ~10−2 and 7×10−1 PAL), and high, near modern oxygen conditions (state V: 9×10−1 to 1.1×10° with median ~10° and minimum and maximum values of 7×10−1 and 1.3×10° PAL) (Fig. 1b). To explore the extent of stability in each stable state, we calculated the energy corresponding to each steady state point. The energy was calculated as the surface area under the curve (Fig. 1a) which is the integral of dpO2/dt over the range of pO2 (∫ pO2,f pO2,i dpO2/dt dO2). The amount of energy depends on the extent of the so-called basin of attraction (a region in which the initial conditions (pO2) lead to the eventual convergence of the system to a specific stable steady state or attractor). For instance, state I with the largest basin of attraction has the minimum energy implying large stability around this state (Fig. 2).

According to our results, the only stable states occur under states I, III, and V. The first stable state, state I, emerges under anoxic conditions, followed by an unstable state (state II) with an atmospheric oxygen value around ~10−5 PAL. State III is another stable state which is followed by an unstable state (state IV). Under the unstable states, a small perturbation in the global oxygen mass balance would lead to a shift toward adjacent stable states (states I, III, and state V). The expression of the state shifts would likely be high levels of variance in atmospheric oxygen levels around a relatively stable baseline condition. The results also show that modern-like high oxygen conditions of state V can only be achieved if the atmospheric oxygen level passes the pO2 threshold for the unstable state IV (~10−1 PAL) (Figs. 1, 2). This means that there is a rather large “basin of attraction” for State III which suggests that under a wide range of parameter values (Table S1) that can represent a wide range of environmental conditions and a wide range of initial conditions (pO2 values), the system moves dynamically toward State III. The basin of attraction for State V starts from the maximum pO2 at State IV, which means that reaching State V requires an increase in the pO2 level above this value.

To further explore the influence of each feedback on the global oxygen mass balance, we conducted a sensitivity analysis. Specifically, to test the sensitivity of our results for the steady-state pO2 levels to the strength of each feedback, we fixed the value of each oxygen feedback to its maximum value used in the stochastic analysis and recorded the change in the median of dpO2/dt vs pO2 (Figure S2). The magnitude of change in the value for a steady-state point for pO2 in response to a change in the rate of oxygen production or consumption indicates the relative importance of that feedback in driving the steady-state point for pO2. For instance, according to Fig. 2S, a change in the rate of hydrogen escape would result in a major shift in the value of the first steady-state point for pO2 with no substantial change on other steady-state points for pO2. This result suggests an important role of hydrogen escape in driving the first pO2 steady-state point. According to the results from the sensitivity analysis, the rate of hydrogen escape and the atmospheric oxidation of methane have the largest influence on the emergence of the first and second steady-state for pO2 under anoxic conditions (Figure S2), consistent with the previous work30. Specifically, a change in the rate of hydrogen escape and atmospheric methane oxidation can result in a shift or disappearance of the first and second steady-state points for atmospheric oxygen levels (Figure S2).

The major role of the rate of hydrogen escape on the first stable state for the atmospheric oxygen level (state I) indicates the key role of hydrogen escape in stabilizing the atmospheric oxygen level on the early Earth. According to our results, an increase in the rate of atmospheric methane oxidation influences the emergence of the second steady-state point (state II) around 10−5 PAL, consistent with the modeling results from previous studies (e.g., 30).

Based on our results from the sensitivity test, the emergence of the third and fourth steady-state points under low-oxygen conditions (states III & IV) is controlled by the oceanic oxygen production (JNPP) and the sink of oxygen in the ocean and on land. An increased rate of oceanic oxygen production (JNPP) would shift the imbalance value of oxygen (dpO2/dt) to positive values, which would influence the appearance of steady-state points under low-oxygen conditions. At the same time, according to our results, while under anoxic conditions (<10−5 PAL) the oxidative weathering of pyrite and organic matter is essentially ineffective in controlling the global oxygen mass balance (Figures S2, S3), an increased rate of oxidative weathering under low-oxygen condition would dramatically change the oxygen imbalance value which would impact the emergence of steady-state points under low-oxygen conditions. Of the two steady-state points under low-oxygen conditions (states III & IV), the first one is the result of an increase in oxygen production in the wake of the evolution of oxygenic photosynthesis while the rate of oxidative weathering is low. An increase in the atmospheric oxygen level would result in an increase in the rate of oxidative weathering which would result in a negative imbalance for the oxygen. A continued increase in the production of oxygen in the ocean would have resulted in a decrease in anoxia in the ocean which would have acted as positive feedback in the global oxygen mass balance and have resulted in the emergence of the state IV. Such a state is, however, unstable, which means that a change in the oceanic oxygen mass balance (e.g., a change in the oxygen production through net primary production or marine oxygen consumption) would shift the system to the adjacent stable state (state III or state V), depending on the directionality of the change. Our results also suggest a minor role of enhanced phosphorous recycling under anoxic conditions in driving the emergence of steady-state points under low-oxygen conditions, caused by low-sulfate conditions under low-oxygen conditions.

Our new results on the potential for multiple stable states can be placed in the context of widely held views of atmospheric oxygen evolution. For instance, some foundational modeling work has suggested the presence of a bi-stability in the atmospheric oxygen level, with one stable state less than 10−5 PAL, another stable state above 5×10−3 PAL, and an unstable state in between these two stable states30. Our results are consistent with this work—where there are two stable states (states I and III) before and after around 10−5 PAL with an unstable state around 10−5 PAL (state II; Fig. 1). Our results for the range of atmospheric oxygen at unstable state II (1×10−5 to 1×10−4 with median ~3×10−5 and minimum and maximum values of 3×10−7 and 2×10−4 PAL) are largely consistent with the previous modeling results. The results from the photochemical modeling suggest instability of atmospheric oxygen between 5×10−8 and 5×10−4 PAL that are largely consistent with our results for the range of atmospheric oxygen level at unstable state II23.

Our results for the range of stable states under low-oxygen conditions is also consistent with the some previous results. For instance, the results from a one-dimensional time-dependent photochemical model suggest a stable range of atmospheric oxygen above ~5×10−4 to 5×10−3 PAL23. This range is consistent with the range of atmospheric oxygen resulting from our stochastic simulation for state III (state III: ~1×10−3 to 1×10−2 with median ~3×10−3). The results from another modeling study also indicated the presence of three stable states, one under anoxic conditions, one under low-oxygen conditions, and one under high-oxygen conditions27, which is roughly consistent with our results here. However, these studies did not focus on the possibility and extent of the fluctuation in the atmospheric oxygen level under low-oxygen conditions, which is increasingly supported by recent lines of evidence (e.g., 14,31). Our results suggest the Earth system can reach a steady state at high oxygen conditions (state IV; ~10−1 PAL), as suggested by previous work27, but that such a state is largely unstable (Fig. 2). Our results are consistent with the range of steady-state oxygen levels suggested by the previous studies. However, our analysis can provide further insight into the nature of such steady-state points, the range of variations in atmospheric conditions within these Earth system states, and the transitions between the states.

We can directly compare these stable states to estimates of atmospheric oxygen levels from the rock record. Prior to the rise of cyanobacteria, the escape of hydrogen atoms to space was the only source of oxygen in the atmosphere32. This oxygen production pathway would have been significant under very low oxygen conditions, especially if there were high atmospheric methane levels32. The increased pO2 level caused by significant hydrogen escape flux would result in a decrease in the atmospheric methane level (pCH4) which would, in turn, result in a reduction in the hydrogen escape flux and lead to the first stable state of atmospheric oxygen (State I; Fig. 1). The wide range of values for state I also suggests the possibility of large fluctuations and dynamic conditions in the atmospheric oxygen level under anoxic conditions.

The evolution of cyanobacteria would have added another source of oxygen to the global oxygen mass balance, resulting in a shift to a higher atmospheric oxygen level. Increased oxygen production is likely to have transitioned the Earth system to the second steady state point of atmospheric oxygen (State II; Fig. 1). Such a state, however, unlike the first state, is unstable, and the system would only occupy such state transiently (e.g., a tipping point) and a small change in the rate of oxygen production or consumption would result in a transition into the adjacent stable states (state I and state III).

According to our results, oscillations in MIF-S production, suggested by the recent results from geochemical proxies14,33, can be explained by a fluctuation in the atmospheric oxygen level in both states I and II. Specifically, the maximum value in state I is ~2×10−7 PAL which is higher than the minimum value suggested for the threshold of atmospheric oxygen for caseation in the production of MIF-S (e.g., 23) which implies that it is possible that a fluctuation in the atmospheric oxygen level under state I result in an oscillation in the production of MIF-S signal. Suggested higher values for the threshold of atmospheric oxygen level for MIF-S formation (~10−6 PAL) are also consistent with the range of atmospheric oxygen level at state II which suggests that the apparent fluctuation in MIF-S signal can alternatively be attributed to an oscillation in the atmospheric oxygen level at unstable state II. Regardless, according to this framework, the GOE can then be defined as a departure from unstable State II and the transition to the second stable state for atmospheric oxygen level (state III: Fig. 1).

Possible large fluctuation under low-oxygen conditions

The Proterozoic, post-GOE world would have been trapped under low oxygen conditions with the possibility of large swings in atmospheric oxygen level. The absence of evidence for MIF-S for Proterozoic rock record (early Paleoproterozoic excepted) and the absence of evidence for modern-like oxygen conditions, suggest that atmospheric oxygen during the Proterozoic would have been between 10−5 PAL to modern atmospheric oxygen level (Fig. 3). Within this wide possible range for atmospheric oxygen levels, our results show only two steady states with only one stable state (State III: ~10−3 – 10−2 PAL) (Fig. 3B). Our stochastic simulation suggests that at stable State III atmospheric oxygen level can vary approximately more than two orders of magnitude (~minimum and maximum of ~8×10−5 and 3×10−2 PAL), which suggests that atmospheric oxygen levels can take any values within this range without falling out of the stability zone (Fig. 1). This result opens the possibility for orders of magnitude swings in the atmospheric oxygen level under low-oxygen conditions. According to our results, the evidence for higher oxygen values beyond the ~10−2 PAL, as the upper threshold for oxygen level at State III, can only be explained by reaching State IV. This state, however, is unstable, meaning that relatively high oxygen conditions at State IV would have been ephemeral and a small change in the rate of oxygen production or consumption in the ocean and on land would force the system to move dynamically to adjacent stable states (state III or state V; Figs. 1, 2). Taken together, our finding suggests that there could have been large fluctuations in atmospheric oxygen levels during the Proterozoic, likely within the range of about 10−4 to 10−2 PAL (Fig. 3B). This conclusion is, arguably, in accordance with some geological evidence (e.g., sulfur isotope systematics14 and Cr isotope values34, among other redox proxies1).

According to our results, high-oxygen conditions during the Lomagundi-Jatuli event (LJE) can be attributed to two stable states of III and V. Specifically, considering that there is no strong evidence to date for large oscillation of atmospheric oxygen during the LJE that can be linked to an unstable state of atmospheric oxygen (e.g., state IV), the suggested high oxygen level during LJE can be explained by either an atmospheric oxygen level in state III or V. The maximum value of atmospheric oxygen in state III is around 10−2 PAL which suggests that an increase in the rate of oxygen production during LJE can be explained by a rise in atmospheric oxygen level to the maximum value of state III. Alternatively, a high atmospheric oxygen level during this time can be explained by reaching a modern-like atmospheric oxygen level (state V). Regardless, our results suggest large shifts in atmospheric oxygen levels across the LJE transition.

Late Proterozoic and transition to the modern-like pO2 level

The late Proterozoic and early Phanerozoic times marked the end of the large fluctuation in the atmospheric oxygen level and the transition into the last stable state of atmospheric oxygen. Specifically, our results suggest that an increase in pO2 levels would have resulted in moving from the stable State III to the stable State V with the modern-like pO2 value. This increase in oxygen levels could have been caused by either a decrease in oxygen sinks in the ocean or land or an increase in the oceanic oxygen production caused by enhanced marine primary production35,36, which is consistent with the results from geochemical proxies that are suggestive of a transition into roughly modern-like atmospheric oxygen level around this time37.

Strong stabilizing feedback of plant evolution on pO2

The evolution of plants in the Phanerozoic likely helped stabilize atmospheric oxygen around modern value. To explore the impact of plant evolution on the atmospheric oxygen level, we conducted two stochastic simulations with and without plants. We explored the effect of the rise of plants by removing the flux of carbon burial from land plants from the total oxygen production flux (JNPP; Eq. 2 in the “Methods” section). According to our results, the median atmospheric oxygen level without plants would be around 0.4 atmospheric oxygen level (PAL). This result shows that the rise of plants resulted in an increase in the oxygen production flux and helped with stabilizing the atmospheric oxygen level near modern value (Fig. 4). Similar to the results from previous studies (e.g., 13), our results also suggest that, by exerting negative feedback such as fire, plants may also help prevent large swings in atmospheric oxygen levels.

Taken together, our stochastic results from dynamical analysis of global oxygen mass balance suggest three stable states for atmospheric oxygen level: one under anoxic conditions ~10−8 PAL, one around 3×10−3 PAL, and one near modern value (~10° PAL). According to our results, the atmospheric oxygen level could have fluctuated greatly for most of the Earth’s history until the transition into the last stable state in the late Proterozoic and early Phanerozoic, where the oxygen level finally reached its modern value. The transition into state V with the modern-like oxygen level is coincident with the major rise and radiation of complex life, reinforcing the need to explore the role of instability in the global oxygen cycle played in triggered biological innovations in early animals38.

Methods
Global oxygen mass balance
To explore the range of stable states for the atmospheric oxygen level, we used a simple expression for the global oxygen mass balance:

d(pO2)/dt = JNPP + JHydrogen − JCH4,O2 − Jvolc − Jocean,sink − Jorg,weath − Jpyr,weath (2)

Here d(pO2)/dt is the change in atmospheric oxygen level in the unit of Tmol per year, JNPP, JHydrogen, are the sources of oxygen respectively from burial of organic carbon on land and in the ocean, and the rate of escape of hydrogen atoms to space. The JCH4,O2, Jvolc, Jocean,sink, Jorg,weath, Jpyr,weath, are sinks of oxygen which respectively represent methane oxidation in the atmosphere, volcanic outgassing, the oceanic sink of oxygen (oxidation with sulfide, iron, and methane), rate of oxidative weathering of organic matter and pyrite. Below we detail the parameterization scheme used for each parameter in Eq. (1). It should also be noted that to account for the uncertainty associated with the choice of values of the model input parameters, we took a stochastic approach wherein the values for each modeling input parameter were selected randomly assuming a uniform distribution and the most probable range of results were obtained (Fig. 1b).

JNPP. The JNPP as a total carbon burial on land and in the ocean can be expressed as:

JNPP = JNPP,ocean + JNPP,land (3)

The flux of carbon burial in the ocean, JNPP,ocean, can be parameterized as:

JNPP,ocean = ε * NPP * area * Rc/p (4)

Where ε is the efficiency of the carbon pump, NPP is the net primary production, and the last term is the area of the ocean. The value of NPP is multiplied by a productivity factor to account for a lower productivity at lower oxygen levels39,40:

NPP = NPPmodern * f (5)

We used high burial efficiency under anoxic conditions, consistent with the modeling and experimental results as well as the modern observations that indicate an increase in the burial efficiency of carbon under anoxic conditions41,42. The fraction of oxic condition is calculated as one minus the fraction of anoxia (1- fanoxic). We considered full anoxia (fanoxic=1) under very low oxygen conditions (pO2 = 10−13 PAL), consistent with results from geochemical proxies, and modern-like anoxia for modern-like pO2 (fanoxic ~0.0001 to 0.001). The fraction of anoxia for pO2 less than 10−5 PAL was less than 0.05, consistent with the results from geochemical proxies and the biogeochemical modeling that suggest weak ocean oxygenation under the atmospheric oxygen threshold for production of mass-independent fractionation of sulfur (e.g., 43,44). The fraction of anoxia in the ocean for other initial conditions (pO2) was varied between these two end members. To account for the uncertainty involved in choosing the modeling parameters, the value for each parameter was selected in a stochastic manner with the assumption of uniform distribution for the range of each parameter.

We also assume lower rates of oxygen production at lower atmospheric oxygen levels, following previous studies (e.g., 27). Considering a simple mass balance framework for atmospheric oxygen which considers oxygen production through organic matter burial and oxygen consumption through the reaction of oxygen with reduced species as well as the oxidative weathering of pyrite and organic matter (JO2,production = JO2,sink). Given that the rate of oxygen consumption is a function of atmospheric oxygen level and at lower atmospheric oxygen level, the rate of oxygen consumption would be lower too (e.g., oxidative weathering of organic matter which is further discussed below), to maintain a low-atmospheric condition, the rate of oxygen production would also need to be lower at lower atmospheric oxygen level. This relationship is like shifts in nutrient availability under low-oxygen conditions (e.g., 27).

The term Rc/p is the ratio of organic carbon to organic phosphorous normalized to modern (oxic) conditions. The term Rc/p accounts for the effect of change in redox condition and expansion of anoxia on phosphorous recycling and organic matter burial. Specifically, under anoxic conditions, the release of phosphorous would be enhanced, which would result in an increase in marine primary production and organic matter burial11,45. The link between anoxia and enhanced phosphorous availability however depends on the availability of sulfide which controls the availability of dissolved iron (e.g., 46). Given the interaction between phosphorous and dissolved iron, an increase in sulfide availability would result in a decrease in the iron sink of phosphorous (e.g., 46). Thus, an expansion in the anoxia can foster a higher rate of organic matter burial only if the rate of sulfide production through sulfate reduction is high which requires high sulfate levels. We followed the previous studies11,45 to parameterize the term Rc/p and included a factor that accounts for sulfate availability as well:

Rc/p = [ (C/P)anoxic * (fanoxic) + (C/P)oxic * (1 - fanoxic) ] / (C/P)oxic * RSO4 (6)

Here (C/P)anoxic, and (C/P)oxic, are the C to P ratios under anoxic, and oxic conditions suggested by the previous studies ((C/P)anoxic = 1100; (C/P)oxic = 250)11,45. The [O2] and [O2]modern are the seawater oxygen concentration and modern level, respectively. We used the previous earth system modeling that relate the seawater oxygen concentration to various pO2 levels47.

The term RSO4 accounts for the effect of sulfate level on phosphorous availability, which is parameterized as a Monod-type formulation, similar to the parameterization scheme suggested for the rate of sulfate reduction:

RSO4 = [SO4] / ([SO4] + km,SO4) (7)

Where [SO4] and km,SO4 denote the seawater sulfate concentration and the half-saturation constant for sulfate reduction. We assumed sub-mM levels for seawater sulfate under low oxygen conditions (pO2 < 10−2 PAL) and mM range and super-mM range of seawater sulfate for higher pO2 levels, consistent with the range of seawater sulfate suggested by the previous studies under low pO2 and high pO2 levels.

We utilized a parameterization from previous studies to calculate the flux of carbon burial on land13,48:

JNPP,land = 2 * E * VNPP(pO2) * VNPP(T) * VNPP(CO2) * kfire * fignit (8)

In this formulation, the strength of carbon burial on land is a function of pO2, temperature (VNPP(T)), and carbon dioxide (VNPP(CO2)). The term E is the rate of carbon burial on land, whose value was from the previous works and factor 2 is for normalization13,48. The last term in the equation corresponds to the fire feedback on organic burial on land. Following the previous works, the terms in Eq. (5) were parameterized as follows:

VNPP(pO2) = 1.5 - 0.5 * pO2 (9)

VNPP(T) = 1 - ((T - 25) / 25)^2 (10)

VNPP(CO2) = (pCO2 - pCO2,min) / (pCO2,1 + (pCO2 - pCO2,min)) (11)

Where T is the temperature in the unit of degrees Celsius, pCO2 is CO2 atmospheric level in ppm, pCO2,min denotes the minimum atmospheric CO2 concentration in which plants can grow, and pCO2,1 is a constant13,29,48,49.

In the fire feedback term, kfire is a constant and term “ignit” is a linear function of pO2:

ignit = max(586.2 * pO2 - 122.102, 0) (12)

JHydrogen and JCH4,O2. Following previous works, we considered the escape flux of hydrogen atoms to space to be a function of methane partial pressure (pCH4)4,30,32:

JHydrogen = α * μ * pCH4 (13)

Where α is a constant and μ denotes the number of moles in 1 atmosphere. It should be noted that the number of moles in 1 atm is constant. The atmospheric concentration of methane, pCH4, would be calculated by solving the following methane mass balance equation:

d(pCH4)/dt = JCH4 - JCH4,O2 - JHydrogen (14)

Here JCH4 is the net flux of biogenic methane from the oceans, JCH4,O2 is the rate of atmospheric methane oxidation, and JHydrogen, as explained above, is the hydrogen escape flux.

The net flux of biogenic methane from the ocean is considered to be a function of the net primary production which regulates the amount of organic matter available for methanogenesis, the extent of anoxia in the ocean, and a coefficient, which reflect the amount of methane being oxidized in the ocean which itself is the function of oxygen availability. At higher atmospheric oxygen, the value of this coefficient increases, reflecting an increase in the oceanic oxidant reservoir (e.g., sulfate and oxygen concentrations) which would impact the rate of methane oxidation in the ocean.

JCH4 = α1 * (1 - ε) * NPP * (fanoxic * area) (15)

Parameters α1, ε, NPP, and fanoxic, respectively, denote the coefficient for methane production, burial efficiency of organic matter, net primary production, and the areal fraction of ocean anoxia.

We used the previous parameterization to simulate the rate of atmospheric methane oxidation30. Specifically, the rate of methane oxidation is assumed to be a function of the partial pressure of both atmospheric oxygen and methane multiplied by an oxidation coefficient.

JCH4,O2 = Ψ * μpO2 * μpCH4^0.7 (16)

Here μ denotes the number of moles in 1 atmosphere. The oxidation parameter, Ψ, is parameterized as logΨ = p1x^4 + p2x^3 + p3x^2 + p4x^1 + p5, where x = log(μ pO2) and coefficients p1 = 0.00286, p2 = −0.1608, p3 = 3.216, p4 = −26.53, p5 = 76.04 are resulted from fitting a curve to the previously published data (Figure S1). The details of the goodness of the fit, along with the range of values and their confidence intervals are shown in Figure S1.

While our parameterization scheme for the rate of atmospheric methane oxidation is similar to previous studies, some assumptions in this parameterization can potentially cause uncertainty in the rate of methane oxidation. For instance, this parameterization makes assumptions about the age of the Sun as well as about the surface temperature, both of which effects can potentially influence the photochemistry of CH4 and O2. The assumption about the age of the sun can impact the Lyman alpha UV flux, an important factor for kinetics of CH4 photolysis and the Earth’s surface temperature can impact the water vapor content of the troposphere with impacts on the abundance of OH and the rate of methane oxidation. To account for the uncertainty associated with these assumptions, we varied the kinetics of atmospheric methane oxidation in our stochastic simulation.

Finally, considering that the methane mass balance equation depends on pO2, we first solved the methane mass balance equation and obtained the fixed points (steady states) of the methane mass balance equation, and then used those fixed points to solve the oxygen mass balance equation.

Jocean,sink. The oceanic sink of oxygen was a function of the strength of marine net primary production and fraction of anoxia, as well as the seawater oxygen concentration (e.g., 27,45,50).

Jocean,sink = k1 * (1 - ε) * NPP * (fanoxic * area) * (1 - [O2] / ([O2] + km,O2)) (17)

Where [O2] and km,O2 denote the seawater oxygen concentration and the half-saturation constant for aerobic respiration. As explained above, we used the results from previous earth system modeling to relate the pO2 to seawater oxygen concentrations. The constant “k1” is a factor that accounts for variability in the pool of reducing compounds (e.g., sulfide, iron, methane) in the ocean. The constant “k1” reflects the relative comparison between the rate of oxygen consumption by reducing species. For instance, a value of 0.1 means that the rate of oxygen consumption by reducing species is about an order of magnitude lower than the maximum rate of anaerobic respiration. Under anoxic and low-oxygen conditions of Archean and Proterozoic, this value would have likely been low. For instance, most of the sulfide resulting from sulfate reduction would have precipitated as pyrite, consistent with the results from sulfur isotope records that suggest a near-one value for a fraction of pyrite (fpyr). Our results for the Jocean,sink for the modern-like condition are consistent with the rate of oxygen uptake in the modern oxygen minimum zone. Specifically, using our formulation for the rate of oxygen consumption (Eq. 16 in the “methods” section), and considering the modern organic matter burial efficiency (~0.005), the modern net primary production (~100 gram/m2/year), and the modern areal fraction of anoxia (~0.005), the order of magnitude rate estimate for the rate of oxygen consumption would be around 0.02 teramol per year. And using the measured rate of sulfate reduction in the modern oxygen minimum zones which is on the order of 0.1 mmol/m2/day, and the areal fraction of anoxia in the ocean and considering that half of the reduced sulfur would become oxidized with oxygen, the rate of oxygen consumption through reaction with sulfide would become around 0.03 tera per mol per year which is comparable with the value resulted from our estimate. Finally, to account for the uncertainty, this value along with other values used in this equation and elsewhere was varied within a wide range in our stochastic simulation (Table S1).

Jorg,weath and Jpyr,weath. Consistent with previous modeling results, the oxidative weathering of pyrite and organic matter was assumed to change with pO2 28,29,48,51. The previous modeling results show that at a certain pO2 level, the oxidative weathering of organic matter and pyrite would be complete and equal to their modern values. This result suggests that above a certain pO2 level, the weathering flux would no longer be controlled by pO2. With this in mind, we assumed that the weathering flux of organic matter and pyrite changes with pO2 followed a Michaelis-Menten equation:

Forg,weath = Fweath1 * (pO2 / (kweath1 + pO2)) * (NPP / NPPmodern) (18)

Fpyr,weath = Fweath1 * (pO2 / (kweath1 + pO2)) (19)

Here Fweath1 denotes the modern values for pyrite and organic matter weathering, and kweath1 is the atmospheric oxygen level above which the weathering of pyrite and organic matter is complete (Table S1). To account for the effect of the amount of organic matter produced in the ocean, the oxidative rate of organic matter would be multiplied by the ratio of NPP to its modern value (NPP/NPPmodern). Finally, to account for the variability in choosing the parameters, we took a stochastic approach where modeling parameters were randomly sampled within a specified range supported by literature (Table S1), and the most probable range of stable states of pO2 was calculated.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "## Summary\nThe paper constructs a dynamic stochastic model of the oxygen concentration on Earth. Their\nmodel is a simple ordinary differential equation describing the rate of change of O2 levels\n($dp_\\mathrm{O_2}/dt$). For a wide range of input parameters, the find the steady-state O2\nconcentrations between 1e-13 and 1 present atmospheric levels (PAL). The nominal, median model\nfind 5 steady-states. Three of the steady-states are stable, while 2 are not stable. The authors then\nargue that the history of oxygen on Earth might be understood in terms of these steady-states. One\nstable steady state occurs for Archean amounts of O2 (~3e-6 PAL), Another stable steady state\noccurs at ~1e-3 PAL O2, perhaps representative of the Proterozoic eon. Finally, the modern Earth is\nat the final stable steady state at ~ 1 PAL O2.\n## Overall comments\nOverall, the paper is fairly well written, concise, and generally makes sense. It is a valuable addition\nto the literature because considers a new stochastic dynamic approach to determine what O2\nconcentrations are most likely to be stable. This is a new contribution to the literature, based on my\nknowledge. I have a few relatively minor comments that I would like to see addressed before\npublication.\nThe article argues for five steady-state O2 concentrations based on Figure 1a. However, Figure 1a\nonly plots the median solution the stochastic analysis. Do all five the steady-states always exist for\nall sampled parameter combinations? Or are there sometimes only three steady-states (or only one\nsteady-state)? It would be useful if the article elaborated on this, and made it clear that all five\nsteady-states only exist for a subset for the possible input parameters (if this is the case).\nAn additional comment is that stable steady State I is not compatible with the Archean geologic\nrecord and photochemical modeling. State I occurs between 1e-7 and 1e-5 PAL O2. This range of\nO2 levels is about the same O2 concentration where S-MIF ceases to occur (~5e-8 to ~1e-6 PAL\naccording to photochemical modeling in Wogan et al. (2022), PNAS). Therefore, if Earth occupied\nState I during the Archean, we would expect S-MIF to occasionally turn on and off in the geologic\nrecord throughout the eon. However, we do not observe oscillating S-MIF during much of the\nArchean. I recommend pointing out this inconsistency in the text, or alternatively consider\nadjustments to the model so that lower O2 concentrations are allowed for O2. For example, if State\nI occurred between 1e-9 and 1e-5 PAL O2 in your stochastic analysis, then the lower range would\nbe compatible with the Archean geologic record.\nPlease also see the detailed line-by-line comments below. Perhaps some of them can improve the\nmanuscript.\n## Detailed comments\n- Line 31: The sentence starting with \"After the GOE's...\" is awkwardly worded and long. Maybe it\nshould be two sentences?\n- Line 40: Missing closed parentheses.\n- Line 41: Replace \"oxygenation\" with \"oxygen\". Replace \"was\" with \"has been\".\n- Line 60: Consider rewording \"The dynamical system analysis in its simplest definition is the study\"\nto \"Dynamical system analysis, in its simplest definition, is the study\"\n- Line 70: Why doesn't Equation (1) include sinks of oxygen from volcanic outgassing (i.e. outgassed\nH2, CO, etc.)?\n- Line 72: Include the units of pO2 (is it a partial pressure in bars?). Include units for terms like\nJ_NPP\n- Line 89: Typo. Replace \"has\" with \"have\"\n- Figure 1a: consider including the 68% and 98% confidence intervals on this plot. A question that\njumps to mind is that do all combinations of parameters have the 5 listed steady states? Or are\nthere some parameter sets (which are plausible) that tell a very different story compared to the\nmedian line given in the figure.\n- Line 155: The paper states \"Such a state, however, unlike the first state, is unstable, allowing for\nthe pO2 level to oscillate around this point in response to changes in the strength of global oxygen\nsinks or sources.\" This would not be my interpretation of unstable steady-state II. I would expect\npO2 to oscillate around a stable steady-state, but not an unstable steady state. For a stable steady-\nstate (e.g. State I), natural perturbations in the O2 flux to the atmosphere will cause O2 levels to\nincrease and decrease (i.e. oscillate) but ultimately be attracted back to stable steady state. The\nunstable State II is a tipping point which I would expect the atmosphere to occupy only transiently.\n- Line 159: The paper states \"Therefore, fluctuations in the atmospheric oxygen level around State II\ncould have resulted in an oscillation in the production of MIF-S signals as well—consistent with\nrecent proxy records.\" I do not like this explanation for reasons given in my comment to Line 155. A\npossible alternative explanation is that stable steady-state I shifted to higher values (for whatever\nreason), which allowed for S-MIF to turn on and off as O2 oscillated about the stable steady state. The\nrange of State I is between 1e-7 and 1e-5 PAL. S-MIF does not occur for O2 levels bigger than\nabout ~5e-7 PAL (there is uncertainty in this value). For example, perhaps State I was at 1e-7 for a\nlong time in the Archean. At the end of the Archean, for some reason (e.g. cyanobacteria emerging),\nState I shifted to be at 5e-7 PAL. At this point oscillations in S-MIF are expected.\n- Line 159: A further additional comment is that State I seems to be fairly centered on the S-MIF\ntransition, which might be unrealistic. Photochemical model suggest that S-MIF will stop above\n~5e-8 to ~1e-6 PAL (e.g. See Figure 1 in Wogan et al. (2022), PNAS). State I is between 1e-7 and 1e-5\nPAL. Therefore, if indeed stable O2 concentrations are between 1e-7 and 1e-5 PAL, then I would\nexpect natural perturbations/oscillations to Earth's redox to turn on and off S-MIF during much of\nthe Archean. However, we know that this didn't happen from the geologic record. In summary, the\nbig overlap between State I and the O2 levels where S-MIF ceases to occur seems inconsistent with\nthe persistent S-MIF in the Archean.\n- Line 159: Finally, the stable State I (between 1e-7 and 1e-5 PAL) does not agree atmospheric\nstability derived from photochemical modeling. For example, photochemical modeling in Gregory\net al. (2021), EPSL, suggests that O2 is unstable between 3e-6 and 1.0e-2 PAL. Furthermore,\nphotochemical modeling in Wogan et al. (2022), PNAS argues for instability between 5e-8 and 5e-4\nPAL. This inconsistency is OK and does not need to be explained or justified, however I recommend\nstating this inconsistency in the Discussion section of the article.\n- Line 202: This is a very long sentence. I suggest splitting it into two sentence for clarity.\n- Line 204: Typo. Earth does not have a upper case \"e\"\n- Line 254: Include units for pO2 here.\n- Line 291: The sentence starting with \"The fraction of oxic...\" is awkward\n- Line 322: T is temperature in what units? What are the units of pCO2?\n- Line 332: Typo. \"atom\" should be plural.\n- Line 356: Cite the papers where you acquired this parameterization for CH4 photochemistry.\n- Line 359: This parameterization (wherever it is derived) probably assumes a particular age of the\nSun (e.g. 2.0 Ga). This parameterization also probably assumes a particular climate of Earth (e.g.\neffect CH4 and O2 photochemistry. The age of the Sun determines the Lyman alpha UV flux, which\nis important for photolysis of CH4, which ultimately catalyzes the mutual destruction of CH4 and\nO2 (Maybe see Claire et al. (2012), ApJ for how UV flux changes with stellar age). Climate is\nimportant too because it controls the water vapor content of the troposphere (hotter means more\nH2O in troposphere). H2O photolysis makes OH. OH reacts with CH4, ultimately driving mutual\nannihilation of CH4 and O2 (see Ruiz et al. (2023), EPSL). It would be useful to state these caveats\nof the parameterization in the methods section.\n- Line 364: I'm confused at how you solve Equation (12). What are the knowns, and what are the\nunknowns? Below I write out the equation\n$$0 = \\frac{d(p_\\mathrm{CH_4})}{dt} = J_\\mathrm{CH_4} - J_\\mathrm{CH_4,O_2} - J_\\mathrm{Hydrogen}$$\n$$0 = \\left( \\alpha_1 (1-\\epsilon) (NPP) (f_\\mathrm{anoxic} (area)) \\right) - \\left( \\Psi (\\mu p_\\mathrm{O_2}) (\\mu p_\\mathrm{CH_4})^{0.7} \\right) - ( \\alpha \\mu p_\\mathrm{CH_4}) $$\nWhat are the knowns and unknowns in the above equation?\n- Line 382: Seems like there may be a typo in the sentence beginning with \"With this in mind...\"\n- Line 394: I recommend the authors publish the code for reproducing a subset of their result on\nZenodo. For example, upload a script that reproduces Figure 1a. The calculations should run readily\non a laptop in a few seconds at most.\n- Table S1: It seems strange to sample atmospheric O2 and CH4 levels independently. For example,\nO2 at 20% of the atmosphere is not photochemically compatible with methane at 3000 ppm (the\nrequired fluxes of both gases are way to big). A better model might do a sampling that treats O2 and\nCH4 as dependent on each other.\n- Table S1: The number of moles in an atmosphere that has 1 atm pressure at the surface, depends\non the mean molecular weight of the atmosphere. An atmosphere of pure CO2 needs less moles\ncompared to an atmosphere of pure H2 to make a surface pressure of 1 atm.\n- Table S1: Typo with \"Atmospehric\" in the last row.\n- Figure S2: The caption does not explain all the different shadings. What sort of confidence\nintervals are these?\n- Figure S3: What does the shading mean in these figures? The figure caption could contain this\ninformation.", "reviewer_2": "This paper proposes that three broad stable states exist for atmospheric O2 on Earth, as a\nconsequence of two unstable states that lie at around 1e-4 and 1e-1 PAL. The first of these unstable\nstates, which is implicated in the GOE, is already well documented by e.g. Goldblatt et al. (2006) in\nNature, and the present work uses approximately the same equations to represent it. The second\nunstable state at ~0.1 PAL is therefore the key focus of this work, and is an intriguing idea that\nwould help us better explain the dynamics of atmospheric oxygen during the Meso-Neoproterozoic\nand earlier Phanerozoic.\nI began this review with interest, and I wish I could be more positive about the work, but\nunfortunately there is a real lack of analysis in the paper which means that I cannot verify the\nconclusions.\nThere is no text about what drives the stable and unstable states in this model and whether they are\na reasonable approximation of the real system behavior. Even the instability around the Great\nOxidation, which has been covered in detail before, is not given any mechanistic explanation in this\npaper. Going back to the original work, it was driven by the effect of ozone shielding on the O2-CH4\nreaction, but this is not noted in the present work - although I can see the equation for the process.\nOverall the key question this work raises remains unanswered: what is the\nphysical/biogeochemical reason for the new instability that the authors have found around 0.1 PAL\n(state IV)? There is so much space in the nature comms format to get into this. It should be possible\nto find an answer by looking at which terms are driving the instability and by running alternative\nensembles with processes omitted.\nUltimately, it is hard to give a verdict on suitability of this work for publication, because I do not\nthink the work is finished. Maybe some simple analysis will provide a reasonable mechanistic\nexplanation for this interesting instability. But also maybe it won’t. There are some worrying kinks in\nfigure 1 at precisely 0.1 PAL and 0.3 PAL which hint at the model being prescribed certain shifts in\nbehavior at these points, right where the critical instability lies. This makes a realistic physical\nreasoning behind the model behavior seem unlikely to me.\nLine by line comments:\nTitle: I personally like a title to give some hint as to the conclusion of the paper. The current title\nfeels like a review paper to me.\nLine 31: GOE rather than GOE’s?\nLine 43: “although support for this rise has come more from biogeochemical modeling than proxy\nrecords”. There is proxy support here through Fe speciation (Sperling), Ce anomalies (Wallace), C\nisotope fractionation (Edwards) and fossil charcoal (Glasspool). In my opinion there is better proxy\nsupport for a Paleozoic stepwise shift in O2 than there is for a Neoprotrozoic one.\nLine 55: “with the traditional or commonly held view of largely invariant atmospheric oxygen levels\nduring the Proterozoic”. I agree, but is this view actually commonly held anymore? The authors\nthemselves give several high profile examples of more dynamic O2. Perhaps this can be rephrased\nmore confidently?\nThe section describing the model does not introduce what is actually in the model and on what\nmodels it builds. Having equations in the methods is fine but we need to know what you have done\nhere in simple terms – much of this is taken from previous work and can be easily explained as\nsuch.\nWhat are the wide range of environmental and initial conditions tested? There is no obvious list or\ntable here? And no discussion at all in the main text.\nLine 89: have instead of has?\nLine 100: Earth system states instead of Earth’s system states\nLine 149: “the increased”\nLine 154: Earth system rather than Earth’s system\nLine 157: The minimum\nLine 171: during the Proterozoic\nLine 173: what is a ‘fully stable state’ as opposed to a stable state?\nLine 179: “A higher oxygen value”?\nIs figure 1 produced with plants or without? Surely there should be two lines here for the plants and\nno plants system? Ah ok it is in the SI. This should probably be brought into the main figure, it’s\nimportant!\nMore on figure 1. Many other things have changed over earth history that should change this curve.\nE.g. emergence of continents, reducing CO2 levels, increasing solar flux, changing mantle inputs.\nHow do you take account of all this?\nLine 203: what is a “fully flat stable state”?\nFinal paragraph “According to our results, the atmospheric oxygen level could have fluctuated\ngreatly for most of the Earth’s history until the transition into the last stable state in the late\nProterozoic and early Phanerozoic, where the oxygen level finally reached its modern value”. This is\nconsistent with your model but so are a great many other hypotheses. Indeed the abstract seems at\nodds here as it favors the low background with a few notable excursions.", "reviewer_3": "Review of ‘Dynamics of Atmospheric Oxygenation’ by Fakhraee & Planavsky\nOverview/key results\nThis study utilizes a dynamical analysis approach to constrain O levels over Earth history. The\n2\nresults show three stable steady state solutions, at pO values of ~10-6 PAL, ~10-3 PAL, and ~1\nconstrain oxygen levels over Earth history, comparing them to the values summarized in an\nexisting review paper (see Figure 3 in particular). Particularly striking implications of this part of\nthe study are modern-like O values coincident with the Lomagundi-Jatuli carbon isotope\n2\nexcursion (~2.3—2.0 Ga) and low pO during the mid-Proterozoic (though the former is not\n2\nhighlighted in the study). The manuscript claims that the results indicate that pO could have\n2\nfluctuated during the mid-Proterozoic between 10-4 and 10-2 PAL.\nThis is a really interesting study, which applies a powerful approach to an important research\narea.\nReviewer’s summary using Nature Communications review criteria:\nPlease see ‘Comments/suggestions’ section for further details.\nValidity\nFurther clarity and detail are required to confirm the robustness of the data interpretation. For\nexample, Figure 1b contains important information necessary to judge the validity of the work\nbut is not adequately described.\nSignificance\nWith sufficient justification, the results could be of significance to the field. Constraining pO\n2\nvalues, particularly for the Mesoproterozoic, is an area of high interest and debate, in which\nvarious kinds of modeling and geochemical studies can come to considerably different\nconclusions. However, I think that such implications could be justified further and discussed\nfurther in the text. For example, Figure 3 suggests that oxygen levels were roughly modern\nduring the Lomagundi-Jatuli event, but this is not really mentioned in the text, despite being a\nstriking implication of the results, and it is insufficiently justified. I have noted some other\nsimilar suggestions in the ‘Comments/suggestions’ section, including small corrections to better\ntake into account previous work. In addition, the differences between this study and previous\nstudies should be explained more clearly, because the significance of this study as a scientific\nadvance in the context of other work is not communicated clearly.\nData, methodology, and analytical approach\nThe dynamical analysis appears to be a good approach to understand oxygen levels over Earth\nhistory. Modeling of atmospheric steady states can be an important tool to understand this\nproblem. The approach takes into account relevant sources and sinks for atmospheric oxygen,\nincluding atmospheric oxidation of methane. The ‘Methods’ section is detailed and contains\nmost of the important information. There are a few small details that could be included, and I\nhave made a note of these in the ‘Comments/suggestions’ section below. The figures are\nrelevant and well-chosen, but some important information required to understand a couple of\nthe figures fully is required, and some figures contain small errors that could be misleading (see\nbelow). The main results/output are shown in Figure 1b, but (as mentioned above) there is an\ninsufficient description of this information, which makes it difficult to assess how well the\nresults support the conclusions. Figure 1b is particularly important, as I think it displays the\nspread in possible steady state solutions and therefore gives an indication of the uncertainty,\nbut it is insufficiently explained.\nClarity, context, and references\nThe article is concisely written and very logically organized. The abstract is clear and accessible.\nThe content is pitched at a level suitable for the broad readership of Nature Communications.\nThe manuscript discusses some previous modeling work that also finds stable and unstable\nsteady state solutions and places this work in that context.\nPlease find more detailed comments and suggestions below.\nComments/suggestions\nMain text:\nThere are some grammatical errors throughout that occasionally obscure the clarity of the\narguments.\nL54-56: “it is becoming increasingly difficult to reconcile insights from the rock record with the\ntraditional or commonly held view of largely invariant atmospheric oxygen levels during the\nProterozoic.” I suggest that this sentence is reworded, as it is not clear whether references 17\nand 24 are given to support the “it is becoming increasingly difficult to reconcile insights” claim\nor given as references for the “traditional or commonly held view.” If the latter, the view of\nhigher oxygen levels during the Proterozoic also counts as “insights from the rock record” so I\nthink this sentence should be rephrased to acknowledge this. In other words, there is evidence\nfor higher oxygen levels as well as lower oxygen levels; it isn’t just a matter of evidence in the\nrock record vs. the “traditional” view.\nL96: “Our results were fairly stable… choosing each modeling input parameter.” What is meant\nby “fairly stable,” and why does this mean that “the simulation number is sufficient”? Also, the\nnumber of models run should be stated more clearly.\nL113: “According to our results, the only stable states occur under States I, II, and V.” There is a\nmistake in this sentence---it should be “under States I, III, and V.” Also, it would be helpful to\npoint out that stable steady states are defined as those when d2(pO )/dt2 < 0 and that unstable\n2\nsteady states are defined as those when d2(pO )/dt2 > 0 (perhaps using equivalent qualitative\n2\nlanguage instead of these equations).\nL128: “While our results are in accordance with these modeling results, they further suggest a\npossibility of an unstable state around 10-4 PAL which allows for a possible fluctuation in the\natmospheric oxygen level around this level.” The manuscript includes comparison to Goldblatt\net al. (2006), which is very important. However, I think this sentence should be clarified,\nbecause it is not clear to me how the results here differ from Goldblatt’s finding of two stable\nsteady states with an unstable steady state solution between. Sorry if I’m just missing\nsomething.\nL185: “This conclusion is in accordance with geological evidence (e.g., sulfur isotope\nsystematics, Cr isotope, among other redox proxies.” I suggest this is reworded as “in\naccordance with some geological evidence” or equivalent, as some geological evidence\nsuggests pO above 1% PAL.\nexplained. The parameter changed between these tests should be mentioned more explicitly\nhere and/or in the Methods.\nL203: “possible fully flat stable state for oxygen levels higher than modern value.” The meaning\nof “fully flat stable state” should be clarified. Also, looking at Figure S2 (left hand), it looks like\nthe d(pO )/dt is negative at pO > 1 PAL, rather than zero, so I don’t understand why this is\n2 2\nconsidered a stable state.\nL210: “… but these results suggest plants may also help limit large swings in atmospheric\noxygen levels.” This is not a novel result, so this clause should be slightly reworded, as it\ncurrently sounds like this study was the first to discover this.\nMethods:\nL255: J is repeated. One should be J .\nocean,sink CH4,O2\nL297: “To account for uncertainty… for the range of each parameter.” This is a really helpful\nsentence. If the uniform distribution applies to all parameters in Table S1, rather than just the\nparameters related to J , then this sentence could be put further towards the top of the\nNPP,ocean\nMethods section, rather than here where it is more hidden.\nL322: A sentence could be added to this section to explain how values for T are obtained. Does\nT change, or is a fixed value chosen?\nL322: pCO is defined twice, and pCO is not defined.\n2 2,1\nFigures:\nFigure 1: The caption does not sufficiently explain Figure 1b. For example, the meaning of the\nhorizontal width of the shapes is not explained. The vertical black lines and the thin black lines\nat the edges of the shapes are not defined. I assume that the white circles show the median\nvalues, but this should be stated explicitly.\nFigure 2: The caption states that the vertical axis label “Energy” is the potential energy, but as a\nreader I am unclear about what that means. I do not understand why the different stable and\nunstable states have different potential energies nor how that is calculated. I suggest that more\ndetail is given here, both in the caption and in the text.\nFigure 3b: This is a summary/implications figure, so it is important that it accurately reflects the\nresults.\n• The caption does not sufficiently explain the meaning of the blue band and the black-\noutlined box – e.g., do they encapsulate all the models run, or 95% of them, etc.? A really\nhelpful thing to note would be the relationship of the width of these bands to the shapes\nand/or black vertical lines on Figure 1b, as it is difficult to work out where the width of the\nbands come from.\n• Relatedly, it is important that the components of Figure 3b are drawn accurately. For\nexample, the lower bound of the mid-Proterozoic band (this study) is ~8e-5 PAL, but the\nlower band of the State III shape is 2e-4, and the lower limit of the black vertical line is even\nhigher. Though the difference between 8e-5 and 2e-4 PAL is not huge, it could be a little\nmisleading as drawn, as the mid-Proterozoic box appears to suggest slightly lower pO\n2\nvalues than the results support.\n• Very minor: Similarly, the medians and black-outlined boxes of the 2nd (Lomagundi) and 4th\n(Phanerozoic) boxes should be at exactly the same pO , but a careful examination shows\n2\nthat they are a little offset. This is extremely minor, but should be tidied up before\npublication. The lower limit of the mid-Proterozoic curve from Lyons et al. (2014) is 1e-4\nPAL, while it is drawn slightly higher here.\n• The figure as drawn suggests that oxygen levels were ~modern during the time of the\nLomagundi excursion, shown in the figure between 2.3 and 2 Ga. However, the State III\nconstraints (pO ~1e-3 PAL) could have alternatively been used for this period of Earth\n2\nhistory. I assume that they were not used in the figure because of evidence mentioned in\nthe text of pO > 1% PAL (lines 33-34), but since the blue band lies outside the Lyons et al.\n2\n(2014) constraints box anyway, it is not clear why the possibility of the State III levels during\nthis time is not included. There should be some discussion of this particular period of Earth\nhistory in the text. If leaving Figure 3 as it is to show only the State V levels, the discussion\nshould comprise either: (1) a justification of ~modern pO levels for this period, explicitly\n2\nstating why pO < 1% can definitely be dismissed, as it would be a really important\n2\nimplication of the current work if this can really be robustly defended, or (2) the limitations\nand caveats of the results’ support of State V oxygen levels only for this time.\nSupplementary Information:\nTable S1: “Coefficient for ocean anoxic sink.” If there is no reference for the range of numbers\nchosen for k (0.01—0.1), where did this range come from? I see that in L373 in the Methods, it\n1\nis stated that it is varied “within a wide range” to account for uncertainty, but I think some\njustification for this particular range should be included.\nFigure S3: The plots presented here contain a lot of useful information, but I can’t see any\ndiscussion of what they show, either in the main text or the SI, and the figure is not referred to in\nthe text. It would be beneficial to add some sentences describing the main takeaways from\nthese plots, because there is some interesting information there. In addition, a longer caption\nwould be helpful.\nThank you very much for an interesting paper to review. I hope to see it published in due\ncourse.", "reviewer_author_response": "Response to Reviewers\nReviewers’ comments are in black.\nOur responses are in blue.\n-------------------------------------------------------------------------------------------------------------------------------\nReviewer #1 (Remarks to the Author):\nSummary\nThe paper constructs a dynamic stochastic model of the oxygen concentration on Earth. Their model is a\nsimple ordinary differential equation describing the rate of change of O levels (dpO2 /dt). For a wide range\n2\nof input parameters, the find the steady-state O2 concentrations between 1e-13 and 1 present atmospheric\nlevels (PAL). The nominal, median model find 5 steady-states. Three of the steady-states are stable, while\n2 are not stable. The authors then argue that the history of oxygen on Earth might be understood in terms\nof these steady-states. One stable steady state occurs for Archean amounts of O2 (~3e-6 PAL), Another\nstable steady state occurs at ~1e-3 PAL O2, perhaps representative of the Proterozoic eon. Finally, the\nmodern Earth is at the final stable steady state at ~ 1 PAL O2.\nOverall comments\nOverall, the paper is fairly well written, concise, and generally makes sense. It is a valuable addition to the\nliterature because considers a new stochastic dynamic approach to determine what O2 concentrations are\nmost likely to be stable. This is a new contribution to the literature, based on my knowledge. I have a few\nrelatively minor comments that I would like to see addressed before publication.\nThanks for the time and the feedback on our work.\nThe article argues for five steady-state O2 concentrations based on Figure 1a. However, Figure 1a only\nplots the median solution the stochastic analysis. Do all five the steady-states always exist for all sampled\nparameter combinations? Or are there sometimes only three steady-states (or only one steady-state)? It\nwould be useful if the article elaborated on this, and made it clear that all five steady-states only exist for a\nsubset for the possible input parameters (if this is the case).\nThanks for this point. The results from our stochastic simulation were filtered so that they would be\nconsistent with the evidence of a steady state at a high oxygen level which is the near-modern value, and a\nsteady state under widely anoxic conditions (pO < 10-5 present atmospheric level (PAL)), consistent with\n2\nthe geochemical proxies (notably mass-independent fractionation of sulfur MIF-S) and a steady state\nbetween these two end members. For instance, we did not consider the results from the stochastic simulation\nthat show no steady state under anoxic conditions or no steady state under high modern-like oxygen levels\nas they both would be inconsistent with accepted evidence from the geochemical proxies. Nonetheless,\nmost of the results from the stochastic simulation do show an emergence of the five steady states two of\nwhich would be unstable.\nWe have now clarified this point in the text:\n“Our results from the stochastic simulation were filtered to be consistent with the geochemical proxies.\nSpecifically, we only considered the results from the stochastic simulation that include a steady state point\nunder high oxygen modern-like conditions, a steady state point under anoxic conditions (pO < 10-5 present\n2\natmospheric level (PAL)) which is consistent with the presence of anoxic atmosphere during Earth’s early\nstage (e.g., mass-independent fractionation of sulfur), and a steady-state point between the two end\n1\nmembers of high and anoxic conditions that would generally correspond to the geologic times where there\nis no evidence for anoxic nor high oxygen conditions.”\nAn additional comment is that stable steady State I is not compatible with the Archean geologic record and\nphotochemical modeling. State I occurs between 1e-7 and 1e-5 PAL O2. This range of O2 levels is about\nthe same O2 concentration where S-MIF ceases to occur (~5e-8 to ~1e-6 PAL according to photochemical\nmodeling in Wogan et al. (2022), PNAS). Therefore, if Earth occupied State I during the Archean, we would\nexpect S-MIF to occasionally turn on and off in the geologic record throughout the eon. However, we do\nnot observe oscillating S-MIF during much of the Archean. I recommend pointing out this inconsistency in\nthe text, or alternatively consider adjustments to the model so that lower O2 concentrations are allowed for\nO2. For example, if State I occurred between 1e-9 and 1e-5 PAL O2 in your stochastic analysis, then the\nlower range would be compatible with the Archean geologic record.\nThanks for this point. Following your suggestion, we have conducted another stochastic simulation with a\nwider range of parameters for hydrogen escape flux which resulted in a shift in the first steady state point.\nThe updated results are copied below. The range of value for the first steady state is now between ~8x10-9\n(25th percentile) – 2x10-7 (75thpercentile) with a median of ~ 3x10-8 and minimum and maximum values of\n~1x10-13 and 4x10-7 PAL.\nUpdated figure 1. The results of stochastic dynamical analysis of global oxygen mass balance.\nmanuscript.\nDetailed comments\n• Line 31: The sentence starting with “After the GOE’s. . . ” is awkwardly worded and long. Maybe it should\nbe two sentences?\nThanks. The text is now modified as:\n“After the GOE, atmospheric oxygen levels are generally accepted to have passed the threshold (~ 10-6 -\n10-5 present atmospheric level (PAL); 4) for the formation of mass-independent fractionation of sulfur (MIF-\nS), as a geochemical proxy for an anoxic atmosphere 5,6. And almost certainly reached values >1% PAL or\nhigher during the geological interval known as the Lomagundi-Jatuli Event that records the highest\nexcursion in the isotopic composition of carbon in marine carbonate (δ13C)1,3,7,8.”\n• Line 40: Missing closed parentheses.\nCorrected. Thanks.\n• Line 41: Replace “oxygenation” with “oxygen”. Replace “was” with “has been”.\nCorrected. Thanks.\n• Line 60: Consider rewording “The dynamical system analysis in its simplest definition is the study” to\n“Dynamical system analysis, in its simplest definition, is the study”\nModified as suggested. Thanks.\n• Line 70: Why doesn’t Equation (1) include sinks of oxygen from volcanic outgassing (i.e. outgassed H2,\nCO, etc.)?\nThanks for this point. The sink from volcanic outgassing was incorporated into the rate of oxidative\nweathering as well as the hydrogen escape in the stochastic simulation. We have now explicitly included\nthis term in the equation. Since this effect was already accounted for in the model, this change did not affect\nour results nor the conclusion presented in the main text.\nThe corrected equation is now:\n“\n𝑑(𝑝𝑂2)\n= 𝐽 +𝐽 −𝐽 − 𝐽 − 𝐽 − 𝐽 − 𝐽 (1).”\n𝑁𝑃𝑃 𝐻𝑦𝑑𝑟𝑜𝑔𝑒𝑛 𝐶𝐻4,𝑂2 𝑣𝑜𝑙𝑐 𝑜𝑐𝑒𝑎𝑛,𝑠𝑖𝑛𝑘 𝑜𝑟𝑔,𝑤𝑒𝑎𝑡ℎ 𝑝𝑦𝑟,𝑤𝑒𝑎𝑡ℎ\n𝑑𝑡\n• Line 72: Include the units of pO2 (is it a partial pressure in bars?). Include units for terms like J_NPP\nThanks, we included a unit for the J_NPP flux. Also, since we solve for the dpO2/dt not pO2, the unit for\npO may not be directly relevant. We have now included a unit for dpO /dt as well.\n2 2\n• Line 89: Typo. Replace “has” with “have”\nCorrected. Thanks.\n3\n• Figure 1a: consider including the 68% and 98% confidence intervals on this plot. A question that jumps\nto mind is that do all combinations of parameters have the 5 listed steady states? Or are there some 1\nparameter sets (which are plausible) that tell a very different story compared to the median line given in the\nfigure.\nThanks for this comment. We have updated Figure 1 to include information on the 25 and 75 percentile\nrange of values for each steady state point as well as the minimum and maximum value for each steady\nstate. We have also included this information in the main text. The updated text and the updated Figure 1\nnow appear as:\nUpdated text:\n“The results from the stochastic analysis of the global oxygen mass balance suggest that atmospheric\noxygen can be in balance under essentially anoxic conditions (state I; ~8x10-9 (25th percentile) – 2x10-7\n(75thpercentile) with median ~ 3x10-8 and minimum and maximum values of ~1x10-13 and 4x10-7 PAL),\nvery low oxygen conditions (state II; 1x10-5 – 1x10-4 with median ~ 3x10-5 and minimum and maximum values of\n3x10-7 and 2x10-4 PAL), low oxygen condition (state III: ~ 1x10-3 – 1x10-2 with median ~ 3x10-3 and\nminimum and maximum values of ~8x10-5 and 3x10-2 PAL), relatively high oxygen condition (state IV: ~\n1x10-1 – 4x10-1 with median ~ 2x10-1 and minimum and maximum values of ~10-2 and 7x10-1 PAL), and high,\nnear modern oxygen conditions (state V: 9x10-1 – 1.1x100 with median ~ 100 and minimum and maximum\nvalues of 7x10-1 and 1.3x100 PAL) (Figure 1b).”\nUpdated Figure:\nUpdated figure 1. The results of stochastic dynamical analysis of global oxygen mass balance.\nConcerning your second point on the presence of 5 steady states, as we explained above, we only considered\nthe results of the stochastic simulation that are consistent with the geochemical proxies that point toward at\nleast one steady state point under anoxic condition (e.g., MIF-S signal), a modern-like steady state point,\nand a steady state point between these two end-members that are consistent with the absence of evidence\nfor anoxic (e.g., MIF-S signal) and high oxygen conditions during Earth’s middle age.\n4\n• Line 155: The paper states “Such a state, however, unlike the first state, is unstable, allowing for the pO2\nlevel to oscillate around this point in response to changes in the strength of global oxygen sinks or sources.”\nThis would not be my interpretation of unstable steady-state II. I would expect pO2 to oscillate around a\nstable steady-state, but not an unstable steady state. For a stable steady-state (e.g. State I), natural\nperturbations in the O2 flux to the atmosphere will cause O2 levels to increase and decrease (i.e. oscillate)\nbut ultimately be attracted back to stable steady state. The unstable State II is a tipping point which I would\nexpect the atmosphere to occupy only transiently.\nThanks. Yes, you are right. We have now rephrased the text as follows:\n“Such a state, however, unlike the first state, is unstable, and the system would only occupy such state\ntransiently (e.g., tipping point) and a small change in the rate of oxygen production or consumption would\nresult in a transition into the adjacent stable states (state I and state III).”\n• Line 159: The paper states “Therefore, fluctuations in the atmospheric oxygen level around State II\ncould have resulted in an oscillation in the production of MIF-S signals as well—consistent with recent\nproxy records.” I do not like this explanation for reasons given in my comment to Line 155. A possible\nalternative explanation is that stable steady-state I shifted to higher values (for whatever reason), which\nallowed for S-MIF to turn on and off as O2 oscillated about the stable steady state. The range of State I is\nbetween 1e-7 and 1e-5 PAL. S-MIF does not occur for O2 levels bigger than about ~5e-7 PAL (there is\nuncertainty in this value). For example, perhaps State I was at 1e-7 for a long time in the Archean. At the\nend of the Archean, for some reason (e.g. cyanobacteria emerging), State I shifted to be at 5e-7 PAL.\nAt this point oscillations in S-MIF are expected.\nThanks for this point. We agree that the oscillation in the MIF-S signal can be explained by a fluctuation\nwithin a stable state. We would like to, however, note that while the minimum value for the threshold of\nMIF-S formation (~10-7PAL; Wogan et al. 2022) is consistent with the maximum values of atmospheric\noxygen under state I (~2x10-7PAL), the range required for the oscillation around the threshold of MIF-S\nformation (~10-6 PAL) lies within the state II. Our range of atmospheric oxygen at the unstable state of II is\nconsistent with the range of atmospheric oxygen suggested by the previous studies which characterized it\nas an unstable state (Wogan et al. 2022; Goldblatt et al. 2006). Oscillation in the MIF-S production in a\nrelatively short geological timescale can alternatively be explained by a fluctuation in a state in which\natmospheric oxygen would be most sensitive to a small change in the rate of oxygen consumption and\nproduction and such condition can only occur at an unstable state. Nonetheless, we have discussed both\npossibilities for the oscillation of MIF-S signal which involves the fluctuation in the atmospheric oxygen\nlevel in a stable state as well as the alternative explanation that involves an oscillation around an unstable\nstate II.\nThe following text is now added to the main text:\n“According to our results, an oscillation in the MIF-S production, suggested by the recent results from\ngeochemical proxies 14,34, can be explained by a fluctuation in the atmospheric oxygen level in both states\nI and II. Specifically, the maximum value in state I is ~2x10-7PAL which is higher than the minimum value\nsuggested for the threshold of atmospheric oxygen for caseation in the production of MIF-S (e.g., 23) which\nimplies that it is possible that a fluctuation in the atmospheric oxygen level under state I result in an\noscillation in the production of MIF-S signal. Suggested higher values for the threshold of atmospheric\noxygen level for MIF-S formation (~10-6 PAL) are also consistent with the range of atmospheric oxygen\nlevel at state II which suggests that the apparent fluctuation in MIF-S signal can alternatively be attributed\nto an oscillation in the atmospheric oxygen level at unstable state II. Regardless, according to this\n5\nframework, the GOE can then be defined as a departure from unstable State II and the transition to the\nsecond stable state for atmospheric oxygen level (state III: Figure 1).”\n• Line 159: A further additional comment is that State I seems to be fairly centered on the S-MIF transition,\nwhich might be unrealistic. Photochemical model suggest that S-MIF will stop above ~5e-8 to ~1e-6 PAL\n(e.g. See Figure 1 in Wogan et al. (2022), PNAS). State I is between 1e-7 and 1e-5 PAL. Therefore, if\nindeed stable O2 concentrations are between 1e-7 and 1e-5 PAL, then I would expect natural\nperturbations/oscillations to Earth’s redox to turn on and off S-MIF during much of the Archean. However,\nwe know that this didn’t happen from the geologic record. In summary, the big overlap between State I and\nthe O2 levels where S-MIF ceases to occur seems inconsistent with the persistent S-MIF in the Archean.\nThanks for this comment. As explained above, we extended the range of parameters in our stochastic\nsimulation so that the range of atmospheric oxygen in the state I under anoxic conditions would be\nconsistent with the results from geochemical proxies.\n• Line 159: Finally, the stable State I (between 1e-7 and 1e-5 PAL) does not agree atmospheric stability\nderived from photochemical modeling. For example, photochemical modeling in Gregory et al. (2021),\nEPSL, suggests that O2 is unstable between 3e-6 and 1.0e-2 PAL. Furthermore, photochemical modeling\nin Wogan et al. (2022), PNAS argues for instability between 5e-8 and 5e-4 PAL. This inconsistency is OK\nand does not need to be explained or justified, however I recommend stating this inconsistency in the\nDiscussion section of the article.\nThanks for pointing this out. Our updated results from the stochastic simulation use a wider range of\nparameter space does exhibit a consistency between our results and the previous modeling works. For\ninstance, the range of atmospheric oxygen level at unstable state II (state II; 1x10-5 – 1x10-4 with median ~\n3x10-5 and minimum and maximum values of 3x10-7 and 2x10-4 PAL) is consistent with the range of unstable\nstate suggested by Wogan et al. (2022).\nWe have now discussed this point in the text:\n“Our results for the range of atmospheric oxygen at unstable state II (1x10-5 (25th percentile) – 1x10-4 (75th\npercentile) with median ~ 3x10-5 and minimum and maximum values of 3x10-7 and 2x10-4 PAL) are largely\nconsistent with the previous modeling results. Notably, the results from the photochemical modeling suggest\ninstability of atmospheric oxygen between 5x10-8 and 5x10-4PAL that are largely consistent with our results\nfor the range of atmospheric oxygen level at unstable state II 23.”\n• Line 202: This is a very long sentence. I suggest splitting it into two sentence for clarity.\nThanks. We have now revised the text for this section. The revised text now reads as:\n“The evolution of plants in the Phanerozoic likely helped stabilize atmospheric oxygen around modern\nvalue. To explore the impact of plant evolution on the atmospheric oxygen level, we conducted two\nstochastic simulations with and without plants. We explored the effect of the rise of plants by removing the\nflux of carbon burial linked to land plants from the total oxygen production flux (J equation 2 in the\nNPP;\nmethods section). According to our results, the median atmospheric oxygen level without plants would be\naround 0.4 atmospheric oxygen level (PAL). This result shows that the rise of plants resulted in an increase\nin the oxygen production flux and helped with stabilizing the atmospheric oxygen level near modern value\n(Figure 4). Similar to the results from previous studies (e.g., 13), our results also suggest that, by exerting\nnegative feedback such as fire, plants may also help limit large swings in atmospheric oxygen levels.”\n6\n• Line 204: Typo. Earth does not have a upper case “e”\nCorrected. Thanks.\n• Line 254: Include units for pO2 here.\nThe unit for pO /dt is included. Sorry. Thanks.\n2\n• Line 291: The sentence starting with “The fraction of oxic. . . ” is awkward\nThanks. The sentence is now modified as:\n“The fraction of oxic condition is calculated as one minus the fraction of anoxia (1- f ).”\nanoxic\n• Line 322: T is temperature in what units? What are the units of pCO2?\nThanks. The text is now modified as:\n“Where T is the temperature in the unit of degrees Celsius, pCO is CO atmospheric level in ppm”.\n2 2\n• Line 332: Typo. “atom” should be plural.\nCorrected. Thanks.\n• Line 356: Cite the papers where you acquired this parameterization for CH4 photochemistry.\nDone. Thanks.\n• Line 359: This parameterization (wherever it is derived) probably assumes a particular age of the Sun (e.g.\n2.0 Ga). This parameterization also probably assumes a particular climate of Earth (e.g. 280 K surface\ntemperature). Yet both the age of the Sun and the climate of Earth should significantly effect CH4 and O2\nphotochemistry. The age of the Sun determines the Lyman alpha UV flux, which is important for photolysis\nof CH4, which ultimately catalyzes the mutual destruction of CH4 and O2 (Maybe see Claire et al. (2012),\nApJ for how UV flux changes with stellar age). Climate is important too because it controls the water vapor\ncontent of the troposphere (hotter means more H2O in troposphere). H2O photolysis makes OH. OH reacts\nwith CH4, ultimately driving mutual annihilation of CH4 and 2 O2 (see Ruiz et al. (2023), EPSL). It would\nbe useful to state these caveats of the parameterization in the methods section.\nThanks. The following text is now added:\n“While our parameterization scheme for the rate of atmospheric methane oxidation is similar to previous studies,\nsome assumptions in this parameterization can potentially cause uncertainty in the rate of methane oxidation.\nFor instance, this parameterization makes assumptions about the age of the Sun as well as about the surface\ntemperature, both of which effects can potentially influence the photochemistry of CH and O . The assumption\n4 2\nabout the age of the sun can impact the Lyman alpha UV flux, an important factor for kinetics of CH photolysis\n4\non the abundance of OH and the rate of methane oxidation. To account for the uncertainty associated with\nthese assumptions we varied the kinetics of atmospheric methane oxidation in our stochastic simulation.”\n• Line 364: I’m confused at how you solve Equation (12). What are the knowns, and what are the unknowns?\nBelow I write out the equation\n7\n0 = d(pCH4 ) dt = JCH4 − JCH4,O2 − JHydrogen 0 = (α1(1 − ϵ)(NP P)(fanoxic(area))) − Ψ(µpO2 )(µpCH4\n) 0.7 − (αµpCH4 )\nWhat are the knowns and unknowns in the above equation?\nThanks for this point. In this equation, both pO and pCH are the only unknowns. By solving this equation,\n2 4\nthe methane steady-state points (fixed points) would be written as a function of pO and we then use this\n2\nrelationship to solve the change in the oxygen imbalance (dpO /dt; equation 1).\n2\n• Line 382: Seems like there may be a typo in the sentence beginning with “With this in mind. . . ”\nCorrected. Thanks.\n• Line 394: I recommend the authors publish the code for reproducing a subset of their result on Zenodo.\nFor example, upload a script that reproduces Figure 1a. The calculations should run readily on a laptop in\na few seconds at most.\nThanks for the suggestion. We have made the code available through a Zenodo link.\n• Table S1: It seems strange to sample atmospheric O2 and CH4 levels independently. For example, O2 at\n20% of the atmosphere is not photochemically compatible with methane at 3000 ppm (the required fluxes\nof both gases are way to big). A better model might do a sampling that treats O2 and CH4 as dependent on\neach other.\nThanks for this point. The methane concentration would be solved from the methane mass balance equation\nand the range of values in Table pertain to an original version of the work that did not include the methane\nmass balance. We have now omitted the range of methane in this table.\n• Table S1: The number of moles in an atmosphere that has 1 atm pressure at the surface, depends on the\nmean molecular weight of the atmosphere. An atmosphere of pure CO2 needs less moles compared to an\natmosphere of pure H2 to make a surface pressure of 1 atm.\nThanks for this point. We have now discussed this potential caveat in the method section.\nThe following text is now added to the methods:\n“It should be noted that the number of moles in 1 atm depends on the relative abundance of gases. For instance,\na CO -rich atmosphere contains less moles compared to an H -rich atmosphere. To account for this caveat\n2 2\nthe value of α was varied in our stochastic simulation.”\n• Table S1: Typo with “Atmospehric” in the last row.\nCorrected. Thanks.\n• Figure S2: The caption does not explain all the different shadings. What sort of confidence\nintervals are these? • Figure S3: What does the shading mean in these figures? The figure caption could contain this\ninformation.\nThanks. We have now included additional information for these figures.\n8\n-------------------------------------------------------------------------------------------------------------------------------\nReviewer #2 (Remarks to the Author):\nThis paper proposes that three broad stable states exist for atmospheric O2 on Earth, as a consequence of\ntwo unstable states that lie at around 1e-4 and 1e-1 PAL. The first of these unstable states, which is\nimplicated in the GOE, is already well documented by e.g. Goldblatt et al. (2006) in Nature, and the present\nwork uses approximately the same equations to represent it. The second unstable state at ~0.1 PAL is\ntherefore the key focus of this work, and is an intriguing idea that would help us better explain the dynamics\nof atmospheric oxygen during the Meso-Neoproterozoic and earlier Phanerozoic.\nI began this review with interest, and I wish I could be more positive about the work, but unfortunately\nthere is a real lack of analysis in the paper which means that I cannot verify the conclusions.\nThanks for the time and the feedback on our work. We have now conducted further analysis to address the\nconcerns. You may find our responses to the concerns below.\nThere is no text about what drives the stable and unstable states in this model and whether they are\na reasonable approximation of the real system behavior. Even the instability around the Great Oxidation,\nwhich has been covered in detail before, is not given any mechanistic explanation in this paper. Going back\nto the original work, it was driven by the effect of ozone shielding on the O2-CH4 reaction, but this is not\nnoted in the present work - although I can see the equation for the process. Overall the key question this\nwork raises remains unanswered: what is the physical/biogeochemical reason for the new instability that\nthe authors have found around 0.1 PAL (state IV)? There is so much space in the nature comms format to\nget into this. It should be possible to find an answer by looking at which terms are driving the instability\nand by running alternative ensembles with processes omitted.\nUltimately, it is hard to give a verdict on suitability of this work for publication, because I do not think the\nwork is finished. Maybe some simple analysis will provide a reasonable mechanistic explanation for this\ninteresting instability. But also maybe it won’t. There are some worrying kinks in figure 1 at precisely 0.1\nPAL and 0.3 PAL which hint at the model being prescribed certain shifts in behavior at these points, right\nwhere the critical instability lies. This makes a realistic physical reasoning behind the model behavior seem\nunlikely to me.\nThanks for this point. We followed your suggestion and have conducted additional analysis to address your\nconcern. Briefly, we have conducted an additional sensitivity analysis to explore the impact that each\nfeedback would have on the emergence of pO steady-state points (fixed points). The more recent sensitivity\n2\nanalysis of the pO steady-state points to the change in the strength of various oxygen feedbacks in the\n2\nsupplement. Specifically, we increased the strength of each feedback to its maximum value used in the\nstochastic simulation and tracked the change in the median of the dpO /dt at different pO values (initial\n2 2\nconditions). Our results from sensitivity analysis are shown in Figure S2. Our results indicate that the\nemergence of the first steady-state point (state I) under widely anoxic conditions (pO <10-5 present\n2\natmospheric level (PAL)) is regulated by the flux of hydrogen escape, consistent with the previous work\n(e.g., Goldblatt et al. 2006). Specifically, the first steady-state point is regulated by the balance between\noxygen production through the escape of hydrogen and oxygen consumption by reducing volcanic gases.\nAccording to our results, a change in the rate of atmospheric methane oxidation would influence the\nemergence of the second steady state point (state II) around 10-5 PAL, which is consistent with the modeling\nresults from previous studies. The steady-state points under low-oxygen conditions (state III and IV) are\nmainly controlled by the balance between oxygen production in the ocean and the strength of the oxygen\n9\nsink in the ocean (J ) and on land (J and J ). Specifically, while the rate of oxidative\nanoxia pyr,weath og,weath\nweathering is essentially ineffective in impacting the global oxygen mass balance under widely anoxic\nconditions (pO <10-5PAL), it would have a large impact on the emergence of steady-state points under\n2\nlow-oxygen conditions. Precisely, the first steady state point under low-oxygen conditions (state III) is the\nresult of the interplay between the rate of oxygen production in the ocean, caused by the rise of oxygenic\nphotosynthesis, as a source and the rate of oxidative weathering and the rate of reaction of oxygen with\nreducing species in the ocean as the sink of oxygen. At higher oxygen levels, the rate of oxidative\nweathering would increase which would have resulted in a shift of dpO /dt to the negative zone. According\n2\nto our result from the sensitivity analysis (Figure S2), an increased rate of oxidative weathering can result\nin the disappearance of state IV and state V. The emergence of state IV can then be explained as the result\nof an increase in the rate of oceanic oxygen production (J ) which would be accompanied by a decrease\nNPP\nin the rate of oxygen consumption (J ) and, in turn, would have outcompeted the rate of oxidative\nanoxia\nweathering on land. Given that this state is an unstable state and emerges after and prior to a stable state,\na\nsmall change in the ocean oxygen mass balance (Figure S2) would cause the system to move toward one\nof the stable states.\nAdditionally, the unstable state around 10-1 PAL (state IV) can be explained in the context of the dynamical\nanalysis. The presence of two stable states forces the system to have an unstable state that acts as a transition\nbetween two stable states (Figure R1). As shown in the schematic figure below (Figure R1) with two stable\nstates at high and low oxygen (solid circles), the transition between these two (the dashed line) is inevitable\nand would result in the emergence of an unstable state (open circle). This means that an unstable point at a\nlow pO level (~0.01 PAL) is the consequence of having a low pO stable and a high modern-like pO stable\n2 2 2\npoint and acts as a transition state from low-oxygen to high-oxygen conditions and its magnitude, as\nexplained above, would be controlled by the strength of sink of oxygen in the ocean and on land.\nRegarding your point on the model behavior around the unstable points, the model is not forced by any\nspecific parameter space and, as laid out in the model description in the main text as well as in the methods\nsection, we took a stochastic approach and no fixed value for any model parameter and the simulation were\nrun under a continuous random sampling of the model input parameters assuming a uniform distribution.\nNeedless to say the model main code and the code used to produce the results presented in Figure 1 is\navailable with the posted Zenodo link.\nWe have now included the newly made figure from the sensitivity analysis along with a discussion on the\ncontrols on various pO steady-state points in the text.\nMain text:\n“Our results from the sensitivity analysis suggest that the rate of hydrogen escape and the atmospheric\noxidation of methane have the largest influence on the emergence of the first and second steady-state for\npO under anoxic conditions (Figure S2), consistent with the previous study. Specifically, our results suggest\n2\nthat the first steady-state point is controlled by the balance between the sink of oxygen from reaction with\nreduced volcanic gases and the source of oxygen through the escape of the hydrogen to the atmosphere. An\nincrease in the rate of hydrogen escape would lead to the disappearance of the first steady-state point (state\nI)), indicating the key role of hydrogen escape in stabilizing the atmospheric oxygen level before the rise of\noxygenic photosynthesis. According to our results, an increase in the rate of atmospheric methane oxidation\nwould have influenced the emergence of the second steady-state point (state II) around 10-5 PAL, consistent\nwith the modeling results from previous studies.\nlow-oxygen conditions (states III & IV) is controlled by the oceanic oxygen production (J ) and the sink\nNPP\nof oxygen in the ocean and on land. An increased rate of oceanic oxygen production (J ) would shift the\nNPP\nimbalance value of oxygen (dpO /dt) to positive values which would influence the appearance of steady-\n2\nstate points under low-oxygen conditions. At the same time, according to our results, while under anoxic\nconditions (<10-5 PAL) the oxidative weathering of pyrite and organic matter is essentially ineffective in\ncontrolling the global oxygen mass balance (Figure S2 and S3), an increased rate of oxidative weathering\nunder low-oxygen condition would dramatically change the oxygen imbalance value which would impact\nthe emergence of steady-state points under low-oxygen conditions. Of the two steady-state points under\nlow-oxygen conditions (states III & IV), the first one is the result of an increase in oxygen production in the\nwake of the evolution of oxygenic photosynthesis while the rate of oxidative weathering is not high. An\nincrease in the atmospheric oxygen level would result in an increase in the rate of oxidative weathering\nwhich would result in a negative imbalance for the oxygen. A continued increase in the production of oxygen\nin the ocean would have resulted in a decrease in anoxia in the ocean which would have acted as positive\nfeedback in the global oxygen mass balance and have resulted in the emergence of the state IV. Such a state\nis, however, unstable which means that a change in the oceanic oxygen mass balance (e.g., a change in the\noxygen production through net primary production or marine oxygen consumption) would shift the system\nto the adjacent stable state (state III or state V), depending on the directionality of the change.”\nMethods:\n“To further explore the influence of each feedback on\nthe global oxygen mass balance, we conducted a\nsensitivity analysis. Specifically, to test the sensitivity\nof our results for the steady-state pO levels to the\n2\nstrength of each feedback, we fixed the value of each\noxygen feedback to its maximum value used in the\nstochastic analysis and recorded the change in the\nmedian of dpO /dt vs pO (Figure S2). Our results\n2 2\nsuggest that the rate of hydrogen escape and\natmospheric oxidation of methane exert the strongest\nleverage on the emergence of the first and second\nsteady-state points for pO . For steady-state points\n2\nunder low-oxygen conditions, the rates of oxygen\nconsumption in the ocean and on land (e.g., pyrite\nweathering) have the strongest influence. Notably,\nan increase in the net primary production would\nresult in the disappearance of the steady-state points\nunder low-oxygen and high-oxygen conditions. This\nresult suggests that a change in the marine primary\nproduction (e.g., nutrient delivery and or amount of\navailable oxygen-reducing species) can have a large\nimpact on the global oxygen mass balance.”\nFigure S2. The results from the sensitivity analysis of the pO2\nsteady-state point to the change in the strength of various\noxygen feedback.\nstable points. The stable points are shown in solid circles and the open circle denotes the unstable\npoint.\nLine by line comments:\nTitle: I personally like a title to give some hint as to the conclusion of the paper. The current title feels like\na review paper to me.\nThanks for the suggestion. We have now changed the title to:\n“Large dynamic in the geological atmospheric oxygen level”\nLine 31: GOE rather than GOE’s?\nCorrected. Thanks.\nLine 43: “although support for this rise has come more from biogeochemical modeling than proxy records”.\nThere is proxy support here through Fe speciation (Sperling), Ce anomalies (Wallace), C isotope\nfractionation (Edwards) and fossil charcoal (Glasspool). In my opinion there is better proxy support for a\nPaleozoic stepwise shift in O2 than there is for a Neoprotrozoic one.\nThanks for this point. We have now modified the text as follows:\n“ The rise of land plants has also been proposed to drive a step jump in atmospheric oxygen levels (e.g.,\n13), consistent with the results from both biogeochemical modeling and proxy records.”\nLine 55: “with the traditional or commonly held view of largely invariant atmospheric oxygen levels during\nthe Proterozoic”. I agree, but is this view actually commonly held anymore? The authors themselves give\nseveral high profile examples of more dynamic O2. Perhaps this can be rephrased more confidently?\nGood point. We have now modified the text as:\n“Although our picture of Earth’s redox evolution is ever-developing, it is becoming increasingly difficult to\nexplain some of the recent insights from the rock record that suggest variable atmospheric oxygen (e.g., 14)\n12\nand are inconsistent with the traditional view of largely invariant atmospheric oxygen levels during the\nProterozoic (e.g., 1).”\nThe section describing the model does not introduce what is actually in the model and on what models it\nbuilds. Having equations in the methods is fine but we need to know what you have done here in simple\nterms – much of this is taken from previous work and can be easily explained as such.\nThanks for this point. We have now expanded the model description in the text. The following texts are\nnow added:\n“The production of oxygen through organic burial is the sum of oxygen production in the ocean and on\nland. The oceanic oxygen production is calculated by considering the strength of net primary production\n(NPP (Tmol year-1)) and the burial efficiency of organic matter. Similar to the previous studies, the\nterrestrial oxygen production includes the fire feedback which is parameterized as a function of pO carbon\n2,\ndioxide level (pCO ), and temperature. The rate of hydrogen escape is a function of atmospheric methane\n2\nlevel (pCH ), following the previous studies.”\n“The rate of atmospheric methane oxidation (J ) is similar to previous studies and it is a function of\nCH4,O2\npO , pCH ,and an oxidation parameter (Ψ). The rate of oceanic sink of oxygen is parameterized to be a\n2 4\nfunction of NPP and the amount of anoxia in the ocean, which is a function of pO . Finally, consistent with\n2\nprevious modeling results, the oxidative weathering of pyrite and organic matter was assumed to change\nwith pO .”\nhere? And no discussion at all in the main text.\nThanks for this point. We meant a wide range of parameter values that can represent various environmental\nconditions. The range of parameter values is described in Table S1. We have now clarified this point in the\ntext. The text is now modified as:\n“ … that under a wide range of parameter values (Table S1) that can represent a wide range of environmental\nconditions and a wide range of initial conditions (pO values)…”\nCorrected. Thanks.\nLine 100: Earth system states instead of Earth’s system states\nCorrected. Thanks.\nLine 149: “the increased”\nCorrected. Thanks.\nLine 154: Earth system rather than Earth’s system\nCorrected. Thanks.\nCorrected. Thanks.\nLine 171: during the Proterozoic\nCorrected. Thanks.\nLine 173: what is a ‘fully stable state’ as opposed to a stable state?\nThanks for this point. We meant a stable state. It has now changed to “stable state”.\nLine 179: “A higher oxygen value”?\nCorrected. Thanks. It has now changed to “higher oxygen values”.\nIs figure 1 produced with plants or without? Surely there should be two lines here for the plants and\nno plants system? Ah ok it is in the SI. This should probably be brought into the main figure, it’s\nimportant!\nThanks for the suggestion. We have now added the following figure to the text.\nFigure 4. The impact of plant evolution on the steady-state atmospheric oxygen level. According to\nour results, the rise of plants on land would have helped with stabilizing the atmospheric oxygen level near\nmodern value.\nMore on figure 1. Many other things have changed over earth history that should change this curve.\nE.g. emergence of continents, reducing CO2 levels, increasing solar flux, changing mantle inputs.\nHow do you take account of all this?\nThanks for this point. Yes, we fully agree that other factors could have potentially influenced the geologic\nhistory of atmospheric oxygen levels. Nonetheless, we consider the change in other factors would have\nresulted in a change in the strength of the main feedbacks for the global oxygen mass balance. For instance,\n14\nthe change in the tectonic activity or continental configuration can result in a change in the rate of oxygen\nconsumption on land (e.g., by modulating the amount of weatherable material). And that’s why we took a\nstochastic approach that accounts for the variation in the strength of the global sources and sinks of oxygen\nin response to changes in environmental factors. For instance, the fire feedback parameterization used in\nour model (which is similar to the previous studies) accounts for the change in pCO level, and in our\n2\nstochastic simulation, we considered the full range of pCO suggested for the Phanerozoic time when the\n2\nrise of plants would have given a rise to the fire feedback.\nWe have now discussed this point in the text. The following text is now added:\n“Finally, while our model accounts for the major feedbacks on the global oxygen mass balance, there are\nother environmental factors (e.g., tectonic activity, continent configuration, etc.) that likely impact\natmospheric oxygen levels and have varied throughout Earth’s history. However, these factors would have\nlikely impacted some of the major feedbacks on the global oxygen mass balance that we so considered.\nFurther, our stochastic approach would likely capture, at least in part, the effect of these factors in altering\nthe feedback strengths.”\nLine 203: what is a “fully flat stable state”?\nThanks. We have revised the text for this section. The revised text now reads as:\n“The evolution of plants in the Phanerozoic likely helped stabilize atmospheric oxygen around modern\nvalue. To explore the impact of plant evolution on the atmospheric oxygen level, we conducted two\nstochastic simulations with and without plants. We explored the effect of the rise of plants by removing the\nflux of carbon burial from land plants from the total oxygen production flux (J equation 2 in the methods\nNPP;\nsection). According to our results, the median atmospheric oxygen level without plants would be around\n0.4 atmospheric oxygen level (PAL). This result shows that the rise of plants resulted in an increase in the\noxygen production flux and helped with stabilizing the atmospheric oxygen level near modern value (Figure\n4). Similar to the results from previous studies (e.g., 13), our results also suggest that, by exerting negative\nfeedback such as fire, plants may also help limit large swings in atmospheric oxygen levels.”\nFinal paragraph “According to our results, the atmospheric oxygen level could have fluctuated greatly for\nmost of the Earth’s history until the transition into the last stable state in the late Proterozoic and early\nPhanerozoic, where the oxygen level finally reached its modern value”. This is consistent with your model\nbut so are a great many other hypotheses. Indeed the abstract seems at odds here as it favors the low\nbackground with a few notable excursions.\nWe have modified the abstract to better account for the point discussed in the main text. Thanks for the\nsuggestion. The final sentence of the abstract is now modified as:\n“Our results, collectively, suggest a possible large dynamics in the atmospheric oxygen level throughout\nEarth’s history.”\n15\n-------------------------------------------------------------------------------------------------------------------------------\nReviewer #3\nReview of ‘Dynamics of Atmospheric Oxygenation’ by Fakhraee & Planavsky\nOverview/key results\nThis study utilizes a dynamical analysis approach to constrain O2 levels over Earth history. The results\nshow three stable steady state solutions, at pO2 values of ~10-6 PAL, ~10-3 PAL, and ~1 PAL, and two\nunstable steady state solutions between. The study applies these solutions to constrain oxygen levels over\nEarth history, comparing them to the values summarized in an existing review paper (see Figure 3 in\nparticular). Particularly striking implications of this part of the study are modern-like O2 values coincident\nwith the Lomagundi-Jatuli carbon isotope excursion (~2.3—2.0 Ga) and low pO2 during the mid-\nProterozoic (though the former is not highlighted in the study). The manuscript claims that the results\nindicate that pO2 could have fluctuated during the mid-Proterozoic between 10-4 and 10-2 PAL.\nThis is a really interesting study, which applies a powerful approach to an important research area.\nThanks for the time and feedback on our work.\nReviewer’s summary using Nature Communications review criteria:\nPlease see ‘Comments/suggestions’ section for further details.\nValidity\nFurther clarity and detail are required to confirm the robustness of the data interpretation. For\nexample, Figure 1b contains important information necessary to judge the validity of the work\nbut is not adequately described.\nThanks for this point. We have modified this figure and discussed the results of this figure in the text. The\nmodified figure and its modified caption now appear as:\nModified figure:\n“Figure 1. The results of stochastic dynamical analysis of global oxygen mass balance. Panel (a) shows\nthe median of the results of the stochastic analysis for a wide range of environmental and initial conditions.\nPanel (b) corresponds to the full range of results for different pO steady states (dpO /dt = 0). Our results\n2 2\nsuggest five steady states for the atmospheric oxygen level. The box plot in panel (b) is the result of the\nstochastic simulation where the middle line in the box corresponds to the median value, the edges of the\nbox correspond to 25th (lower edge) and 75th (upper edge) percentile, and the ends of the whiskers\ncorrespond to the minimum and maximum values for each state. The results of the stochastic simulation\nwere filtered to be consistent with the results from geochemical proxies. We only considered results that\ninclude at least three steady-state points: a modern-like steady state, a steady state under anoxic condition,\nconsistent with geochemical proxies that indicate an anoxic atmosphere (e.g., mass-independent\nfractionation of sulfur (MIF-S)) before the GOE, and one steady-state point under the low-oxygen condition\nthat is largely reflective of an environmental conditions that is neither anoxic (lack of MIF-S signal) and\nnor is modern-like (~ 1 present atmospheric oxygen level (PAL)).”\nUpdated text (1):\n“The results from the stochastic analysis of the global oxygen mass balance suggest that atmospheric\noxygen can be in balance under essentially anoxic conditions (state I; ~8x10-9 (25th percentile) – 2x10-7\n(75thpercentile) with median ~ 3x10-8 and minimum and maximum values of ~1x10-13 and 4x10-7 PAL),\nvery low oxygen conditions (state II; 1x10-5 – 1x10-4 with median ~ 3x10-5 and minimum and maximum values of\n3x10-7 and 2x10-4 PAL), low oxygen condition (state III: ~ 1x10-3 – 1x10-2 with median ~ 3x10-3 and\nminimum and maximum values of ~8x10-5 and 3x10-2 PAL), relatively high oxygen condition (state IV: ~\n1x10-1 – 4x10-1 with median ~ 2x10-1 and minimum and maximum values of ~10-2 and 7x10-1 PAL), and high,\nnear modern oxygen conditions (state V: 9x10-1 – 1.1x100 with median ~ 100 and minimum and maximum\nvalues of 7x10-1 and 1.3x100 PAL) (Figure 1b).”\nUpdated text (2):\n“Our results indicate five steady states for the atmospheric oxygen levels (dpO /dt = 0) (Figure 1).\n2\nwhere the sign of dpO /dt changes from positive to negative. In these states, the system tends to remain in\n2\nequilibrium because the positive-to-negative transition indicates a decrease in pO values when starting\n2\nfrom values larger than the steady state, eventually reaching stability. Conversely, when starting from\nvalues smaller than the steady state, the positive-to-negative transition signifies an increase in pO values,\n2\npulling the system back toward equilibrium (states I, III, V; Figure 1 & 2). Unstable states are also steady-\nstate points (dpO /dt = 0) but differ in that the sign of dpO /dt changes from negative to positive (states II,\n2 2\nIV; Figure 1 & 2). In these states, the system tends to move away from equilibrium, as the negative-to-\npositive transition implies an increase in pO values when starting from values larger than the steady state,\n2\nleading the system away from stability. Conversely, when starting from values smaller than the steady state,\nthe negative-to-positive transition indicates a decrease in pO values, pushing the system further from\n2\nequilibrium.\nOur results from the stochastic simulation were filtered to be consistent with the geochemical proxies.\nSpecifically, we only considered the results from the stochastic simulation that include a steady state point\nunder high oxygen modern-like conditions, a steady state point under anoxic conditions (pO < 10-5 present\n2\natmospheric level (PAL)) which is consistent with the presence of anoxic atmosphere during Earth’s early\nstage (e.g., mass-independent fractionation of sulfur), and a steady-state point between the two end\n17\nmembers of high and anoxic conditions that would generally correspond to the geologic times where there\nis no evidence for anoxic nor high oxygen conditions.”\nSignificance\nWith sufficient justification, the results could be of significance to the field. Constraining pO2 values,\nparticularly for the Mesoproterozoic, is an area of high interest and debate, in which various kinds of\nmodeling and geochemical studies can come to considerably different conclusions. However, I think that\nsuch implications could be justified further and discussed further in the text. For example, Figure 3 suggests\nthat oxygen levels were roughly modern during the Lomagundi-Jatuli event, but this is not really mentioned\nin the text, despite being a striking implication of the results, and it is insufficiently justified. I have noted\nsome other similar suggestions in the ‘Comments/suggestions’ section, including small corrections to better\ntake into account previous work.\nThanks for this point. We have now discussed the implication of our results for Lomagundi-Jatuli event.\nThe following text is now added to the main text:\n“According to our results, high-oxygen conditions during the Lomagundi-Jatuli event (LJE) can be\nattributed to two stable states of III and V. Specifically, considering that there is no evidence of a large\noscillation of atmospheric oxygen during the LJE that can be linked to an unstable state of atmospheric\noxygen (e.g., state IV), the suggested high oxygen level during LJE can be explained by either an\natmospheric oxygen level in state III or V. The maximum value of atmospheric oxygen in state III is around\n10-2 PAL which suggests that an increase in the rate of oxygen production during LJE can be explained by\na rise in atmospheric oxygen level to the maximum value of state III. Alternatively, a high atmospheric\noxygen level during this time can be explained by reaching a modern-like atmospheric oxygen level (state\nV). Regardless, our results suggest large shifts in atmospheric oxygen levels across the LJE transition.”\nIn addition, the differences between this study and previous studies should be explained more clearly,\nbecause the significance of this study as a scientific advance in the context of other work is not\ncommunicated clearly.\nThanks for the suggestion. Following your suggestion, we have also extended our discussion on the\ncomparison of our results with the previous studies.\nThe modified text for the comparison of our results with previous works now reads as:\n“Our results on the multiple stable states are comparable with widely held views of atmospheric oxygen\nevolution. For instance, some foundational modeling work has suggested the presence of a bi-stability in\nthe atmospheric oxygen level, with one stable state less than 10-5 PAL, another stable state above 5×10-\nthese results where there are two stable states (states I and III) before and after around 10-5 PAL with an\nunstable state around 10-5 PAL (state II; Figure 1). Our results for the range of atmospheric oxygen at\nunstable state II (1x10-5 – 1x10-4 with median ~ 3x10-5 and minimum and maximum values of 3x10-7 and\n2x10-4 PAL) are largely consistent with the previous modeling results. Notably, the results from the\nphotochemical modeling suggest instability of atmospheric oxygen between 5x10-8 and 5x10-4PAL that are\nlargely consistent with our results for the range of atmospheric oxygen level at unstable state II 23.\nOur results for the range of stable state under low-oxygen condition is also consistent with the previous\nmodeling results. For instance, the results from a one-dimensional time-dependent photochemical model\nsuggest a stable range of atmospheric oxygen above ~ 5x10−4 to 5x10−3 PAL 23. This range is consistent\n18\nwith the range of atmospheric oxygen resulting from our stochastic simulation for state III (state III: ~\n1x10-3 – 1x10-2 with median ~ 3x10-3). The results from another modeling study also indicated the presence\nof three stable states, one under anoxic conditions, one under low-oxygen conditions, and one under high-\noxygen conditions28, which is roughly consistent with our results here. However, these studies did not focus\non the possibility and extent of the fluctuation in the atmospheric oxygen level under low-oxygen conditions,\nwhich is increasingly supported by recent lines of evidence (e.g.,14,32). Our results suggest the Earth system\ncan reach a steady state at high oxygen conditions (state IV; ~10-1 PAL), as suggested by previous work28,\nbut that such a state is largely unstable (Figure 2). Our results are consistent with the range of steady-state\noxygen levels suggested by the previous studies. However, our analysis can provide further insight into the\nnature of such steady-state points, the range of variations in atmospheric conditions within these Earth\nsystem states, and the transitions between the states.”\nData, methodology, and analytical approach\nThe dynamical analysis appears to be a good approach to understand oxygen levels over Earth history.\nModeling of atmospheric steady states can be an important tool to understand this problem. The approach\ntakes into account relevant sources and sinks for atmospheric oxygen, including atmospheric oxidation of\nmethane. The ‘Methods’ section is detailed and contains most of the important information. There are a few\nsmall details that could be included, and I have made a note of these in the ‘Comments/suggestions’ section\nbelow. The figures are relevant and well-chosen, but some important information required to understand a\ncouple of the figures fully is required, and some figures contain small errors that could be misleading (see\nbelow). The main results/output are shown in Figure 1b, but (as mentioned above) there is an\ninsufficient description of this information, which makes it difficult to assess how well the results support\nthe conclusions. Figure 1b is particularly important, as I think it displays the spread in possible steady state\nsolutions and therefore gives an indication of the uncertainty, but it is insufficiently explained.\nThanks for this point. We have now, as explained above, updated the figure to better communicate the\nmodeling results with the audience as well as provided additional information in the text and the caption of\nthe figure.\nClarity, context, and references\nThe article is concisely written and very logically organized. The abstract is clear and accessible. The\ncontent is pitched at a level suitable for the broad readership of Nature Communications. The manuscript\ndiscusses some previous modeling work that also finds stable and unstable steady state solutions and places\nthis work in that context.\nThank you for the feedback and the time to review our work.\nPlease find more detailed comments and suggestions below.\nComments/suggestions\nMain text:\nThere are some grammatical errors throughout that occasionally obscure the clarity of the\narguments.\nL54-56: “it is becoming increasingly difficult to reconcile insights from the rock record with the\ntraditional or commonly held view of largely invariant atmospheric oxygen levels during the\nProterozoic.” I suggest that this sentence is reworded, as it is not clear whether references 17\nand 24 are given to support the “it is becoming increasingly difficult to reconcile insights” claim\nor given as references for the “traditional or commonly held view.” If the latter, the view of\nhigher oxygen levels during the Proterozoic also counts as “insights from the rock record” so I\nthink this sentence should be rephrased to acknowledge this. In other words, there is evidence\nfor higher oxygen levels as well as lower oxygen levels; it isn’t just a matter of evidence in the\nrock record vs. the “traditional” view.\nThanks for this point. We meant that the recent evidence is hard to explain with the traditionally accepted\nview of invariant atmospheric oxygen levels during the Proterozoic. We have now rephrased this sentence\nand modified the references.\n“ Although our picture of Earth’s redox evolution is ever-developing, it is becoming increasingly difficult\nto explain some of the recent insights from the rock record that suggest variable atmospheric oxygen (e.g.,\n14.) with the traditional view of largely invariant atmospheric oxygen levels during the Proterozoic (e.g.,\n1).”\nL96: “Our results were fairly stable… choosing each modeling input parameter.” What is meant\nby “fairly stable,” and why does this mean that “the simulation number is sufficient”? Also, the\nnumber of models run should be stated more clearly.\nThanks for this point. We meant that the number of runs for the stochastic simulation was 1000 and changing\nthe run numbers beyond 1000 did not influence the outcome of the stochastic simulation which suggests\nthat our number of runs is sufficient to capture the full range of values for each model input parameters.\nWe have now further clarified this point in the text. The modified text now reads as:\n“The number of runs for our stochastic simulation was 1000. Nonetheless, no noticeable difference in the\nmodel results was observed for a higher number of runs, which suggests that our selected number of runs\n(1000) is sufficiently large to capture the variability in the range of model input parameters.”\nL113: “According to our results, the only stable states occur under States I, II, and V.” There is a\nmistake in this sentence---it should be “under States I, III, and V.” Also, it would be helpful to\npoint out that stable steady states are defined as those when d 2 (pO2)/dt2 < 0 and that unstable\nsteady states are defined as those when d2 (pO2)/dt2 > 0 (perhaps using equivalent qualitative\nlanguage instead of these equations).\nCorrected. Thanks. The corrected sentence now reads as:\n“According to our results, the only stable states occur under states I, III, and V.”\nAnd thanks for the suggestion for making further clarification regarding the definition of a stable vs an\nunstable state in the text. We have now added the following text to the results section:\n“Generally, stable states (states I, III, V; Figure 1 & 2) are defined as the steady-state points (dpO /dt = 0)\n2\nwhere the sign of dpO /dt changes from positive to negative. In these states, the system tends to remain in\n2\nequilibrium because the positive-to-negative transition indicates a decrease in pO values when starting\n2\nfrom values larger than the steady state, eventually reaching stability. Conversely, when starting from\nvalues smaller than the steady state, the positive-to-negative transition signifies an increase in pO values,\n2\npulling the system back toward equilibrium (states I, III, V; Figure 1 & 2). Unstable states are also steady-\nstate points (dpO /dt = 0) but differ in that the sign of dpO /dt changes from negative to positive (states II,\n2 2\nIV; Figure 1 & 2). In these states, the system tends to move away from equilibrium, as the negative-to-\npositive transition implies an increase in pO values when starting from values larger than the steady state,\n2\nleading the system away from stability. Conversely, when starting from values smaller than the steady state,\n20\nthe negative-to-positive transition indicates a decrease in pO values, pushing the system further from\n2\nequilibrium.”\nL128: “While our results are in accordance with these modeling results, they further suggest a\npossibility of an unstable state around 10-4 PAL which allows for a possible fluctuation in the\natmospheric oxygen level around this level.” The manuscript includes comparison to Goldblatt\net al. (2006), which is very important. However, I think this sentence should be clarified,\nbecause it is not clear to me how the results here differ from Goldblatt’s finding of two stable\nsteady states with an unstable steady state solution between. Sorry if I’m just missing\nsomething.\nThanks for this point. Our results are similar to the results of Goldblatt et al. (2006). The text is now\nmodified as:\n“For instance, some foundational modeling work has suggested the presence of a bi-stability in the\natmospheric oxygen level, with one stable state less than 10-5 PAL, another stable state above 5 × 10-3 PAL,\nand an unstable state in between these two stable states 31. Our results demonstrate consistency with\nthese results where there are two stable states (states I and III) before and after around 10-5 PAL with an\nunstable state around 10-5 PAL (state II; Figure 1).”\nL185: “This conclusion is in accordance with geological evidence (e.g., sulfur isotope\nsystematics, Cr isotope, among other redox proxies.” I suggest this is reworded as “in\naccordance with some geological evidence” or equivalent, as some geological evidence\nsuggests pO2 above 1% PAL.\nThanks for this suggestion. We have followed the suggestion, and the text is now modified as:\n“This conclusion is, arguably, in accordance with some geological evidence (e.g., sulfur isotope\nsystematics14, Cr isotope 35, among other redox proxies 1).”\nL201: The difference between the model without plants and the model with plants is not clearly\nexplained. The parameter changed between these tests should be mentioned more explicitly\nhere and/or in the Methods.\nThanks for the suggestion. The discussion on the effect of the rise of plants on atmospheric oxygen levels\nis now modified as:\n“The evolution of plants in the Phanerozoic likely helped stabilize atmospheric oxygen around modern\nvalue. To explore the impact of plant evolution on the atmospheric oxygen level, we conducted two\nstochastic simulations with and without plants. We explored the effect of the rise of plants by removing the\nflux of carbon burial on land from the total oxygen production flux (J equation 2 in the methods section).\nNPP;\nAccording to our results, the median atmospheric oxygen level without plants would be around 0.4\natmospheric oxygen level (PAL). This result shows that the rise of plants resulted in an increase in the\noxygen production flux and helped with stabilizing the atmospheric oxygen level near modern value (Figure\n4). Similar to the results from previous studies (e.g., 13), our results also suggest that, by exerting negative\nfeedback such as fire, plants may also help limit large swings in atmospheric oxygen levels.”\nWe have also added a new figure (below) to the main text that shows the effect of rise of plants on the\natmospheric oxygen level:\nour results, the rise of plants on land would have helped with stabilizing the atmospheric oxygen level near\nmodern value.\nL203: “possible fully flat stable state for oxygen levels higher than modern value.” The meaning\nof “fully flat stable state” should be clarified. Also, looking at Figure S2 (left hand), it looks like the d(pO2)/dt is negative at pO2 > 1 PAL, rather than zero, so I don’t understand why this is\nconsidered a stable state.\nThanks for this point. We meant a stable state by a fully stable state. We have now removed the term “fully”\nto avoid confusion. Regarding your point about Figure S2, we have now replaced this figure with a new\nfigure in the main text (shown above).\nL210: “… but these results suggest plants may also help limit large swings in atmospheric\noxygen levels.” This is not a novel result, so this clause should be slightly reworded, as it\ncurrently sounds like this study was the first to discover this.\nThanks for this point. Yes, this is not a new result and is a confirmatory result that has been shown by\nprevious studies. We have now modified the text as follows:\n“Similar to the results from previous studies (e.g., 13), our results also suggest that, by exerting negative\nfeedback such as fire, plants may also help limit large swings in atmospheric oxygen levels.”\nMethods:\nL255: Jocean,sink is repeated. One should be JCH4,O2.\nCorrected. Thanks.\nL297: “To account for uncertainty… for the range of each parameter.” This is a really helpful\nsentence. If the uniform distribution applies to all parameters in Table S1, rather than just the\nparameters related to JNPP,ocean, then this sentence could be put further towards the top of the\nMethods section, rather than here where it is more hidden.\nwith regard to the stochastic simulation:\n“It should also be noted that to account for the uncertainty associated with the choice of values of the model input\nparameters, we took a stochastic approach wherein the values for each modeling input parameter were selected\nrandomly assuming a uniform distribution and the most probable range of results were obtained (Figure 1b).”\nL322: A sentence could be added to this section to explain how values for T are obtained. Does\nT change, or is a fixed value chosen?\nThanks for this point. The value of “T” was changed within the suggested range of temperature during\nPhanerozoic.\nL322: pCO2 is defined twice, and pCO2,1 is not defined.\nCorrected. Thanks.\nThe modified text now reads as:\n“Where T is the temperature in the unit of degrees Celsius, pCO is CO atmospheric level in ppm, pCO denotes\n2 2 2,min\nthe minimum atmospheric CO concentration in which plants can grow, and pCO is a constant13,30,46,47.”\n2 2,1\nFigures:\nFigure 1: The caption does not sufficiently explain Figure 1b. For example, the meaning of the\nhorizontal width of the shapes is not explained. The vertical black lines and the thin black lines\nat the edges of the shapes are not defined. I assume that the white circles show the median\nvalues, but this should be stated explicitly.\nWe have now updated Figure 1 so that the results of the stochastic analysis would be better communicated\nto the reader. We also provided additional information regarding the results in Figure 1b (box plots) as well\nas details regarding the stochastic analysis.\nFigure 2: The caption states that the vertical axis label “Energy” is the potential energy, but as a\nreader I am unclear about what that means. I do not understand why the different stable and\nunstable states have different potential energies nor how that is calculated. I suggest that more\ndetail is given here, both in the caption and in the text.\nThanks for this point. We provided additional details regarding the calculation of energy. The following\ntexts are now added to the caption of Figure 2 as well as to the results section which includes information\nabout the calculation of energy.\nUpdate caption:\n“Figure 2. The schematic presentation of the different oxygen steady states. The result from the dynamical\nanalysis of a simple oxygen mass balance suggests four steady states (dpO /dt = 0) for the atmospheric\n2\noxygen level. The first one is a stable state (State I) under anoxic conditions, the second one (State II) is an\nunstable point, and the third one is the first stable state for oxygen (State III) with an oxygen range of ~10-\n3 present atmospheric level (PAL). The fourth steady state point is unstable is ~ 10-1 PAL (State IV). The last\nstable state is State V with a modern-like atmospheric level. The “Energy” corresponds to the potential\n23\nenergy. The energy was calculated as the surface area under the curve (Figure 1a) which is the integral of\ndpO /dt over the range of pO (-∫\n𝑝𝑂2,𝑓 𝑑𝑝𝑂2\n𝑑𝑝 ). The amount of energy depends on the extent of the basin\n2 2 𝑝𝑂2,𝑖 𝑑𝑡 𝑂2\nof attraction (a region in which the initial conditions (pO ) lead to the eventual convergence of the system\n2\nto a specific stable steady state or attractor). State I with the largest basin of attraction has the minimum\nenergy implying a large stability around this state.”\nAdded text to the results section:\n“To explore the extent of stability in each stable state, we calculated the energy corresponding to each\nsteady state point. The energy was calculated as the surface area under the curve (Figure 1a) which is the\nintegral of dpO /dt over the range of pO (-∫\n𝑝𝑂2,𝑓 𝑑𝑝𝑂2\n𝑑𝑝 ) (Figure 2). The amount of energy depends on\n2 2 𝑝𝑂2,𝑖 𝑑𝑡 𝑂2\nthe extent of the basin of attraction (a region in which the initial conditions (pO ) lead to the eventual\n2\nconvergence of the system to a specific stable steady state or attractor). For instance, state I with the largest\nbasin of attraction has the minimum energy implying a large stability around this state (Figure 2).”\nFigure 3b: This is a summary/implications figure, so it is important that it accurately reflects the results.\n• The caption does not sufficiently explain the meaning of the blue band and the blackoutlined box – e.g.,\ndo they encapsulate all the models run, or 95% of them, etc.? A really helpful thing to note would be the\nrelationship of the width of these bands to the shapes and/or black vertical lines on Figure 1b, as it is difficult\nto work out where the width of the bands come from.\nThanks for this point. We have now updated this figure and the caption of the figure. The range of the blue\nband is the same as the range of values presented in Figure 1b in the box plot. The lower and upper edges\nof the blue band correspond to the 25th and 75th percentile of the results from the stochastic simulation.\nWe have now included this information in the figure.\nThe updated figure now appears as:\nThe updated caption now reads as:\n“Figure 3. The geological history of atmospheric oxygen level. Panel (a) corresponds to the records of\nmass-independent fractionation of sulfur (MIF-S)5,34, as a proxy for an anoxic atmosphere. And (b) shows\n24\nthe change in atmospheric oxygen through time. The blue boxes with a line denote the trajectory for\natmospheric oxygen level based on the results in this study and the grey boxes correspond to the traditional\ntrajectory of atmospheric oxygen level 1,50. The range of values in blue boxes is the same as the range\npresented in Figure 1b where the lower and upper edge of each blue box correspond to the 25th and 75th\npercentile of the data that resulted from the stochastic simulation. Our results suggest a possibility of a\nlarge fluctuation in the atmospheric oxygen level for most of the Earth’s history. The rise of atmospheric\noxygen in the late Proterozoic and early Phanerozoic, coincident with the rise and radiation of the complex\nlife marks the end of the fluctuations in the atmospheric oxygen level and transition into stable State V with\nmodern-like oxygen level. The high atmospheric oxygen during Lomagundi, after GOE, can be explained\nby either atmospheric oxygen in state V (~100 PAL) or reaching the maximum value for the atmospheric\noxygen level in stable state III (~10-2 PAL).”\n• Relatedly, it is important that the components of Figure 3b are drawn accurately. For example, the lower\nbound of the mid-Proterozoic band (this study) is ~8e-5 PAL, but the lower band of the State III shape is\n2e-4, and the lower limit of the black vertical line is even higher. Though the difference between 8e-5 and\n2e-4 PAL is not huge, it could be a little misleading as drawn, as the mid-Proterozoic box appears to suggest\nslightly lower pO2 values than the results support.\nThanks for this suggestion. We have now modified the figure to better represent the range of results\npresented in Figure 1b.\n• Very minor: Similarly, the medians and black-outlined boxes of the 2nd (Lomagundi) and 4th\n(Phanerozoic) boxes should be at exactly the same pO2, but a careful examination shows that they are a\nlittle offset. This is extremely minor, but should be tidied up before publication. The lower limit of the mid-\nProterozoic curve from Lyons et al. (2014) is 1e-4 PAL, while it is drawn slightly higher here.\nThanks for this point. We have now modified the figure to address this point.\n• The figure as drawn suggests that oxygen levels were ~modern during the time of the\nLomagundi excursion, shown in the figure between 2.3 and 2 Ga. However, the State III constraints\n(pO2~1e-3 PAL) could have alternatively been used for this period of Earth history. I assume that they\nwere not used in the figure because of evidence mentioned in the text of pO2 > 1% PAL (lines 33-34),\nbut since the blue band lies outside the Lyons et al. (2014) constraints box anyway, it is not clear why the\npossibility of the State III levels during this time is not included. There should be some discussion of this\nparticular period of Earth history in the text. If leaving Figure 3 as it is to show only the State V levels,\nthe discussion should comprise either: (1) a justification of ~modern pO2 levels for this period, explicitly\nstating why pO2 < 1% can definitely be dismissed, as it would be a really important\nimplication of the current work if this can really be robustly defended, or (2) the limitations\nand caveats of the results’ support of State V oxygen levels only for this time.\nThanks for this point. While the community consensus around the atmospheric oxygen level during\nLomagundi is yet to be achieved, we suggest two possible explanations for the atmospheric oxygen during\nthis time. One possible explanation could be the atmospheric oxygen level during time has reached the\nmaximum value in stable state III (~10-2 PAL) and the alternative explanation would be to near modern\natmospheric oxygen level (state V). We have now added a discussion about these two possibilities for\natmospheric oxygen level during Lomagundi. We have also added a question mark for the atmospheric\noxygen during Lomagundi to reflect the uncertainty surrounding the atmospheric oxygen level during time.\nWe have also included a mention of the uncertainty for the atmospheric oxygen level in the caption of the\nfigure.\nUpdated caption:\n“Figure 3. The geological history of atmospheric oxygen level. Panel (a) corresponds to the records of\nmass-independent fractionation of sulfur (MIF-S)5,34, as a proxy for an anoxic atmosphere. And (b) shows\nthe change in atmospheric oxygen through time. The blue boxes with a line denote the trajectory for\natmospheric oxygen level based on the results in this study and the grey boxes correspond to the traditional\ntrajectory of atmospheric oxygen level 1,50. The range of values in blue boxes is the same as the range\npresented in Figure 1b where the lower and upper edge of each blue box correspond to the 25th and 75th\npercentile of the data that resulted from the stochastic simulation. Our results suggest a possibility of a\nlarge fluctuation in the atmospheric oxygen level for most of the Earth’s history. The rise of atmospheric\noxygen in the late Proterozoic and early Phanerozoic, coincident with the rise and radiation of the complex\nlife marks the end of the fluctuations in the atmospheric oxygen level and transition into stable State V with\nmodern-like oxygen level. The high atmospheric oxygen during Lomagundi, after GOE, can be explained\nby either atmospheric oxygen in state V (~100 PAL) or reaching the maximum value for the atmospheric\noxygen level in stable state III (~10-2 PAL).”\nAdded text on Lomagundi to the main text:\n“According to our results, high-oxygen conditions during the Lomagundi-Jatuli event (LJE) can be\nattributed to two stable states of III and V. Specifically, considering that there is no evidence of a large\noscillation of atmospheric oxygen during the LJE that can be linked to an unstable state of atmospheric\noxygen (e.g., state IV), the suggested high oxygen level during LJE can be explained by either an\natmospheric oxygen level in state III or V. The maximum value of atmospheric oxygen in state III is around\n10-2 PAL which suggests that an increase in the rate of oxygen production during LJE can be explained by\na rise in atmospheric oxygen level to the maximum value of state III. The end of the LJE can, then, be\n26\nmarked as a transition into a lower atmospheric oxygen level in state III. Alternatively, a high atmospheric\noxygen level during this time can be explained by reaching a modern-like atmospheric oxygen level (state\nV). Regardless, our results suggest a possibly large dynamic in atmospheric oxygen levels across the LJE\ntransition.”\nSupplementary Information:\nTable S1: “Coefficient for ocean anoxic sink.” If there is no reference for the range of numbers\nchosen for k1 (0.01—0.1), where did this range come from? I see that in L373 in the Methods, it\nis stated that it is varied “within a wide range” to account for uncertainty, but I think some\njustification for this particular range should be included.\nThanks for this point. We have now included additional discussion for this coefficient in the methods.\nThe modified text now reads as:\n“The constant “k1” is a factor that accounts for variability in the pool of reducing compounds (e.g., sulfide, iron,\nmethane) in the ocean. The constant “k1” reflects the relative comparison between the rate of oxygen consumption by\nreducing species. For instance, the value of 0.1 means that the rate of oxygen consumption by reducing species is\nabout an order of magnitude lower than the maximum rate of anaerobic respiration. Under anoxic and low-oxygen\nconditions of Archean and Proterozoic, this value would have likely been low. For instance, most of the sulfide\nresulting from sulfate reduction would have precipitated as pyrite, consistent with the results from sulfur isotope\nrecords that suggest a near-one value for the fraction of pyrite (f ). Our results for the J for the modern-like\npyr ocean,sink\ncondition are consistent with the rate of oxygen uptake in the modern oxygen minimum zone. Finally, to account for\nthe uncertainty, this value along with other values used in this equation and elsewhere was varied within a wide range\nin our stochastic simulation (Table S1).”\nFigure S3: The plots presented here contain a lot of useful information, but I can’t see any\ndiscussion of what they show, either in the main text or the SI, and the figure is not referred to in\nthe text. It would be beneficial to add some sentences describing the main takeaways from\nthese plots, because there is some interesting information there. In addition, a longer caption\nwould be helpful.\nThanks for this point. We have now cited this figure in the text and provided additional details regarding\nthis figure in the caption.\nThe caption is now modified as:\n“ Figure S3. Results of stochastic analysis for oxygen sources and sinks. The blue band corresponds to\n95% of the results, and each corresponds to roughly 10% of the results. The results show a change in the\nstrength of various feedbacks on the global oxygen mass balance as a function of different pO levels. The\n2\nflux of oxygen production through net primary production (J ) would increase with pO and reach a\nNNP 2\nnear-modern value around 100 present atmospheric oxygen level (PAL). The rate of oxidative weathering\nwould only reach its maximum at pO above 10-4 PAL, consistent with the previous modeling results. High\n2\nmethane concentration under a widely anoxic condition (<10-6 PAL) results in a high hydrogen escape and\nwith an increase in pO level, its magnitude becomes smaller and its role in the global oxygen mass balance\n2\nbecomes less important.”\nThank you very much for an interesting paper to review. I hope to see it published in due course.\nThank you very much for the helpful suggestions and comments.\nReviewer #1 (Remarks to the Author):\nThe updated draft addresses most of my concerns with the submitted manuscript. I have no more\nmajor comments, except the Zenodo archive of the source code is not very well documented. I\nrecommend the archive improve the \"README\" file with a full description for how to reproduce the\nresults. Besides this concern, I recommend the article for publication.\nReviewer #1 (Remarks on code availability):\nI reviewed the code repository and tried to run the uploaded script. Unfortunately, the code is\nwritten in Matlab which is a proprietary language, and not easily accessible. A free partial Matlab\ninterpreter exists, called Octave. I installed Octave, then tried to run main script, and the code\nworked up until the plotting routing, then errored out because some of the plotting required a\nlegitimate proprietary Matlab installation.\nAlso, the \"README\" in the code archive gives almost no information about how to run the script. A\nbetter \"README\" would describe exactly the Matlab version to install in order to reproduce the\nwork.\nOverall, the repository is not very well documented/reproducible, and it would be ideal if it could be\nimproved.\nReviewer #2 (Remarks to the Author):\nMy concern with this work was that the processes that drove the newly-identified O2 instability\naround 0.1 PAL in their model were not investigated. The authors have now added some extra\nanalysis in response to my questions and I also appreciate them attaching the model code. The\nnew text was not marked in the revised ms which made it difficult to read in context, but it was\nsufficient – alongside the code - for me to understand what is going on in the model. The authors\nhave a term in their model for ‘ocean oxygen consumption’ which is a linear function of the degree\nof marine anoxia. The assumption here is that more marine anoxia leads to a greater concentration\nof reduced compounds to react with O2. The overall result in the model, as far as I understand it, is\nthat when atmospheric oxygen concentration starts to (for example) decrease, marine anoxia\nexpands, and this results in further consumption of O2 in a positive feedback loop - which creates\nthe instability around 0.1 PAL and allows stable states aroubnd 1 PAL and 0.01 PAL. This is\nsomewhat contrary to established ideas. Although marine anoxia may promote a greater\nconcentration of reduced compounds that consume O2, it also promotes a greatly increased burial\nrate of organic carbon – which should increase O2 production. This relationship is quite clear from\nthe geological record, most obviously in highly organic rich sediments deposited during OAEs. To be\nfair to the authors, they have included an increased burial efficiency of organic carbon under\nmarine anoxia, but the simplified equations used mean that the scaling of the O2 consumption\neffect is much stronger in the model than the scaling for burial efficiency. This is because today’s\nanoxic area is so small and the authors assume that this area already conmsumes a significant\namount of O2, so even small anoxic areas rapidly increase the O2 consumption flux. Moreover, they\nauthors have not included a phosphorus cycle – and it is well understood that anoxia promotes P\nrecycling from sediments which drives greater productivity and atmospheric O2 production. \nVanCappellen and Ingall (1996, Science) proposed that this P recycling feedback stabilized\natmospheric and marine oxygen and prevented runaway anoxia, which is the opposite of the\nprocess proposed here by the authors.\nPersonally, I would tend to believe that the overall feedback aorund these O2 levels is stabilizing (as\nper Vancappellen and Ingall). There is tentative direct evidence for an increase in wildfire activity –\nhinting at rising O2 - during the Toarcian oceanic anoxic event (Baker et al., 2017 – Nat Comm),\nwhich supports their hypothesis (marine anoxia -> rising atmos O2) rather than the model proposed\nhere. That said, we are comparing between two highly simplified box models and very limited direct\nevidence, so it is not easy to decide which is better. Overall, I think this paper can be further\nmodified to give a clearer explanation of what is going on in the model around this transition. I think\nit is important to justify the chosen rate of consumption of O2 by reduced marine compounds at the\npressent day, given the small extent of anoxia area. It can also be clearly caveated that the key\nresults here depend on this process dominating O2 production through increased burial efficiency\nwhen anoxia rises, and also there being no significant phosphorus recycling in anoxic oceans. Of\ncourse, in a low-sulfate conditions there may well be minimal phosphorus recycling anyway, so\nthere is some leeway there. Overall I think the idea that positive feedbacks could outweigh negative\nones under moderately low O2 levels is something worth thinking about, so I would support\npublication of a paper revised as above, unless I have some of those points wrong of course in\nwhich case I am happy to be corrected.\nMinor points:\n- I would suggest modifying the title a little to “large-scale dynamics…”\n- Figure 1: states V and IV are labelled the wrong way around it seems?\nReviewer #3 (Remarks to the Author):\nReview of revision of “Large dynamic in the geological atmospheric oxygen level”\nThe revisions made have greatly improved the manuscript. I have a few further comments and\nsuggestions, described below. Thank you very much again for an interesting manuscript to review.\nComments and suggestions\nL79: “following the previous studies” – I suggest removing “the.”\nL88: Repeated sentence here.\nL107—108: “will be sampled randomly within the range suggested in the literature (Table S1), assuming a\nuniform distribution.” It is great that the code is now accessible. However, now that I have been able to\nlook at the file ‘OxygenDynamic.m,’ I am a little concerned about the way the parameters are chosen. I’m\nnot particularly fluent in Matlab so may have made a mistake, but it looks like some parameters are not\nsampled randomly with a uniform distribution, as claimed in the text. For example, the code samples pyrite\nweathering factors from the array ‘Weathering_pyrite_MC,’ which has values of [1500 1700 2000 2200\n2500], which are not evenly distributed. Similarly, ‘pCO2_1_MC’ is [100 200 280 300 500 600 700 1000].\nSee also ‘Temp_MC.’ Apologies if I have misunderstood the code or am looking at the wrong file, but I think\neither the models should be run with uniform distributions of these variables, or the non-uniform\ndistribution should be described and justified in the text. This method of sampling could affect where the\nmedian line lies and therefore where the steady states lie. See also L352.\nAnother apparent inconsistency between the code and the model description in the text is that the variable\nFw_oc comes from randomly sampling values from the list [1 1.5 2 3 3.5 4 4.5 5 7 10], but I think this is\nthe same as F in Table S1, which is said to be randomly sampled from a uniform distribution between the\norg\nvalues of 3.75 and 7.75. Again, apologies if I’ve misunderstood.\nL177—207: Two paragraphs from “Our results from the sensitivity analysis suggest…” to “…depending\non the directionality of the change.” The addition of discussion of these results in the main text is\nimportant and very interesting. However, I have a few suggestions about this section:\n• L177: A sentence or two introducing the sensitivity analysis and describing what was done should\nbe given here, as it hasn’t been mentioned earlier in the paper. This could begin with the first\ncouple of sentences from the section in the Methods (see comment for L396-506 below).\n• L177-182: “the rate of hydrogen escape and the atmospheric oxidation of methane have the largest\ninfluence on the emergence of the first and second steady-state for pO2 under anoxic conditions” and\n“our results suggest that the first steady-state point is controlled by the balance between the sink of\noxygen from reaction with reduced volcanic gases and the source of oxygen through the escape of the\nhydrogen to the atmosphere.” Also L197-199: “Of the two steady-state points under low-oxygen\nconditions (states III & IV), the first one is the result of an increase in oxygen production in the wake of\nthe evolution of oxygenic photosynthesis while the rate of oxidative weathering is not high.” It’s not\nclear to me how these statements are obtained from the results presented in Figure S2. More\ndetails about what Figure S2 shows and how the results support these three statements would\nstrengthen the text here and help the reader follow the argument from the figures and results to\nthese implications.\n• “An increase in the rate of hydrogen escape would lead to the disappearance of the first steady-state\npoint (state I).” I don’t follow how this is seen from the results in Fig. S2, so suggest more explanation is\ngiven here to help the reader.\nL179: “…consistent with the previous study.” The previous study in mind should be referenced here.\nMethods:\nL347: J is defined as volcanic outgassing in L81 of the main text, but it would be useful if it was also\nvolc\ndefined here with the other O sinks. \nH2 hydrogen\nmethane oxidation.\nL354—494: There is no section detailing how J (in Equation 1) was parameterized.\nvolc\nL396—506:\n• This paragraph is a good addition to the Methods. I suggest that the content of this paragraph\n(particularly the first two sentences) is put in the main text instead. The first two sentences could\ngo at the beginning of the paragraph at L177 before “Our sensitivity analysis…”\n• L497-499: This is a helpful sentence, which should be in the main text as mentioned above. It would\nbe helpful to add whether 1000 runs were done for each case, and if not, how many were run.\nFigures:\nFigure 1b: The caption for this figure has been greatly improved. However, a description of what the\ncircular points represent is missing.\nSI:\nFigure S1 caption: It would be helpful to reference the “previous photochemical modeling” explicitly\nhere, either by using the numbered reference(s) or an in-text citation (e.g., “The circles represent data\nfrom Name et al. (20XX)”).\nFigure S2: A sentence describing the blue line could be added to the caption for clarity, e.g., “The blue\nline shows the…”\nFigure S3: Is J the same as J ? The title on the y axis of the upper right panel should be either\nanoxia ocean,sink\ndefined or changed.\nGeneral comments:\nSome grammatical errors remain, which can obscure the clarity of the argument.\nReviewer #3 (Remarks on code availability):\nI did not run the code myself, but I read the scripts.\nAs noted in my main review file, my main concern with this revision of the manuscript is that there\nappear to be inconsistencies between the code file and the description of the model in the text. In\nparticular, it doesn't look like the parameter values are sampled randomly with a uniform\ndistribution. I may have made a mistake in reading the code, but recommend that the authors\ncheck this and make any necessary changes before publication.\nResponse to Reviewers\nReviewers’ comments are in black.\nOur responses are in blue.\n-------------------------------------------------------------------------------------------------------------------------------\nReviewer #1 (Remarks on code availability):\nI reviewed the code repository and tried to run the uploaded script. Unfortunately, the code is written in\nMatlab which is a proprietary language, and not easily accessible. A free partial Matlab interpreter exists,\ncalled Octave. I installed Octave, then tried to run main script, and the code worked up until the plotting\nrouting, then errored out because some of the plotting required a legitimate proprietary Matlab installation.\nAlso, the \"README\" in the code archive gives almost no information about how to run the script. A better\n\"README\" would describe exactly the Matlab version to install in order to reproduce the work.\nOverall, the repository is not very well documented/reproducible, and it would be ideal if it could be\nimproved.\nThank you for this point. We have now updated the README file to include further details regarding the\nversion of the MATLAB needed as well as an additional function to plot the results.\nThe updated “README” file now reads as:\n“# OxygenDynamics\nThe code calculates the range of steady-state solutions for the global oxygen mass balance taking a\nstochastic approach.\nThe code contains a main loop to simulate the range of steady-state in a stochastic manner. The loop\ncontains the range of values used to calculate the different oxygen sinks and sources.\nThe code can be run without requiring high-performance computing. The results from the stochastic\nsimulation would be plotted using a “fanchart” function. The function for “fanchart” is also included in\nthe file named “fanChart.m”. The “fanchart” function needs to be saved in the same folder as the main\ncode in order to plot the results.\nThe MATLAB code is written in MATLAB R2022b version.”\nMy concern with this work was that the processes that drove the newly-identified O2 instability\naround 0.1 PAL in their model were not investigated. The authors have now added some extra\nanalysis in response to my questions and I also appreciate them attaching the model code. The\nnew text was not marked in the revised ms which made it difficult to read in context, but it was\nsufficient – alongside the code - for me to understand what is going on in the model. The authors\nhave a term in their model for ‘ocean oxygen consumption’ which is a linear function of the degree\nof marine anoxia. The assumption here is that more marine anoxia leads to a greater concentration\nof reduced compounds to react with O2. The overall result in the model, as far as I understand it, is\nthat when atmospheric oxygen concentration starts to (for example) decrease, marine anoxia\nexpands, and this results in further consumption of O2 in a positive feedback loop - which creates\nthe instability around 0.1 PAL and allows stable states aroubnd 1 PAL and 0.01 PAL. This is\nsomewhat contrary to established ideas. Although marine anoxia may promote a greater\nconcentration of reduced compounds that consume O2, it also promotes a greatly increased burial\nrate of organic carbon – which should increase O2 production. This relationship is quite clear from\nthe geological record, most obviously in highly organic rich sediments deposited during OAEs. To be\nfair to the authors, they have included an increased burial efficiency of organic carbon under\nmarine anoxia, but the simplified equations used mean that the scaling of the O2 consumption\neffect is much stronger in the model than the scaling for burial efficiency. This is because today’s\nanoxic area is so small and the authors assume that this area already conmsumes a significant\namount of O2, so even small anoxic areas rapidly increase the O2 consumption flux. Moreover, they\nauthors have not included a phosphorus cycle – and it is well understood that anoxia promotes P\nrecycling from sediments which drives greater productivity and atmospheric O2 production. \nVanCappellen and Ingall (1996, Science) proposed that this P recycling feedback stabilized\natmospheric and marine oxygen and prevented runaway anoxia, which is the opposite of the\nprocess proposed here by the authors.\nPersonally, I would tend to believe that the overall feedback aorund these O2 levels is stabilizing (as\nper Vancappellen and Ingall). There is tentative direct evidence for an increase in wildfire activity –\nhinting at rising O2 - during the Toarcian oceanic anoxic event (Baker et al., 2017 – Nat Comm),\nwhich supports their hypothesis (marine anoxia -> rising atmos O2) rather than the model proposed\nhere. That said, we are comparing between two highly simplified box models and very limited direct\nevidence, so it is not easy to decide which is better. Overall, I think this paper can be further\nmodified to give a clearer explanation of what is going on in the model around this transition. I think\nit is important to justify the chosen rate of consumption of O2 by reduced marine compounds at the\npressent day, given the small extent of anoxia area. It can also be clearly caveated that the key\nresults here depend on this process dominating O2 production through increased burial efficiency\nwhen anoxia rises, and also there being no significant phosphorus recycling in anoxic oceans. Of\ncourse, in a low-sulfate conditions there may well be minimal phosphorus recycling anyway, so\nthere is some leeway there. Overall I think the idea that positive feedbacks could outweigh negative\nones under moderately low O2 levels is something worth thinking about, so I would support\npublication of a paper revised as above, unless I have some of those points wrong of course in\nwhich case I am happy to be corrected.\nThanks for this feedback. We have now revised the parameterization scheme used for the sink of oxygen in\nthe ocean and followed the previous parameterization scheme (building foremost from Slomp and Van\nCappellen, 2007). In the revised version, we have explicitly considered the effect of seawater oxygen\nconcentration on the rate of oxygen consumption.\nparameterization scheme for organic matter burial to account for the enhanced phosphorous release under\nanoxic conditions. We followed the previous studies to include the impact of anoxia on phosphorous\nrecycling (e.g., Slomp and Van Cappellen, 2007; Alcott et al., 2019). We also included a factor that reflects\nthe impact of sulfate availability on the phosphorous release under anoxic conditions. This modification,\nhowever, did not result in a major change in the results presented in the text. The minor impact of anoxia\non phosphorous availability, as you also mentioned, is due to the low-sulfate conditions prevalent under an\nanoxic and weakly oxygenated atmosphere, which is also consistent with more limited phosphorous\nrecycling in a modern low-sulfate system (e.g., Xiong et al., 2019).\nAccording to our results, the emergence of an unstable state around 10-1 PAL is controlled by both limited\noxygen production and organic matter burial and a high rate of organic matter and pyrite weathering on\nland. Our results, as shown in Figure S3, suggest only a minor role in oxygen consumption compared to the\nsink of oxygen through pyrite and organic matter weathering as well as the source of oxygen from organic\nmatter burial. The emergence of an unstable state is the direct result of limited oxygen production in the\nabsence of land plants as well as an enhanced phosphorous recycling, caused by a low-sulfate condition,\nand increased rate of oxidative weathering on land (Figure S2).\nIn the context of the dynamical system approach, the presence of two consecutive stable states implies the\nemergence of an unstable state between these two stable states that acts as a transition between these two\nstable states. In the current case, the two stable states are modern-like atmospheric oxygen (stable V) and a\nstable state under low-oxygen conditions (stable III).\nWe have now modified the text both in the main text and the method section that addresses the points raised\nherein.\nModified text 1 (Main text):\n“The oceanic oxygen production is calculated by considering the strength of net primary production (NPP\n(Tmol year-1)), the burial efficiency of organic matter, and a factor that accounts for the effect of anoxia on\nphosphorous cycling and organic matter burial (see Methods).”\nModified text 2 (Main text):\n“Based on our results from the sensitivity test, the emergence of the third and fourth steady-state points\nunder low-oxygen conditions (states III & IV) is controlled by the oceanic oxygen production (J ) and the\nNPP\nsink of oxygen in the ocean and on land. An increased rate of oceanic oxygen production (J ) would shift\nNPP\nthe imbalance value of oxygen (dpO /dt) to positive values which would influence the appearance of steady-\n2\nstate points under low-oxygen conditions. At the same time, according to our results, while under anoxic\nconditions (<10-5 PAL) the oxidative weathering of pyrite and organic matter is essentially ineffective in\ncontrolling the global oxygen mass balance (Figure S2 and S3), an increased rate of oxidative weathering\nunder low-oxygen condition would dramatically change the oxygen imbalance value which would impact\nthe emergence of steady-state points under low-oxygen conditions. Of the two steady-state points under\nlow-oxygen conditions (states III & IV), the first one is the result of an increase in oxygen production in the\nwake of the evolution of oxygenic photosynthesis while the rate of oxidative weathering is not high. An\nincrease in the atmospheric oxygen level would result in an increase in the rate of oxidative weathering\nwhich would result in a negative imbalance for the oxygen. A continued increase in the production of oxygen\nin the ocean would have resulted in a decrease in anoxia in the ocean which would have acted as positive\nfeedback in the global oxygen mass balance and have resulted in the emergence of the state IV. Such a state\nis, however, unstable which means that a change in the oceanic oxygen mass balance (e.g., a change in the\n3\noxygen production through net primary production or marine oxygen consumption) would shift the system\nto the adjacent stable state (state III or state V), depending on the directionality of the change. Our results\nalso suggest a minor role of enhanced phosphorous recycling under anoxic conditions in driving the\nemergence of steady-state points under low-oxygen conditions, caused by low-sulfate conditions under low-\noxygen conditions.”\nModified text 3 (Methods):\n“The flux of carbon burial in the ocean, J , can be parameterized as:\nNPP,ocean\n𝐽 = 𝜀∗𝑁𝑃𝑃∗𝑎𝑟𝑒𝑎∗ 𝑅 (3)\n𝑁𝑃𝑃,𝑜𝑐𝑒𝑎𝑛 𝑐/𝑝\nWhere ε is the efficiency of the carbon pump, NPP is the net primary production, and the last term is the\narea of the ocean. The value of NPP is multiplied by a productivity factor to account for a lower productivity\nat lower oxygen levels40", "reviewer_41": "n𝑁𝑃𝑃 = 𝑁𝑃𝑃 ∗ 𝑓 (4)\n𝑚𝑜𝑑𝑒𝑟𝑛 𝑝𝑟𝑜𝑑\nWe used high burial efficiency under anoxic conditions, consistent with the modeling and experimental\nresults as well as the modern observations that indicate an increase in the burial efficiency of carbon under\nanoxic conditions 42,43. The fraction of oxic condition is calculated as one minus the fraction of anoxia (1-\nf ). We considered full anoxia (f = 1) under very low oxygen condition (pO = 10-13 PAL), consistent\nanoxic anoxic 2\nwith results from geochemical proxies, and modern-like anoxia for modern-like pO (f ~ 0.0001 –\n2 anoxic\n0.001). The fraction of anoxia for pO less than 10-5 PAL was less than 0.05, consistent with the results from\n2\ngeochemical proxies and the biogeochemical modeling that suggest weak ocean oxygenation under the\natmospheric oxygen threshold for production of mass-independent fractionation of sulfur (e.g., 44,45). The\nfraction of anoxia in the ocean for other initial conditions (pO ) was varied between these two end members.\nwas selected in a stochastic manner with the assumption of uniform distribution for the range of each\nparameter.\nThe term R is the ratio of organic carbon to organic phosphorous normalized to modern (oxic) conditions.\nc/p\nThe term R accounts for the effect of change in redox condition and expansion of anoxia on phosphorous\nc/p\nrecycling and organic matter burial. Specifically, under anoxic conditions, the release of phosphorous\nwould be enhanced which can result in an increase in the marine primary production and organic matter\nburial11,46. The link between anoxia and enhanced phosphorous availability however depends on the\navailability of sulfide which controls the availability of dissolved iron (e.g.,47). Given the interaction\nbetween phosphorous and dissolved iron, an increase in sulfide availability would result in a decrease in\nthe iron sink of phosphorous(e.g.,47). Thus, an expansion in the anoxia can foster a higher rate of organic\nmatter burial only if the rate of sulfide production through sulfate reduction is high which requires high\nsulfate levels. We followed the previous studies11,46 to parameterize the term R and included a factor that\nc/p\naccounts for sulfate availability as well:\n𝐶 𝐶\n(\n𝑃\n)𝑎𝑛𝑜𝑥𝑖𝑐∗(\n𝑃\n)\n𝑜𝑥𝑖𝑐\n𝑅 =\n[𝑂2]/[𝑂2]𝑚𝑜𝑑𝑒𝑟𝑛∗(\n𝑃\n𝐶\n)\n𝑎𝑛𝑜𝑥𝑖𝑐\n+(1− [𝑂2]/[𝑂2]𝑚𝑜𝑑𝑒𝑟𝑛)∗(\n𝑃\n𝐶\n)\n𝑜𝑥𝑖𝑐 ∗ 𝑅 (5)\n𝑐/𝑝 𝐶 𝑆𝑂4\n( )\n𝑃 𝑜𝑥𝑖𝑐\n( )\nanoxic oxic\nprevious studies ((C/P) = 1100; (C/P) = 250)11,46. The [O ] and [O ] are, respectively, the\nanoxic anoxic 2 2 modern\nseawater concentration of oxygen at different pO levels and modern conditions. We used the previous earth\n2\nsystem modeling that relate the seawater oxygen concentration to various pO levels48.\nSO4\na Monod-type formulation, similar to the parameterization scheme suggested for the rate of sulfate\nreduction:\n𝑅 =\n[𝑆𝑂4]\n(6)\n𝑆𝑂4\n[𝑆𝑂4]+𝑘𝑚,𝑆𝑂4\nWhere [SO ] and k denote the seawater sulfate concentration and the half-saturation constant for\n4 m,SO4\nsulfate reduction. We assumed sub-mM levels for seawater sulfate under low oxygen conditions (pO <10-2\n2\nrange of seawater sulfate suggested by the previous studies under low pO and high pO levels.”\n2 2\nModified text 4 (Methods):\n“The oceanic sink of oxygen was a function of the strength of marine net primary production and fraction\nof anoxia, as well as the seawater oxygen concentration (e.g., 28,46,51).\n𝐽 = 𝑘1∗ (1−𝜀)∗𝑁𝑃𝑃∗ (𝑓 ∗𝑎𝑟𝑒𝑎)∗(1−\n[𝑂2]\n)\n𝑜𝑐𝑒𝑎𝑛,𝑠𝑖𝑛𝑘 𝑎𝑛𝑜𝑥𝑖𝑐\n[𝑂2]+𝑘𝑚,𝑂2\n(16)\nWhere [O ] and k denote the seawater oxygen concentration and the half-saturation constant for aerobic\n2 m,O2\nrespiration. As explained above, we used the results from previous earth system modeling to relate the pO\n2\nto seawater oxygen concentrations. The constant “k1” is a factor that accounts for variability in the pool\nof reducing compounds (e.g., sulfide, iron, methane) in the ocean. The constant “k1” reflects the relative\ncomparison between the rate of oxygen consumption by reducing species. For instance, a value of 0.1 means\nthat the rate of oxygen consumption by reducing species is about an order of magnitude lower than the\nmaximum rate of anaerobic respiration. Under anoxic and low-oxygen conditions of Archean and\nProterozoic, this value would have likely been low. For instance, most of the sulfide resulting from sulfate\nreduction would have precipitated as pyrite, consistent with the results from sulfur isotope records that\nsuggest a near-one value for the fraction of pyrite (f ). Our results for the J for the modern-like\npyr ocean,sink\ncondition are consistent with the rate of oxygen uptake in the modern oxygen minimum zone. Finally, to\naccount for the uncertainty, this value along with other values used in this equation and elsewhere was\nvaried within a wide range in our stochastic simulation (Table S1).”\nMinor points:\n- I would suggest modifying the title a little to “large-scale dynamics…”\nWe have now modified the title to:\n“ Insights from a dynamical analysis into the history of atmospheric oxygenation”\n5\n- Figure 1: states V and IV are labelled the wrong way around it seems?\nCorrected. Thanks.\nReviewer #3 (Remarks to the Author):\nReviewer #3 (Remarks on code availability):\nI did not run the code myself, but I read the scripts.\nAs noted in my main review file, my main concern with this revision of the manuscript is that there\nappear to be inconsistencies between the code file and the description of the model in the text. In\nparticular, it doesn't look like the parameter values are sampled randomly with a uniform\ndistribution. I may have made a mistake in reading the code, but recommend that the authors\ncheck this and make any necessary changes before publication.\nThanks for the feedback. We have checked the code and modified it accordingly to ensure a random\nsampling with a uniform distribution.\nThe revisions made have greatly improved the manuscript. I have a few further comments and suggestions,\ndescribed below. Thank you very much again for an interesting manuscript to review.\nThanks again for the feedback. You may find below our response to the raised points.\nComments and suggestions\nL79: “following the previous studies” – I suggest removing “the.”\nCorrected. Thanks.\nL88: Repeated sentence here.\nCorrected. Thanks.\nL107—108: “will be sampled randomly within the range suggested in the literature (Table S1), assuming a\nuniform distribution.” It is great that the code is now accessible. However, now that I have been able to\nlook at the file ‘OxygenDynamic.m,’ I am a little concerned about the way the parameters are chosen. I’m\nnot particularly fluent in Matlab so may have made a mistake, but it looks like some parameters are not\nsampled randomly with a uniform distribution, as claimed in the text. For example, the code samples pyrite\nweathering factors from the array ‘Weathering_pyrite_MC,’ which has values of [1500 1700 2000 2200\n2500], which are not evenly distributed. Similarly, ‘pCO2_1_MC’ is [100 200 280 300 500 600 700 1000].\nSee also ‘Temp_MC.’ Apologies if I have misunderstood the code or am looking at the wrong file, but I\nthink either the models should be run with uniform distributions of these variables, or the non-uniform\ndistribution should be described and justified in the text. This method of sampling could affect where the\nmedian line lies and therefore where the steady states lie. See also L352.\nFw_oc comes from randomly sampling values from the list [1 1.5 2 3 3.5 4 4.5 5 7 10], but I think this is\nthe same as Forg in Table S1, which is said to be randomly sampled from a uniform distribution between\nthe values of 3.75 and 7.75. Again, apologies if I’ve misunderstood.\nThanks for this point. We have modified the code to ensure that the parameters would be sampled randomly\nassuming a uniform distribution. In the case of the variable for oxidative weathering of organic (Fw_oc),\nto ensure that in our stochastic simulation, the full range of values would be covered we chose an increment\nunit of one. This choice of increment does not impact our result but would rather help with ensuring a\nsampling for a full range of values in the stochastic simulation.\nL177—207: Two paragraphs from “Our results from the sensitivity analysis suggest…” to “…depending\non the directionality of the change.” The addition of discussion of these results in the main text is\nimportant and very interesting. However, I have a few suggestions about this section:\n• L177: A sentence or two introducing the sensitivity analysis and describing what was done should be\ngiven here, as it hasn’t been mentioned earlier in the paper. This could begin with the first couple of\nsentences from the section in the Methods (see comment for L396-506 below).\n• L177-182: “the rate of hydrogen escape and the atmospheric oxidation of methane have the largest\ninfluence on the emergence of the first and second steady-state for pO2 under anoxic conditions” and\n“our results suggest that the first steady-state point is controlled by the balance between the sink of\noxygen from reaction with reduced volcanic gases and the source of oxygen through the escape of the\nhydrogen to the atmosphere.” Also L197-199: “Of the two steady-state points under low-oxygen\nconditions (states III & IV), the first one is the result of an increase in oxygen production in the wake of\nthe evolution of oxygenic photosynthesis while the rate of oxidative weathering is not high.” It’s not\nclear to me how these statements are obtained from the results presented in Figure S2. More\ndetails about what Figure S2 shows and how the results support these three statements would\nstrengthen the text here and help the reader follow the argument from the figures and results to\nthese implications.\n• “An increase in the rate of hydrogen escape would lead to the disappearance of the first steady-state\npoint (state I).” I don’t follow how this is seen from the results in Fig. S2, so suggest more explanation is\ngiven here to help the reader.\nL179: “…consistent with the previous study.” The previous study in mind should be referenced here.\nThanks for your comments on this paragraph. We have now revised this paragraph as follows:\n“To further explore the influence of each feedback on the global oxygen mass balance, we conducted a\nsensitivity analysis. Specifically, to test the sensitivity of our results for the steady-state pO levels to the\n2\nstrength of each feedback, we fixed the value of each oxygen feedback to its maximum value used in the\nstochastic analysis and recorded the change in the median of dpO /dt vs pO (Figure S2). The magnitude of\n2 2\nchange in the value for a steady-state point for pO in response to a change in the rate of oxygen production\n2\nor consumption indicates the relative importance of that feedback in driving the steady-state point for pO .\nin the value of the first steady-state point for pO with no substantial change on other steady-state points\n2\nfor pO . This result suggests an important role of hydrogen escape in driving the first pO steady-state point.\n2 2\nAccording to the results from the sensitivity analysis, the rate of hydrogen escape and the atmospheric\noxidation of methane have the largest influence on the emergence of the first and second steady-state for\npO under anoxic conditions (Figure S2), consistent with the previous study31. Specifically, a change in the\n2\nrate of hydrogen escape and atmospheric methane oxidation can result in a shift or disappearance of the\nfirst and second steady-state points for atmospheric oxygen levels (Figure S2). The major role of the rate\n8\nof hydrogen escape on the first stable state for the atmospheric oxygen level (state I) indicates the key role\nof hydrogen escape in stabilizing the atmospheric oxygen level before the rise of oxygenic photosynthesis.\nAccording to our results, an increase in the rate of atmospheric methane oxidation influences the emergence\nof the second steady-state point (state II) around 10-5 PAL, consistent with the modeling results from\nprevious studies (e.g.,31).”\nMethods:\nL347: Jvolc is defined as volcanic outgassing in L81 of the main text, but it would be useful if it was also\ndefined here with the other O2 sinks.\nCorrected. Thanks.\nL419: RH2 is not defined. Is it the same as Jhydrogen? See also L445, where it seems to denote the rate of\nmethane oxidation.\nYes, we have now corrected them in the equations. Thanks.\nL354—494: There is no section detailing how J (in Equation 1) was parameterized.\nvolc\nWe used a range of values for J that are now presented in Table S1.\nvolc\nL396—506:\n• This paragraph is a good addition to the Methods. I suggest that the content of this paragraph\n(particularly the first two sentences) is put in the main text instead. The first two sentences could\ngo at the beginning of the paragraph at L177 before “Our sensitivity analysis…”\n• L497-499: This is a helpful sentence, which should be in the main text as mentioned above. It would\nbe helpful to add whether 1000 runs were done for each case, and if not, how many were run.\nThanks for the suggestion. We have now moved this paragraph into the main text. The updated paragraph\nin the main text now reads as:\n“To further explore the influence of each feedback on the global oxygen mass balance, we conducted a\nsensitivity analysis. Specifically, to test the sensitivity of our results for the steady-state pO levels to the\n2\nstrength of each feedback, we fixed the value of each oxygen feedback to its maximum value used in the\nstochastic analysis and recorded the change in the median of dpO /dt vs pO (Figure S2). The magnitude of\n2 2\nchange in the value for a steady-state point for pO in response to a change in the rate of oxygen production\n2\nor consumption indicates the relative importance of that feedback in driving the steady-state point for pO .\nin the value of the first steady-state point for pO with no substantial change on other steady-state points\n2\nfor pO . This result suggests an important role of hydrogen escape in driving the first pO steady-state point.\n2 2\nAccording to the results from the sensitivity analysis, the rate of hydrogen escape and the atmospheric\noxidation of methane have the largest influence on the emergence of the first and second steady-state for\npO under anoxic conditions (Figure S2), consistent with the previous study31. Specifically, a change in the\n2\nrate of hydrogen escape and atmospheric methane oxidation can result in a shift or disappearance of the\nfirst and second steady-state points for atmospheric oxygen levels (Figure S2). The major role of the rate\n8\nof hydrogen escape on the first stable state for the atmospheric oxygen level (state I) indicates the key role\nof hydrogen escape in stabilizing the atmospheric oxygen level before the rise of oxygenic photosynthesis.\nAccording to our results, an increase in the rate of atmospheric methane oxidation influences the emergence\nof the second steady-state point (state II) around 10-5 PAL, consistent with the modeling results from\nprevious studies (e.g.,31).”\nFigures:\nFigure 1b: The caption for this figure has been greatly improved. However, a description of what the\ncircular points represent is missing.\nThanks. The caption is now modified as:\n“Figure 1. The results of stochastic dynamical analysis of global oxygen mass balance. Panel (a) shows\nthe median of the results of the stochastic analysis for a wide range of environmental and initial conditions.\nPanel (b) corresponds to the full range of results for different pO steady states (dpO /dt = 0). The open\n2 2\ncircles in panel (b) represent points outside the range for each steady-state. Our results suggest five steady states\nfor the atmospheric oxygen level. The box plot in panel (b) is the result of the stochastic\nsimulation where the middle line in the box corresponds to the median value, the edges of the box\ncorrespond to 25th(lower edge) and 75th (upper edge) percentile, and the ends of the whiskers correspond\nto the minimum and maximum values for each state. The results of the stochastic simulation were filtered\nto be consistent with the results from geochemical proxies. We only considered results that include at least\nthree steady-state points: a modern-like steady state, a steady state under anoxic condition, consistent with\ngeochemical proxies that indicate an anoxic atmosphere (e.g., mass-independent fractionation of sulfur\n(MIF-S)) before the GOE, and one steady-state point under the low-oxygen condition that is largely\nreflective of an environmental conditions that is neither anoxic (lack of MIF-S signal) and nor is modern-\nlike (~ 1 present atmospheric oxygen level (PAL)).”\nSI:\nFigure S1 caption: It would be helpful to reference the “previous photochemical modeling” explicitly here,\neither by using the numbered reference(s) or an in-text citation (e.g., “The circles represent data from Name\net al. (20XX)”).\nThanks. We included a citation in the caption. The revised caption now reads as:\n“Figure S1. Relationship between atmospheric oxygen concentration and the oxidation parameter for\nmethane oxidation. The circles represent data from the previous photochemical modeling13, and the blue\nline corresponds to the best curve fit.”\nFigure S2: A sentence describing the blue line could be added to the caption for clarity, e.g., “The blue\nline shows the…”\nThanks. The caption is now modified as:\n“Figure S2. The results from the sensitivity analysis of the pO steady-state point to the change in the\n2\nstrength of various oxygen feedback. The blue line represents the median line of the result of the stochastic\nanalysis for a wide range of environmental and initial conditions.”\neither defined or changed.\nThanks. Yes, the y-axis is Jocean,sink.\nGeneral comments:\nSome grammatical errors remain, which can obscure the clarity of the argument.\nThanks. We read through the text to avoid any grammatical errors. Thanks again for the feedback and for\nthe detailed comments.\nReferences:\nAlcott, L.J., Mills, B.J.W., and Poulton, S.W., 2019, Stepwise Earth oxygenation is an inherent property\nof global biogeochemical cycling: Science, v. 366, p. 1333–1337", "reviewer_doi": "10.1016/j.gca.2019.02.011.\nReviewer #2 (Remarks to the Author):\nI am glad the authors have finally provided a clear explanation for the ~0.1 PAL O2 instability in their\nmodel. The mechanism is already published by Laakso and Schrag in 2019, who are cited in the\nnew description. The authors have also now checked the rates of marine oxygen consumption at\nthe present day. So I have no further technical issues with the paper.", "reviewer_switzerland": "A low sulfate euxinic ocean analogue:\nGeochimica et Cosmochimica Acta, v. 251, p. 116–135", "reviewer_decision_letter": "Decision Letter\n\nDear Dr. Fakhraee,\n\nYour manuscript \"Insights from a dynamical analysis into the history of atmospheric oxygenation\" has been peer-reviewed and I am pleased to inform you that it is now accepted for publication in Nature Communications.\n\nThank you for choosing Nature Communications for your work."}, "author_response": ""}]

Decision letter