Photochemical modelling of atmospheric oxygen levels confirms two stable states

Abstract Various proxies and numerical models have been used to constrain O2 levels over geological time, but considerable uncertainty remains. Previous investigations using 1-D photochemical models have predicted how O3 concentrations vary with assumed ground-level O2 concentrations, and indicate how the ozone layer might have developed over Earth history. These classic models have utilised the numerical simplification of fixed mixing ratio boundary conditions. Critically, this modelling assumption requires verification that predicted fluxes of biogenic and volcanic gases are realistic, but also that the resulting steady states are in fact stable equilibrium solutions against trivial changes in flux. Here, we use a 1-D photochemical model with fixed flux boundary conditions to simulate the effects on O3 and O2 concentrations as O2 (and CH4) fluxes are systematically varied. Our results suggest that stable equilibrium solutions exist for trace- and high-O2/O3 cases, separated by a region of instability. In particular, the model produces few stable solutions with ground O2 mixing ratios between 6 × 10 − 7 and 2 × 10 − 3 ( 3 × 10 − 6 and 1% of present atmospheric levels). A fully UV-shielding ozone layer only exists in the high-O2 states. Our atmospheric modelling supports prior work suggesting a rapid bimodal transition between reducing and oxidising conditions and proposes Proterozoic oxygen levels higher than some recent proxies suggest. We show that the boundary conditions of photochemical models matter, and should be chosen and explained with care.


Introduction
Improved constraints of atmospheric oxygen levels over Earth history are important for an enriched understanding of how life and Earth have co-evolved. Molecular oxygen (O 2 ) currently makes up 21% of the atmosphere, but this has not always been the case. The disappearance of sulphur isotope mass-independent fractionation (S-MIF) in the geological record at 2.4-2.5 Ga (Bekker et al., 2004;Farquhar et al., 2000;Warke et al., 2020) is considered strong evidence for the Great Oxidation Event (GOE). To explain this, models require the partial pressure of oxygen (pO 2 , also referred to as the O 2 mixing ratio) to be less than 2.1 × 10 −6 (10 −5 times the present atmospheric level (PAL)) for the Archaean and more than 2.1 × 10 −6 after the GOE (Pavlov and Kasting, 2002).
(and its interpretation) in the geological record is debated, with some studies predicting much higher pO 2 between 0.8 and 1.1 Ga (Gilleaudeau et al., 2016) and even further back in the Mesoproterozoic (Canfield et al., 2018). Even for the Phanerozoic, where the continuous presence of animal life indicates relatively high oxygen levels, uncertainty remains (Krause et al., 2018). For the last 400 Myr, a continuous charcoal record suggests oxygen levels between 15% and 30% (Bergman et al., 2004;Glasspool and Scott, 2010;Glasspool et al., 2015).
Previous 1-D photochemical models (Kasting and Donahue, 1980;Kasting et al., 1985;Segura et al., 2003;Zahnle et al., 2006) have investigated potential palaeo-atmospheres by explicitly constructing models with a range of pO 2 . These models utilised 'fixed mixing ratio' lower boundary conditions (LBCs) for O 2 , in which the experimenter sets the pO 2 and the model is allowed to vary fluxes across the lower boundary and other species' concentrations to produce a steady-state solution which satisfies the given boundary conditions. In a brief summary overview, these previous studies (reproduced below) found that the ozone column density (the number of ozone molecules in an atmospheric column with a surface area of one square centimetre) initially increases as a power-law with pO 2 , before saturating at higher O 2 concentrations (Fig. 1a).
The results of these 1-D photochemical modelling studies have been simplified, parametrised and incorporated into a number of Earth system evolution models, which include shorter-term biological and atmospheric feedbacks with longer-term planetary redox fluxes, such as volcanic degassing and hydrogen escape. Intriguingly, these models have predicted bimodal behaviour with respect to oxygen concentrations. For example, the 2-box model of Goldblatt et al. (2006) predicted two regions only: pO 2 < 2.1 × 10 −6 and pO 2 > 10 −3 . Similarly, the 3-box biogeochemical model of Laakso and Schrag (2017) found stable solutions with pO 2 < 2.1 × 10 −8 and pO 2 > 2.1 × 10 −4 . Other biogeochemical models predict a 'GOE' in which O 2 mixing ratios change from 10 −6 to 2 × 10 −3  or 10 −2 (Alcott et al., 2019) relatively quickly. While biological re-organisation to oxic conditions provides a positive feedback (Catling et al., 2007;Daines and Lenton, 2016), the primary driver captured is the development of an ozone layer which shields O 2 from photolysis Goldblatt et al., 2006), as captured by parametrisations of 1-D photochemical models.
Here, we revise the classic photochemical modelling calculations which relate O 2 and O 3 , by utilising 'fixed flux' LBCs for O 2 and other gases. The O 2 flux into the atmosphere is specified by the experimenter, and the model predicts the resulting mixing ratio profiles in equilibrium with the radiation field, kinetics, and other physical processes (e.g. lightning, particle formation, condensation). Flux boundary conditions represent a closer conceptual match to what real atmospheres do (e.g. fluxes vary and concentrations adjust), compared to the assumption that the planetary and biogenic fluxes somehow adjust in order to maintain fixed ground-level concentrations. First, we briefly introduce the 1-D photochemical model, Atmos. In Section 3, we use Atmos with fixed mixing ratio LBCs for O 2 to replicate previous results, but observe that these models predict potentially unrealistic fluxes. Our primary results comprise a systematic study of atmospheric chemistry resulting from variable O 2 and methane fluxes across the lower boundary. Specifically, we i) vary O 2 fluxes at several fixed CH 4 :O 2 flux ratios (Section 4.1); ii) vary CH 4 :O 2 flux ratios at several fixed O 2 fluxes (Section 4.2); and explore sensitivity to iii) oxidative weathering (Section 4.4) and iv) other important redox fluxes (Section 4.6).
Our investigation of 2067 flux-driven model atmospheres reveals a bimodal oxygen distribution similar to those observed in biogeochemical models, but one that arises entirely in the atmo- Circles, triangles, squares and diamonds show results using Atmos with the same boundary conditions as previous authors. '1980LBCs,' '1985LBCs,' '2003LBCs' and '2006 LBCs' refer to the boundary conditions used by Kasting and Donahue (1980), Kasting et al. (1985), Segura et al. (2003) and Zahnle et al. (2006), respectively. Crosses show Case 0 results, with fixed mixing ratio boundary conditions for O 2 and fixed flux boundary conditions for CH 4 , H 2 , CO and N 2 O. For the boundary conditions used, see Table 1. (a) Ozone column densities with user-specified ground-level O 2 mixing ratio. The results from previous studies are also shown (orange lines). (b) Predicted oxygen fluxes, and (c) predicted CH 4 :O 2 flux ratios required to produce the fixed mixing ratio model atmospheres, plotted against user-specified O 2 mixing ratio. (For interpretation of the colours in the figures, the reader is referred to the web version of this article.) sphere itself. We explicitly demonstrate that some classic results are unstable equilibrium solutions (Section 4.3) and argue that many intermediate oxygen concentrations are photochemically unstable. In Section 5, we discuss implications for palaeo-O 2 con-straints, provide suggestions for future work in both photochemical and biogeochemical modelling, and argue that the choice of numerical methods used by 1-D photochemical models has implications for how the Earth system science community considers the evolution of atmospheric oxygen.

Model description
We use the one-dimensional photochemical model Atmos (open source code available via https://github .com /VirtualPlanetaryLabora tory /atmos), most recently described in Arney et al. (2016). The 'ModernEarth+Cl' template utilised in this study incorporates 87 atmospheric species involved in 372 reactions, whose rates have been updated to follow recent recommendations (Table D1). An altitude grid of 160 0.5 km layers gives a simulated atmosphere of height 80 km, with the tropopause set at 11 km, suitable for a mid-latitude profile. Each model presented here was run with modern eddy diffusion and temperature profiles (Appendix E), and we analyse only steady-state solutions. The lower boundary of the model represents Earth's surface, across which we prescribe fluxes of volcanic and biogenic gases into the bottom-most model grid layer, in addition to gaseous, aqueous and particulate deposition out of the model atmosphere. LBCs are specified for each atmospheric species (Tables B1 and C1) and discussed further below.

Results from fixed mixing ratio photochemical modelling
We used fixed mixing ratio LBCs for O 2 to compare outputs from Atmos to the results of previous studies (shown in Fig. 1a). Kasting and Donahue (1980), Kasting et al. (1985) and Segura et al. (2003) used fixed flux or deposition velocity (specifying the removal rate from the atmosphere) LBCs for CH 4 , H 2 , CO and N 2 O. The magnitudes of these fluxes were determined by predicted fluxes from a fixed mixing ratio LBC model simulating the modern Earth for these species (with the exception of Kasting and Donahue's (1980) H 2 and CO, which were allocated zero flux LBCs). Resulting negative (i.e. out of the atmosphere into the ocean) H 2 fluxes led Kasting et al. (1985) and Segura et al. (2003) to use deposition velocity and negative fixed flux boundary conditions for H 2 , respectively. Zahnle et al. (2006) used fixed mixing ratio LBCs for CH 4 , and chose large positive fluxes of H 2 and CO. These (and all LBCs used in this study) are summarised in Table 1. We specify units of flux in 'photochemical units' or 'pu' (1 pu = 1 molecule cm −2 s −1 ∼ 2.7 × 10 −10 Tmol yr −1 ).
Using Atmos, we produced suites of steady-state model atmospheres with O 2 ground-level mixing ratios varying between 2.1 × 10 −11 and 0.42, choosing LBCs using the specific methods of the previous four studies. A fifth suite (hereafter referred to as Case 0; see Table 1) also incorporates fixed flux LBCs for CH 4 , CO and N 2 O equal to those computed by modern fixed mixing ratio calculations. Like the models of Kasting et al. (1985) and Segura et al. (2003), Atmos predicts a net flux of H 2 out of the atmosphere when run solely with modern fixed mixing ratio LBCs. This is physically unrealistic and numerically allows models significant leeway to mask redox imbalances by dumping hydrogen into the (assumed) ocean. We incorporated the predicted fluxes as fixed flux LBCs and tuned them, primarily by slightly reducing the larger (biogenic) CH 4 flux, to achieve an atmosphere predicting modern mixing ratios for the five species, but with a net flux of H 2 into the atmosphere within an order of magnitude of the modern volcanic flux (Aiuppa et al., 2011). A fixed mixing ratio LBC for preindustrial CO 2 (280 ppmv) was used. Allowing the model to compute CO 2 fluxes is justifiable here as CO 2 is redox-neutral, whereas O 2 , CH 4 , H 2 , CO and N 2 O are not .
The ozone column density increases with specified ground-level pO 2 but plateaus above a pO 2 of 10 −2 , following a chemical pattern discussed in previous work (Fig. 1a). Differences between our models for trace-O 2 concentrations are primarily due to the different choices of boundary conditions. Our results are very similar to those of previous models, with small changes due to updated kinetic rate coefficients, which is unsurprising given that Atmos shares a common code heritage. However, in what follows we discuss observations that led us to discover that some of the steadystate model atmospheres predicted by fixed mixing ratio LBCs exhibit behaviour of unstable equilibrium solutions.

The restrictions of a fixed mixing ratio boundary condition
By construction, fixed mixing ratio LBCs enable the model to vary lower boundary fluxes to produce the user-specified concentration, but the required fluxes may not be physically realistic. Fig. 1b shows the O 2 fluxes computed by the model to maintain the specified O 2 mixing ratio boundary conditions in Fig. 1a. These fluxes (and others not shown) vary across the experiments, making the direct comparison of simulated atmospheres challenging, and are not necessarily representative of how a real biosphere would function. For example, Fig. 1c shows that though the predicted CH 4 :O 2 flux ratio is ∼0.1 for most high-O 2 atmospheres, consistent with theoretical estimates and modern measurements, the CH 4 :O 2 flux ratio varies significantly for moderate-O 2 atmospheres, and is greater than 0.5 for low-O 2 Archaean-like atmospheres, which exceeds the stoichiometric limit that life can extract from CO 2 and H 2 O substrates. Model-computed imbalances in CO and H 2 fluxes into and out of the model (not shown) further complicate and occasionally mask potential non-physical behaviour in other species. A more intuitive approach is to simulate volcanism and biology as 'fluxes' to which the composition of the atmosphere adjusts, which motivates our use of fixed flux LBCs.

Results from flux-driven photochemical modelling
We start with a 'flux-driven modern model,' which we construct as the best-fit model reproducing modern conditions using fixed flux LBCs for O 2 , along with our 'Case 0' fixed flux LBCs for CH 4 , H 2 , CO and N 2 O ( Table 1). The model-predicted O 2 flux is 8.9 × 10 11 pu, resulting in a CH 4 :O 2 flux ratio of 0.094. While the modern gross O 2 production is much larger (∼4 × 10 13 pu), our O 2 flux ignores very short-term sources and sinks of atmospheric O 2 , while capturing the seasonal imbalance between O 2 production and destruction as well as the most important components of the short/medium-term carbon cycle . Our methane flux (8.4 × 10 10 pu) is similar to the estimated net modern methane flux to the atmosphere from natural sources, after a large proportion of the gross methane production flux has been oxidised in sediments (Schlesinger and Bernhardt, 2013, p. 434).
We produced suites of model atmospheres to explore a range of non-modern biologic flux conditions, primarily by changing the absolute O 2 flux and CH 4 :O 2 flux ratio. We cover a parameter space designed to simulate a range of possible states of primary productivity, methane production, methanotrophy and sulphate reduction. Specifically, we varied the O 2 flux between ∼0.01 to ∼40 times that of the flux-driven modern model, which we conceptually identify with varying levels of primary productivity. We assigned the CH 4 flux at various flux ratios from 0.094 (the CH 4 :O 2 flux ratio in the 'flux-driven modern model') to 0.5 (the stoichiometrically-balanced ratio). The former is representative of the modern Earth system, and the latter is a commonly-assumed simplification for the late-Archaean biosphere (Catling and Claire, 2005;Daines and Lenton, 2016). Intermediate CH 4 :O 2 flux ratios might represent biospheres with increasing anaerobic oxidation of Table 1 Lower boundary conditions for the models described in this study. Units of flux are in photochemical units (1 pu = 1 molecule cm −2 s −1 ∼ 2.7 × 10 −10 Tmol yr −1 ).

Case 1: varying O 2 fluxes at fixed CH 4 :O 2 ratios
For a first test case, we varied O 2 flux between 10 10 pu and 4 × 10 13 pu. We produced 5 suites of up to 372 model atmo- patterns to pO 2 , with clusters of solutions at low and high column densities without intermediates (Fig. 2c).
When a CH 4 :O 2 flux ratio of 0.5 was used (results not shown in Fig. 2), the greater total flux of reductants (including H 2 , CO and N 2 O) than oxidants into the atmosphere resulted in a pO 2 always less than the CH 4 mixing ratio. For the 0.49 CH 4 :O 2 flux ratio suite (Fig. 2a-c), the O 2 flux was eventually increased sufficiently such that it exceeded the total reductant flux, despite co-increasing methane fluxes. For this suite, there are 'Archaean-like' solutions with O 2 concentrations up to 10 −6 and one steady-state atmosphere with an O 2 mixing ratio of 4 × 10 −5 , followed by a series of solutions with pO 2 exceeding 3 × 10 −3 (1.5% PAL).
To summarise, our 1222 'Case 1' models predict that oxygen, methane and ozone levels exist in one of two states -trace-O 2 /O 3 solutions, and high-O 2 /O 3 solutions -separated by a 'gap' in which there is only one (very specific) solution for intermediate O 2 surface mixing ratios. To better understand, we explored the detailed altitude-dependent chemistry of models spanning these states.

Two stable states of atmospheric chemistry?
We present mixing ratio profiles (Fig. 2d) and ground-level UV fluxes (Fig. 2e). Panels (i)-(iv) in Fig. 2d show O 2 , O 3 and CH 4 mix- In the 'high-O 2 /O 3 ' atmospheres (panels (iii) and (iv) in Fig. 2d), O 2 is well-mixed and more abundant than methane. In panel (iii), a modern-like stratospheric ozone layer is emerging, which is fully developed by panel (iv). The significant ozone shielding results in increased methane mixing ratios and low UV fluxes at the surface, with modern-like behaviour where no photons shortward of ∼290 nm reach the ground (Fig. 2e). The other high-O 2 atmosphere (panel (iii) in Fig. 2d) blocks substantial UV within the 200-320 nm range, but some photons at DNA-damaging wavelengths (200-220 nm) will reach the ground. This case is similar to the Segura et al.
(2003) (fixed mixing ratio boundary condition) model of 2.1 ×10 −3 for O 2 , which was argued to be the tipping point at which ozone shielding significantly impacts ground UV fluxes.
The panel (v) atmosphere is similar to those in panels (i) and (ii), while panel (vi) shows an atmosphere with tropospheric O 2 and CH 4 mixing ratios both approximately 80 ppm -on the cusp of the reducing/oxic divide. However, there is not sufficient ozone to provide tropospheric UV shielding of O 2 (Fig. 2e). While being a steady-state solution, it seems unlikely that this kind of atmosphere would be stable against small perturbations to the oxygen flux. The dearth of stable solutions between the trace-O 2 and high-O 2 attractors appears primarily due to positive feedbacks between increasing O 2 and O 3 , which result in decreasing O 2 photolysis rates Goldblatt et al., 2006). O 2 photolysis occurs only at wavelengths of 242 nm and below. The most abundant and effective absorbing molecule in the 200-300 nm wavelength range is O 3 , with peak absorption at 240 nm. Increasing O 3 concentrations (Fig. 2c) therefore shield O 2 in the lower atmosphere from photolytic destruction ( Figure A1), enabling O 2 concentrations to rise for a given flux. We explore this further in Section 4.3.

Case 2: varying CH 4 :O 2 ratios at fixed O 2 fluxes
A second set of experiments designed to expand the parameter space of modelled fluxes investigates the reproducibility of our results under different conceptualisations of biogenic flux evolution. In the Case 2 suites, we fixed the absolute O 2 flux at five values (10 11 , 3 × 10 11 , 5 × 10 11 , 10 12 and 5 × 10 12 pu), while decreasing the CH 4 :O 2 flux ratio systematically from 0.5 to 0.1. These experiments might roughly capture the effects of increasingly sulphaterich waters driving the locus of methane oxidation into sediments at a given primary productivity, but are admittedly a vast simplification of biospheric processes . The computed O 2 and CH 4 mixing ratios and corresponding O 3 column densities are shown in Fig. 3. As the CH 4 :O 2 flux ratio is decreased from 0.5 to a critical value dependent on the absolute O 2 flux, an abrupt transition from reducing to oxidising atmospheres occurs (other than for an O 2 flux of 10 11 pu, which is sufficiently low that this critical value is not reached for any CH 4 :O 2 flux ratio less than 0.1). The Case 2 experiments demonstrate a wider 'gap' in solutions than Case 1, with no predicted model atmospheres with O 2 mixing ratios between 4 × 10 −7 and 10 −3 .

The photochemical (in)stability of atmospheres 'in the gap'
Atmospheres with O 2 mixing ratios between 4 × 10 −7 and 4 × 10 −4 are not produced in the flux-driven experiments described so far in this study, barring one case highlighted. By contrast, the identical numerical model utilising fixed mixing ratio boundary conditions produces results for all surface O 2 mixing ratios between 2.1 × 10 −11 and 0.42 (Fig. 1). As all results in this and previous efforts are steady-state equilibrium solutions, it should be possible, numerically speaking, to carefully specify fixed flux boundary conditions to reproduce model atmospheres with any ground-level O 2 mixing ratio, even within the 'gaps' identified in Cases 1 and 2. To explore this, we chose three of our Case 0 fixed mixing ratio LBC models, with pO 2 values of 2.1 × 10 −6 , 2.1 × 10 −5 and 2.1 × 10 −4 . We entered the predicted fluxes of O 2 (to 12 significant figures) as fixed flux boundary conditions and used the predicted species densities as initial values.
For each case, the flux-driven models did not converge when using the fluxes predicted by the fixed mixing ratio models. As fluxes were perturbed incrementally (at constant CH 4 :O 2 ), no mod- els converged until the O 2 fluxes were ∼6%, ∼6% and ∼3% higher, or ∼1%, ∼1% and ∼4% lower than Case 0, respectively (Fig. 4). Subsequently, model solutions with 3 × 10 −8 < pO 2 < 5 × 10 −4 were not produced. Rounding our input fluxes to three significant figures for the 2.1 × 10 −5 case to test for numerical precision issues made no significant difference to the results. These tests supplement the preliminary conclusions from Cases 1 and 2, demonstrating that some steady-state model solutions computed using fixed O 2 mixing ratio LBCs between 4 × 10 −7 and 5 × 10 −4 are in fact unstable equilibrium solutions. They adjust to new solutions given trivial changes in lower boundary fluxes. This suggests, but does not fully prove, that fixed mixing ratio LBCs produce unstable equilibrium solutions in some (potentially large) regions of parameter space. A full stability analysis requires time-dependent computations beyond the scope of this effort, but will be forthcoming. Users of 1-D photochemical models (who must already carefully choose boundary conditions (e.g. Domagal-Goldman et al., 2014;Harman et al., 2018) are further cautioned to consider the stability of model results against small variations in flux, especially if using fixed mixing ratio boundary conditions.

Case 3: including negative feedbacks from oxidative weathering
Photochemical models to date have ignored a key negative feedback on oxygen fluxes at Earth's surface: as atmospheric O 2 concentrations build up, the removal rate via oxidative weathering increases. To test whether this could stabilise atmospheres with lower pO 2 , we ran suites of models with the inclusion of deposition velocity boundary conditions that remove O 2 from the lower boundary at a rate proportional to the O 2 mixing ratio. Numerically, this requires the photosynthetic O 2 flux to be delivered into the atmosphere at 0.5 km above the surface, but this has a negligible effect on the vertical profiles because large eddy diffusion coefficients rapidly mix near-surface layers.
Firstly, we defined the O 2 lower boundary flux as F out = V dep [O 2 ]N 1 , where V dep is the deposition velocity (cm s −1 ), N 1 is the atmospheric density in the first grid step, and [O 2 ] is the ground-level O 2 mixing ratio. In a second suite, we used F out = V dep [O 2 ] 0.5 N 1 , following pyrite oxidation kinetics experiments that show a half-power law relation (Johnson et al., 2019). We focus on pyrite rather than organic carbon oxidation (which we consider a redox-neutral part of the reversal of photosynthesis in our conceptualisation), but conveniently carbon oxidation rate laws are also approximately half-power in O 2 (e.g. Daines et al., 2017). We tuned values of V dep to 5 × 10 −9 cm s −1 and 2.3 × 10 −9 cm s −1 , respectively, in order to reproduce the modern oxidative weathering flux of ∼2.6 × 10 10 pu (Holland, 2002) and modern pO 2 , when using the fixed-flux modern model. For a third suite, we increased V dep by an order of magnitude for the half-power law case.
The results (Fig. 5) reveal that the two-state behaviour observed in Cases 1 and 2 persists in the Case 3 suites, despite the inclusion of a negative feedback flux. Increasing V dep by an order of magnitude produces several steady-state model atmospheres with O 2 mixing ratios between 10 −6 and 10 −3 , but they are not stable against trivial changes in flux. Fig. 6 illustrates a key result of this work. Predicted O 2 ground mixing ratios are plotted against O 3 column densities for the 2067 flux-driven steady-state model solutions in Cases 1-3. Given the monotonic increase in O 3 column density with O 2 mixing ratio for the classic (Case 0) result, it previously appeared plausible to imagine that Earth's atmosphere could have existed in a wide range of stable mid-oxic states. In contrast, our flux-driven results suggest that, for a broad range of O 2 and CH 4 fluxes, the majority of solutions cluster in two locations, with limited stable solutions between.

Confirmation of two stable states of atmospheric oxygen chemistry?
The model atmospheres show a bimodal distribution, with 'high-' (mode of 0.1-0.2) and 'trace-' (mode of 1-2 × 10 −8 ) O 2 atmospheres (Fig. 6, upper panel). With 'high-O 2 ' atmospheres defined as those with pO 2 above 2.1 × 10 −6 and, as a conservative measure, discarding the atmospheres with O 2 concentrations greater than 0.30 (since these are unlikely to have existed for most of the Phanerozoic), 95% (2σ ) of our model atmospheres have pO 2 greater than 5 × 10 −3 , and 97.5% have pO 2 greater than 2 × 10 −3 (∼1% PAL). The results may be additionally conservative as the bulk of the solutions computed with lower than 1% PAL O 2 were from cases in which oxidative weathering was set one order of magnitude higher than presently constrained values. With 'trace-O 2 ' atmospheres defined as those with pO 2 less than 2.1 × 10 −3 , 95% of our model atmospheres have pO 2 less than 6 × 10 −7 . While by no means a comprehensive statistical treatment, this supports previous studies Daines and Lenton, 2016;Goldblatt et al., 2006;Laakso and Schrag, 2017), which concluded that there is a region of pO 2 parameter space spanning several orders of magnitude in which model atmospheric solutions are particularly sparse. For Cases 1-3, the absolute fluxes of O 2 and CH 4 at which the atmospheres switch from reducing to oxic varies. Fig. 7 normalises this by plotting against K oxy , which is defined as the O 2 source flux divided by the stoichiometrically-balanced average fluxes of all reducing gases (Catling and Claire, 2005). Model atmospheres with K oxy less than 1 are reducing, and those with K oxy greater than 1 are oxic, suggesting that the ratio of total reductant to oxidant flux determines in which of the two states the atmosphere lies. It further illustrates the very narrow range in flux space for which intermediate solutions beyond our 95% confidence bounds exist. If any of these lower pO 2 solutions were to persist in the Fig. 7. Model-predicted ground-level O 2 mixing ratios plotted against K oxy , for all model atmospheres from Cases 1-3.
Model atmospheres with K oxy < 1 are reducing, while those with K oxy > 1 are oxidising.
Earth system, they would require very strong (unspecified) negative feedbacks to maintain biologic fluxes in a tight window.

Exploring flux boundary conditions for other redox-relevant species
As a final set of tests, we examined the sensitivity of model output to the choice of boundary conditions for the reducing gases H 2 and CO. The H 2 flux has been shown to affect the pO 2 of steadystate model atmospheres (Laakso and Schrag, 2017). The two-state behaviour described above persists despite changes in the magnitude of reducing fluxes, and further strengthens our conclusion that the user-specification of boundary conditions must be carefully considered, even for trace species.

Case 4: magnitudes of reducing fluxes
For a fourth test case, we varied the magnitude of the fluxes of H 2 , CO, H 2 S and SO 2 into the model atmosphere, using the 'modern low,' 'modern high' and 'Archaean high' values from  ; Table 1). When we used the 'modern low' fluxes as fixed flux LBCs for a model with fixed mixing ratio (modern-like) boundary conditions for O 2 , CH 4 and N 2 O, the resulting mixing ratios for H 2 and CO were 7.7 × 10 −7 and 2.5 × 10 −8 , which are slightly higher and lower than modern estimates, respectively (Ehhalt et al., 2001). Using the computed N 2 O flux and fixed CH 4 :O 2 flux ratio of 0.16, we produced a model suite with varying O 2 flux. We repeated this method for the 'modern high' and 'Archaean high' fluxes.
For all three scenarios, a jump from reducing to oxic atmospheric solutions within a narrow flux range was observed (Fig. 8).
The 'gap' in O 2 mixing ratio solutions is slightly different to that produced in Cases 1-3, and again varies significantly over a range of O 2 fluxes. Regardless, the absolute magnitude of H 2 and CO fluxes (between end-member parametrisations) does not affect the overall two-state behaviour of the model atmospheres with varying user-specified O 2 flux.

Case 5: additional drawdown fluxes of reductants
Both H 2 and CO are consumed voraciously by microbes. We tested the model sensitivity to additional drawdown fluxes, by implementing a linear feedback flux equal to V dep [ X ]N 1 , where [ X ] is the H 2 or CO mixing ratio. Following Zahnle et al. (2006), we chose deposition velocity values of 2.5 × 10 −4 and 1.2 × 10 −4 after Kharecha et al. (2005). We distributed the upward fluxes (from  Zahnle et al. (2006), for H 2 , CO, H 2 S and SO 2 , respectively. White diamonds show the Case 5 results, for which H 2 and CO are given additional deposition velocities (see Table 1).
the flux-driven modern model) log-normally through the first 10 km of the atmosphere . Using the mixing ratio of CO in the flux-driven modern model, we calculated the resulting downward flux and adjusted the upward flux such that the net flux remained the same as for Cases 1-3. Using the same method for H 2 resulted in a net downward flux, so we increased the upward flux until the net H 2 flux was positive (i.e. into the atmosphere). The resulting CO and H 2 mixing ratios were fairly similar to the global measured averages (6.2 × 10 −8 and 5.5 × 10 −7 , respectively; Ehhalt et al., 2001). Using these combined flux and deposition velocity conditions for H 2 and CO, we reproduced a suite of models with varying O 2 flux, at a constant CH 4 :O 2 flux ratio (0.094). The results (Case 5 on Fig. 8) are our best attempt to properly model CO and H 2 , and also display the two-state behaviour observed in Cases 1-3.

Case 6: mixing ratio boundary conditions for H 2 and CO -a cautionary tale
Using modern-like fixed mixing ratio LBCs for H 2 , CO and N 2 O produced results broadly similar to those described in previous subsections, in that there were reducing and oxic model atmospheres, and a transition between them within a very narrow flux range. However, this model setup (preliminary work presented by Gregory et al., 2019) also produced a third cluster of 'very-low O 2 ' atmospheres (7 × 10 −7 < pO 2 < 3 × 10 −5 ; not shown), occurring only inside the narrow transition flux range, and featuring considerable ozone, but not a fully protective ozone shield. In addition, the 'high-O 2 ' atmospheres only had minimum O 2 mixing ratios of 10 −2 , considerably affecting the application of our model to palaeo-O 2 levels over Earth history.
This demonstrates that the fixed mixing ratio LBCs of some trace gases enables the same numerical issues (potentially unphysical fluxes across the lower boundary) as for the more dominant species. Our results from Cases 1-5 with fixed flux boundary conditions for H 2 and CO allowed us to determine that the original three-state behaviour was driven by instabilities in the modelpredicted lower boundary fluxes of H 2 and CO, which the model was switching between net input and output depending on the CH 4 :O 2 flux ratio. This emphasises our primary conclusion that any photochemical modelling effort must fully describe choices for boundary conditions, even for trace species, in order to allow reproducibility and assess reliability.
Having determined from Cases 4-6 that H 2 and CO flux boundary conditions produce rather different model atmospheres to mixing ratio LBCs, but that the magnitude of the fluxes and the inclusion of negative feedback fluxes do not affect the primary result of pO 2 bistability, we focus our discussion on Cases 1-3. While a full exploration of appropriate H 2 and CO lower boundary fluxes is beyond the scope of this paper, the 'flux-driven modern model' fluxes remain a good choice because of the better prediction of H 2 and CO ground-level mixing ratios, and their ease of comparison with previous 1-D photochemical modelling studies.

Box models
Our results impact on previous work which predicts an abrupt shift from reducing to oxic atmospheres over an extremely narrow oxygen flux range (Alcott et al., 2019;Claire et al., 2006;Goldblatt et al., 2006;Laakso and Schrag, 2017). Despite extrapolating from the fixed-mixing ratio models we have called into question (e.g. Case 0), these box-model studies elucidated that flux-driven feedbacks within the Earth system drive a rapid transition between clusters of trace-O 2 and high-O 2 solutions, while passing through intermediates.
Model atmospheres with these 'intermediate' oxygen compositions have been produced by previous 1-D photochemical modelling efforts, but we have shown that some are unstable equilibrium solutions when explored in flux-driven photochemical models. Minuscule perturbations to lower boundary fluxes seemingly drive any intermediate 'solutions' towards attractors in high-O 2 /O 3 or trace-O 2 /O 3 regions of parameter space, due to positive feedbacks involving the formation of an optically thick ozone column Goldblatt et al., 2006). While issues of kinetics remain to be determined, these atmospheric feedbacks will presumably occur much more rapidly than other biogeochemical feedbacks, which we predict would strongly amplify existing nonlinear behaviour in box models, as well as potentially drive strong feedbacks on microbial ecology. A full exploration of the dynamics of photochemical (in)stability and kinetic timescales for transitions requires time-dependent atmospheric models beyond the scope of this paper. Simultaneous incorporation of these short/mediumterm atmospheric forcings into a long-term biogeochemical model with appropriate Earth system feedbacks would be an interesting target for future work.

Proterozoic pO 2
The lower limit of the high-O 2 /O 3 solutions has interesting applications as a potential constraint on Proterozoic pO 2 . Palaeo-O 2 levels over Earth history compatible with the regions of stability shown in this study are compared to proxy constraints in Fig. 9. Though our model results comment only on the potential stability of palaeo-O 2 levels, as opposed to their temporal evolution, we have truncated our regions of stability to fit with the commonlyaccepted temporal constraints associated with the presence of S-MIF, O-MIF and charcoal. Our results are consistent with the estimates of Proterozoic pO 2 given by the absence of redox-sensitive minerals in Proterozoic sediments (e.g. Johnson et al., 2014), evidence for surface ocean oxygenation , and the presence of variable δ 53 Cr since 1390 Ma (Canfield et al., 2018). They are also consistent with (but stronger than) the lower limit of 2.1 × 10 −4 indicated by the presence of an ozone layer (Segura et al., 2003), as revealed by non-zero O-MIF. However, since Purple regions show the ground-level O 2 mixing ratios of most of our model solutions in the context of existing constraints for pO 2 over Earth history. Grey regions show O 2 mixing ratios compatible with some proxies, from the review by Kump (2008). In addition, existing upper limit constraints are shown with downward arrows and annotations in italics ( * Canfield, 1998* Canfield, , 2005† Farquhar et al., 2000;Pavlov and Kasting, 2002 this constraint arose from models using fixed mixing ratio boundary conditions, our results suggest a revision of the lower limit to 2 × 10 −3 , as our flux-driven models produced very few solutions with both an ozone layer and pO 2 lower than this. Our results are inconsistent with arguments for low mid-Proterozoic pO 2 (e.g. Cole et al., 2016;Planavsky et al., 2014Planavsky et al., , 2020.

Earth system feedbacks and switches between states
Our results also prompt comments on the reversibility of the switch between the trace-O 2 /O 3 and high-O 2 /O 3 states. After an oxic atmosphere is established, is it possible to reduce biospheric fluxes, through a primary productivity crash or ecological reorganisation, to sufficiently low levels (with K oxy < 1) such that the ozone layer collapses and pO 2 falls by several orders of magnitude? The continuous existence of terrestrial eukaryotic life precludes this scenario for most of the Phanerozoic, while the presence of O-MIF (Crockford et al., 2019) after the GOE strongly argues against a return to the trace-O 2 /O 3 state during the bulk of the Proterozoic. Apparently, feedbacks between atmospheric chemistry of ozone formation and biospheric oxygen production are quite strong.
Our model results suggest a previously undescribed stabilising feedback to prevent such reversibility. We have assumed that CH 4 :O 2 flux ratios of 0.5 and 0.1 are representative of reducing biospheres and the modern biogeosystem, respectively. With this in mind, consider a reducing atmosphere supported by a putative low productivity late Archaean biosphere featuring a low O 2 flux and a CH 4 :O 2 flux ratio of 0.49. Following this suite of models (along the arrow labelled (1) in Fig. 2c), an increase in O 2 flux leads to a jump from the trace-O 2 to the oxic state, at an O 2 flux near 8.5 × 10 12 pu in this example. In the resulting oxidising atmospheres, increased methane oxidation in sediments decreases the CH 4 :O 2 ratio towards 0.1 (arrow labelled (2)). Subsequently, the decrease in O 2 flux necessary to plunge back into an anoxic world is far greater than for the higher CH 4 :O 2 flux ratio values (arrow labelled (3)). This constitutes a hysteresis loop preventing small decreases in O 2 flux from reversing atmospheric oxidation. The sulphate-dependent transition to lower CH 4 :O 2 ratios is therefore a positive Earth-system feedback that helps maintain an oxic atmospheric state. In addition to time-dependent photochemical modelling, additional ecosystems modelling work is needed to consider this further, such as a study of the remineralisation dynamics and stability of sedimentary packages during the transition from Archaean-like CH 4 :O 2 flux ratios of 0.5 to modern-like ratios below 0.1.

Conclusions
In this study we have shown that photochemical models utilising fixed flux boundary conditions can contribute towards constraints on atmospheric oxygen over Earth history. Even with the inclusion of long-term negative feedback fluxes of O 2 , atmospheric chemistry drives the atmosphere towards one of two states, between which there are few stable equilibrium solutions. Our results specifically highlight a dearth of stable atmospheric configurations with pO 2 between 6 × 10 −7 and 2 × 10 −3 . Our results support estimates of trace-O 2 levels in the Archaean and O 2 levels greater than 1% PAL during the Proterozoic and Phanerozoic, and highlight the importance of carefully prescribed and described boundary conditions in photochemical models.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.