Atmospheric Oxygen Abundance, Marine Nutrient Availability, and Organic Carbon Fluxes to the Seafloor

The global‐scale oxygenation of Earth's surface represents one of the most fundamental chemical transformations in our planet's history. There is empirical and theoretical evidence for at least two distinct and stable regimes of Earth surface oxygenation—a “low‐O2 world” characterized by pervasively anoxic deep ocean waters, and a “high‐O2 world” with dominantly well‐oxygenated deep ocean waters represented by our modern surface environment. Numerous biogeochemical processes and feedbacks control the redox state of the marine system, particularly when considered globally and on geologic timescales. It has therefore proven challenging to provide quantitative and internally consistent estimates of the atmospheric oxygen levels (and thereby, productivity, nutrient availability, and reductant consumption) necessary to oxygenate the deep seas. Here, we leverage an Earth‐system biogeochemical model that tracks the carbon, nitrogen, oxygen, phosphorus, and sulfur cycles to provide new quantitative constraints on this relationship. We explore ocean biogeochemistry and fluxes of reduced carbon to the seafloor across a wide range of atmospheric oxygen levels from 0.01% to 100% of the present atmospheric level (PAL), and implement a stochastic approach to provide formal estimates of uncertainty on our results. We find that deep ocean waters remain largely anoxic, and ocean productivity remains significantly muted relative to the modern marine biosphere, until pO2 levels reach ∼40% PAL. These results have major implications for quantitative constraints on atmospheric pO2 levels during the latest Proterozoic and Paleozoic, both in terms of environmental habitability for early animals and with respect to potential energetic constraints on growing and diversifying benthic communities.

of the mid-Proterozoic likely did not exceed ∼10% present atmospheric level (PAL), and may have been well below this for much or most of mid-Proterozoic time (e.g., Bellefroid et al., 2018;Cole et al., 2016;Planavsky et al., 2014Planavsky et al., , 2018Planavsky et al., , 2020S. Zhang et al., 2016;X. Liu et al., 2016). Moving forward in geologic time, constraints on atmospheric composition have been tied to the onset of the charcoal record in the Late Silurian (∼420 Ma; Glasspool et al., 2004), with more recent work indicating that atmospheric oxygen levels of ∼75% PAL may have been required for self-sustaining fire propagation (Belcher & McElwain, 2008;Belcher et al., 2010). As a result, the transition of the Earth system from a relatively stable mid-Proterozoic world to something closer to the modern state, perhaps by the late Paleozoic (e.g., Dahl et al., 2010;Lenton et al., 2016;Wallace et al., 2017), represents a vast period of time across which we have a relatively poor understanding of atmospheric oxygen levels.
Our primary source of information about the evolving state of Earth's surface redox conditions comes from a range of proxies that reflect the biogeochemistry of marine environments. This includes geochemical redox proxies (e.g., F. Zhang et al., 2018;Lu et al., 2018;Partin et al., 2013;Reinhard et al., 2013;Sperling et al., 2015;Stockey et al., 2020;Wallace et al., 2017) as well as paleontological evidence for evolving ecosystem structure and environmental habitability (e.g., A. G. Liu et al., 2015;Bowyer et al., 2017;Tarhan, 2018). These records provide substantial information about marine biogeochemical conditions; however, these signals are difficult to quantitatively link to atmospheric oxygen abundance. For example, there is currently no consensus on the concentration of atmospheric oxygen that would be required to observe pervasively "oxic" signatures in the trace metal records of marine systems, nor are there firm quantitative links between the environmental oxygen levels implied by a given proxy or set of proxies and the levels required for a particular degree of organism or ecosystem complexity .
There has been recent work exploring the relationship between primary productivity and proxy-based estimates of mid-Proterozoic oxygen levels Laakso & Schrag, 2019;Ozaki et al., 2019). However, these analyses have been focused on the fairly limited range of atmospheric oxygen levels expected for Earth's middle history, without examining the transition to higher pO 2 . Employing large-scale biogeochemical models across more than ∼2-3 orders of magnitude of atmospheric oxygen levels is challenging, since highly non-linear feedbacks are expected as the ocean becomes oxygenated. Modeling these feedbacks requires explicit coupling of the Earth's carbon, oxygen, nitrogen, phosphorus, and sulfur cycles (CANOPS), all of which play key roles either as nutrients for the biosphere or as redox regulators, resulting in significant computational expense. In addition, explicit quantification of key model uncertainties often requires a large ensemble of model runs. Together, these constraints have prevented explicit investigation of ocean ventilation on productivity across a large oxygen range except in very simple model architectures and forward models (e.g., Alcott et al., 2019;Laakso & Schrag, 2017).
Beyond oxygenation of the water column, the transition to a well-oxygenated deep ocean has major implications for nutrient cycling and dramatic shifts in nutrient availability across the food chain-and is thus likely to impact the viability of evolving benthic ecosystems. As a result, it is important to establish a framework for the possible range of potential Earth system states that would be consistent with observations from the rock record. By applying tested biogeochemical modeling methods to this question using increased computing capacity and a moderate complexity model (Ozaki & Tajika, 2013;Ozaki et al., 2011), we provide a framework to quantitatively tie atmospheric conditions to nutrient availability, global marine redox conditions, and biospheric productivity across the pO 2 space relevant for the last ∼2.3 billion years of Earth's history.

Methods
Here, we utilize a biogeochemical Earth system model (CANOPS) in order to identify biogeochemical regimes stable on geologic time scales at atmospheric oxygen levels between 0.01% and 100% PAL. The CANOPS model couples a biogeochemical model with a diffusion-advection model of the global ocean, a parameterized sediment model, and a stagnant film model for air-sea gas exchange. The ocean circulation model robustly reproduces modern profiles of ocean circulation tracers, and the biogeochemical and marine sediment diagenesis models include explicit representations of photosynthetic primary production, a complete series of heterotrophic respiratory pathways, a series of primary and secondary redox reactions, and the deposition, decomposition, and burial of biogenic material in marine sediments (Ozaki & Tajika, 2013;Ozaki et al., 2019). The physical configuration of the CANOPS model is shown in Figure 1a, while our modern "benchmark" simulation is compared to observations from the modern oceans in Figures 1b and 1c. A detailed description of the model, modern Earth system calibration, and all parameterizations and variables can be found in Ozaki and Tajika (2013) and Ozaki et al. (2019).
We build on the previously employed structure of CANOPS in two ways. First, we extend the model's treatment of pyrite (FeS 2 ) formation to include the formation of pyrite in the water column in addition to formation in marine sediments. We consider this a particularly relevant addition to the model framework given our interest in a very wide range of oxygen and sulphate levels and in the biogeochemical transition between largely anoxic and strongly oxygenated systems, which is likely to be accompanied by an attendant large shift in the dynamics of global sulfur cycling. We parameterize water column FeS 2 formation at a given depth as a function of the local availability of dissolved sulfide and an assumed concentration of dissolved Fe 2+ (Dale et al., 2009(Dale et al., , 2015: where WC py represents the water column pyrite formation flux, k py represents a rate constant for pyrite formation, and brackets denote concentration. Second, we employ a revised parameterization for the efficiency of nutrient P scavenging (σ scav ) by Fe-bearing minerals (Fe-P trap) as a function of the redox state of the ocean interior. Specifically, and in contrast to Ozaki et al. (2019), we specify that the relative efficiency of nutrient P scavenging is dependent on the abundance of dissolved O 2 below the photic zone according to: where max scav denotes a maximum scavenging efficiency, sampled assuming a log uniform distribution from a range between 0.01 and 1 during the stochastic analysis, [O 2 ] j = 1 denotes the dissolved oxygen abundance in the ocean layer below the photic zone, and [O 2 ] 0 represents a reference dissolved O 2 concentration that initiates redox dependence of nutrient P scavenging. This is set by default to [O 2 ] = 1 μmol kg −1 following Reinhard et al. (2017) as a threshold centering value for this continuous function between effective scavenging at low [O 2 ] and no scav- enging at higher [O 2 ]. Mechanistically, this parameterization is meant to describe the scavenging and removal of nutrient P below the photic zone when the ocean interior becomes pervasively anoxic. The P scavenging flux is thus equal to the upwelling flux of P to the photic zone multiplied by σ scav . This removal could be due to scavenging and coprecipitation by Fe-oxide mineral phases at the oxic-anoxic interface (Bjerrum & Canfield, 2002;C. Jones et al., 2015), or the removal of nutrient P as a constituent of other reduced Fe-bearing minerals such as Fe-phosphates, green rust, or Fe-silicate phases (Derry, 2015;Zegeye et al., 2012). It is critical to point out that, the specifics of Fe-associated P removal when the oceans are pervasively oxygen-poor are not well constrained, especially in three key respects; (a) the sensitivity of P removal efficiency to water column oxygen (although see Turnewitsch and Pohl (2010) for exploration of this in a modern setting); (b) the impact of changes in ocean chemistry (i.e., carbonate chemistry, pH, etc.) on the efficiency and viability of each of these mechanisms (e.g., Jahnke, 1984;Reinhard et al., 2017); and (c) the functional form of the parameterization. This highlights the need for a robust statistical approach and examination of the sensitivity of this parameter in our model (Figures 3a  and 3b).
An additional key difference between our analysis and previous work (Ozaki et al., 2019;Reinhard et al., 2017) is that we make the simplifying assumption that the C:P ratio of photosynthetic biomass is constant at the classical "Redfield" value of 106:1. There is evidence to suggest that biomass C:P should respond dynamically to environmental conditions (e.g., Galbraith & Martiny, 2015;Reinhard et al., 2017), such that the C:P ratio of primary producers increases as nutrient P availability drops (Quigg et al., 2003), although the behavior of this response is poorly constrained. As such, for simplicity and transparency we retain the simple "Redfield" assumption in our analysis, but conduct a sensitivity test using a higher set C:P ( Figure 3c). We find a substantial difference in productivity as would be expected, and this sensitivity test represents the most dramatic variation in results. This suggests that understanding the quantitative relationship between the globally integrated C:P ratio of the ocean biosphere and environmental boundary conditions is an important topic for future work, as is the attempt to provide empirical constraints on biomass C:P in ancient oceans.
The centerpiece of our overall approach is to run a very large ensemble of model simulations in which multiple key parameters are simultaneously sampled randomly from an assumed prior distribution, with the model subsequently run with atmospheric pO 2 enforced as a constant boundary condition until the S and P cycles reach steady state, although because S residence time is much longer, model run times are determined primarily by S balance. We focus here on six key parameters, chosen for their potential to significantly impact global oxygen and nutrient cycling, and their likely roles in controlling the mechanistic links between atmospheric pO 2 , nutrient biogeochemistry, marine redox, and energy fluxes to marine sediments. In contrast with Ozaki et al. (2019), we do not vary the size of the crustal sulfur reservoirs for simplicity and because reservoir size was not found to have a strong effect on the results. The key control parameters in our analysis along with their constrained ranges and assumed prior distributions are provided in Table 1. The specified ranges and prior distributions are discussed briefly below. Large ensembles were implemented within the Georgia Institute of Technology Partnership for an Advanced Computing Environment (PACE), with downstream data analysis performed using a custom pipeline built in Python. All model code, output data, and analytical pipeline can be found at the https://doi.org/10.5281/ zenodo.5517932.
The key boundary condition for our analysis is atmospheric pO 2 . We examine atmospheric oxygen levels ranging between 0.01% and 100% PAL, assuming a log-uniform prior distribution. By exploring atmospheric oxygen levels across five orders of magnitude, we are able to provide both a repeated analysis of low pO 2 levels and confirm earlier findings of Ozaki et al. (2019) with our newly updated version of CANOPS, as well as place these results into larger context relative to the modern Earth system with stable biogeochemical solutions up to 100% PAL O 2 . Most significantly, this pO 2 range is significantly expanded beyond that of Ozaki et al. (2019) on the high-O 2 end, and is expected to cover the dynamics of the transition to the well-ventilated and highly productive "modern" ocean state.
The half-saturation constant for microbial sulfate reduction (K MSR ) controls the rate at which organic matter is broken down via MSR at a given concentration of SO4 2− . Because estimates of K MSR in natural environments and pure cultures vary over several orders of magnitude (Pallud & Van Cappellen, 2006;Tarpgaard et al., 2011), we implement a relatively wide range of K MSR values of 0.002-2.0 mM, again assuming a log-uniform prior distribution. Based on previous work (Ozaki et al., 2019), it is expected that our results are not particularly sensitive to any plausible expansion of this range. However, given that in this study we are examining a much larger range and higher concentrations of sulfate we have chosen to maintain this parameter as a component of our stochastic analysis.
The vast majority of P is delivered to the oceans via rivers, making the riverine P flux (R P ) a key control on the size of the marine P reservoir. To account for relatively unconstrained variability in this flux-tied to changes in the composition of weathering crust or the colonization of the continents by land plants, for example-we have varied this flux from 0.2 to 2.0 times the modern value, assuming a uniform normal prior distribution. Our range is similar to that of Ozaki et al. (2019), although we limit our minimum value to 0.2 rather than 0 under the pretense that there should always be some non-trivial riverine reactive P flux once there are significant continental land masses above sea level. While Laakso and Schrag (2014) suggest that this flux could be up to 10 times smaller than the modern value this lower limit remains poorly constrained and we chose to limit our range to a single order of magnitude. We note that this parameter is included separately from the erosion/sedimentation rate (see below) as composition and chemical weathering environments will play an independent (albeit not totally decoupled) role in the release of P from the continents.
Given the centrality of the marine P reservoir size to biospheric productivity and the chemistry of the ocean interior, we also include a parameter designed to explore uncertainty associated with the primary sink of P-burial in marine sediments. This parameter, which can be thought of mechanistically as a nutrient P scavenging efficiency ( scav max ), is a proportionality coefficient between 0.01 and 1 that dictates the efficiency of P scavenging (removal from the water column) by oxidized and reduced iron species in anoxic settings (Bjerrum & Canfield, 2002;Derry, 2015;Laakso & Schrag, 2019;Ozaki et al., 2019;Reinhard et al., 2017). We explore this range assuming a log-uniform prior distribution. It is important to note that, as shown in Equation 2, this efficiency coefficient is implemented within a parameterization that describes the dependence of P scavenging on oxygen availability within the water column. That is, regardless of the randomly selected efficiency coefficient, if the water column is well-oxygenated P scavenging will not be operative. The sensitivity of our results to this parameterization is explored below and in Figures 3a and 3b.
The sinking velocity of marine organic matter (V POM ) is a critical parameter regulating the distribution of O 2 demand and nutrient release through the oceanic water column (e.g., De La Rocha & Passow, 2007;Devol & Hartnett, 2001;Kwon et al., 2009;Meyer et al., 2016). It has also been suggested that this parameter would have changed dramatically over time in step with major changes to the ecological structure of the surface marine biosphere (e.g., Butterfield, 1997;Lenton et al., 2014), although recent work casts doubt on major changes to particle sinking velocities through time (Fakhraee et al., 2020). Nevertheless, we assign the global sinking velocity of particulate organic matter to be between 10 and 100 m d −1 , assuming a uniform normal prior distribution, to quantify the large uncertainty associated with settling rates of organic matter to the deep ocean floor, as well as the idea that sinking velocities during the Proterozoic may have been slower due to the lack of larger eukaryotic cells and limited packaging in zooplankton fecal pellets. Our range is similar to that of Ozaki et al. (2019), although we limit our minimum value to 10 m d −1 rather than 0 m d −1 , as the latter is unlikely to be a more physically reasonable endmember.
We implement a single scaling parameter for continental erosion and marine sedimentation rates (f sr ), as these should be closely linked on a global scale when the rock cycle is at steady state. This parameter is normalized to that of the modern Earth (e.g., sr 0 = 1.0). This is a key variable in our model since these rates modify the rate of oxidative weathering of organic carbon and pyrite, as well as the rate of burial of these species in marine sediments, though we note that these fluxes are also controlled by additional variables within the parameterization including oxygen availability and species concentration. The sediment accumulation rate at the seafloor also affects the burial efficiency of P and organic matter and the O 2 penetration depth in the sediment column. Here, we vary this parameter from 0.5 to 1.5, assuming a uniform normal prior distribution, based on the premise that it is equally reasonable for globally integrated rates of erosion and sediment burial to have been either lower than the modern-particularly during the Precambrian (e.g., Husson & Peters, 2017)-or higher, such as during the Pliocene-Pleistocene (e.g., Herman et al., 2013).
After the generation of our complete data set (n = 20,589 models), the data were subsampled to produce a suite of results that are also consistent with our basic understanding of the Earth system as derived from the geochemical record. This includes constraints on marine SO 4 concentrations and an upper limit on globally integrated rates of N fixation. Specifically, at atmospheric oxygen levels between 0.01% and 10% PAL (low pO 2 group) we subsampled the model ensemble such that 0.05 mM < [SO 4 ] ≤ 10.0 mM, while at atmospheric oxygen levels above 10% PAL (high pO 2 group), we subsampled such that 1.0 mM < [SO 4 ] ≤ 60.0 mM. This upper limit reflects the maximum sulfate concentration the model achieved (i.e., no runs were excluded by this limit), however, this limit should reflect the existence of an "evaporite ceiling" (e.g., Canfield & Farquhar, 2009) which is poorly constrained in natural systems, while the lower limit is constrained by the sulfur stable isotope record of Proterozoic sedimentary rocks (Luo et al., 2015;Lyons & Gill, 2010;Planavsky et al., 2012;Scott et al., 2014). We chose the transition for these two sulfate bins to be 10% PAL because a natural gap in the data occurs at this atmospheric oxygen level-that is, prior to any filtering, our model does not find [SO4] > 1 mM to be feasible until atmospheric oxygen is above 10% PAL ( Figure S1 in Supporting Information S1). These conservative and overlapping ranges were chosen specifically to remove mathematically viable solutions for which there is no evidence in the rock record. In particular, this scheme was designed to remove cases of high atmospheric oxygen and extremely high productivity resulting in extremely low SO 4 and eutrophic oceans. We emphasize that this SO 4 constraint was imposed in an effort to remove unrealistic end-member scenarios, and the overlapping ranges were chosen so as not to impose a false apparent bistability in SO 4 ( Figure S1 in Supporting Information S1).
We also subsample our overall model ensemble for cases in which globally integrated rates of N fixation are less than 10 times the modern level. Although we consider this cutoff reasonable, and note that our primary results are not particularly sensitive to this assumption (for instance, increasing the cutoff to 15 times modern only increases the size of our final SO 4 filtered data set by 0.4%), it is important to point out that this upper limit is almost completely unconstrained for the Earth system. While P is widely considered the ultimate limiting nutrient on geologic timescales (Tyrrell, 1999), the mathematical potential for extremely high N fixation rates indicates that there may be instances where global primary productivity could possibly be N-limited at high atmospheric oxygen levels. This would likely occur via more proximate limitation on the bioavailability of iron which is required as a catalyst in nitrogenase (Falkowski, 1997;Raven, 1988). There is also some potential that in a high oxygen, high productivity world trending towards eutrophic oceans, trace element (Mo and V) limitation of N fixation as a result of expanding anoxic environments could act as a negative feedback on runaway eutrophication. The potential for N-limited productivity tied to availability of cofactors (Fe, Mo, and in some cases V) for nitrogenase or the loss of fixed N has been explored previously (Falkowski, 1997;Fennel et al., 2006;Reinhard et al., 2013) but not in the context of extremely high rates of productivity well beyond the modern. This represents a promising avenue for future work.

Results
Subsampling of our overall model ensemble yields the primary data set used in our analysis (n = 1,672 models). Although the subsampled ensemble is still very large, subsampling results in a significant drop in the overall model ensemble size. Nevertheless, median values of global diagnostics in our subsampled ensemble suggest a stochastic ensemble of sufficient size to attain a stationary distribution (Figure 2b). We thus expect our subsampled results to be statistically representative. In some cases we have further divided the "high pO 2 " group based on observations of the inflection points in model behavior as we consider it useful to distinguish "moderate pO 2 " (10%-40% PAL) from "high pO 2" (40%-100% PAL).
There is negligible dissolved O 2 within the water column at pO 2 below ∼10% PAL ( Figure 5). At abyssal depths we begin to observe appreciable levels of O 2 between 10% and 20% PAL, with the most dramatic increase above 40% PAL. Nevertheless, our model results do not achieve a median [O 2 ] that would be considered fully oxic at abyssal depths until atmospheric oxygen levels are above 60% PAL. Correspondingly, below 1% PAL, [PO4 3− ] is extremely low with a median of ∼0.1 μmol, roughly less than 1% of modern values ( Figure 5). Median deep water phosphate concentrations remain below 1 μmol up to 30% PAL and only approach roughly modern values when atmospheric oxygen levels are above ∼60% PAL.
The depositional flux of organic carbon at pO 2 levels below 10% PAL is about an order of magnitude lower at all water depths than the benthic C flux when pO 2 is above 10% PAL (Figure 6). At 200 m depth, median benthic C org flux is 0.062 ± 0.088 (1σ) mol C m −2 yr −1 , and at 4,000 m depth median C org flux is 0.005 ± 0.006 (1σ) mol C m −2 yr −1 . Above 10% PAL, the median flux increases substantially; when pO 2 is 10%-40% PAL, median benthic C org flux is 0.497 ± 0.137 (1σ) mol C m −2 yr −1 at 200 m depth and 0.032 ± 0.012 (1σ) mol C m −2 yr −1 at 4,000 m depth. These values are within roughly 1σ of medians when pO 2 is above 70% PAL. At these highest pO 2 levels, median benthic C org flux is 0.944 ± 0.466 (1σ) mol C m −2 yr −1 at 200 m depth and 0.066 ± 0.031 (1σ) mol C m −2 yr −1 at 4,000 m depth. We note that these are zonally averaged values and benthic C org flux is very heterogeneous in the modern ocean, however, these values are within range of modern data (e.g., Sweetman et al., 2017).

Ventilation of the Deep Ocean
Our model results indicate that until atmospheric oxygen levels increase above ∼40% PAL, the deep ocean remains largely oxygen-poor, nutrient depleted, and significantly less productive relative to the modern ocean ( Figures 4 and 5). This value is in line with the results of previous studies (e.g., Canfield, 1998), however, the approach-which required modern levels of P in the deep sea-requires an Earth system state for which there is no longer strong support (e.g., Cole et al., 2020). It is important to note that the impacts of oxygen limitation on benthic ecosystems and the expression of this limitation in the geochemical proxy record differ substantially. Specifically, as atmospheric oxygen levels increase and the marine system begins to respond more dramatically, most geochemical proxies will tend to act as an "on-off switch" at the onset of this transition. In contrast, marine fauna-particularly larger, more complex organisms (e.g., fishes) and ecosystems-will likely feel the effects of lower-than-modern atmospheric oxygen concentrations up to a pO 2 of ∼60% PAL at which point globally averaged deep waters will achieve [O 2 ] similar to the modern.
Below ∼10% PAL, deep water [O 2 ] concentrations are essentially negligible, or would on average be low enough to suppress aerobic respiration and result in geochemical proxy signatures diagnostic of reducing environmentsthat is, environments characterized by denitrification, reduced iron, or sulfide and traditionally labeled as anoxic to suboxic ( Figure 5). At these low atmospheric oxygen levels, we also find that water column phosphate is strongly limited, with [PO 4 3− ] concentrations in deep waters similar to or lower than surface waters in the modern ocean ( Figure 5a). We note that the basin-averaged water column chemistry of CANOPS is not equipped to model the possibility of localized and/or temporally variable weakly oxygenated conditions in the ocean interior at relatively low atmospheric pO 2 . Nevertheless, these findings are consistent with the consensus of typical background mid-Proterozoic conditions, at least globally averaged and on long time scales (e.g., Laakso & Schrag, 2019;Ozaki et al., 2019;Partin et al., 2013;Poulton & Canfield, 2011;Reinhard et al., 2013Reinhard et al., , 2017Scott et al., 2008).
As atmospheric oxygen levels increase, we observe the onset of a more dramatic response in the marine biogeochemical system. From 10% to ∼40% PAL, oxygen availability in globally averaged deep waters remains quite limited, with oxygen concentrations that would conventionally be classified as dysoxic Tyson & Pearson, 1991). At these concentrations, there is enough oxygen present in the water column to remove ferrous iron and inhibit anaerobic respiration, but negative impacts on benthic ecology would still be potentially significant (e.g., Boag et al., 2018). While we see the most dramatic increase in [O 2 ] above ∼40% PAL, oxygen availability is likely to have impacted benthic organisms even up to ∼60% PAL ( Figure 5). Within this range of atmospheric oxygen levels, localized conditions such as increased temperature or salinity would not only decrease O 2 solubility, but would also push metabolic rates higher thereby compounding the impacts of oxygen availability (Boag et al., 2018;Cole et al., 2020;Pörtner, 2012;   Reinhard et al., 2016), though the inverse is also true in that localized oxygen oases are equally possible. Similarly, we find [PO4 3− ] remains strongly limited up to 40% PAL, and somewhat less so up to 60% PAL as the strength of the Fe-P trap decreases with increasing water column [O 2 ] (Figure 5a). This indicates that the potential of significant environmental impacts of ocean redox on the habitability of benthic environments for larger, more complex organisms should be expected below ∼60% PAL. We highlight that this is for globally averaged values, Figure 4. (a) Export production and (b) marine phosphate reservoir size as a function of atmospheric oxygen levels. Gray bar denotes range of modern estimates for export production (Dunne et al., 2007;Heinze et al., 2009;Laws et al., 2000;Sarmiento and Gruber, 2006) and the marine phosphate reservoir (Delaney, 1998;Guidry et al., 2000). Red indicates the low O 2 group (pO 2 < 10% present ocean level [PAL]) while blue indicates the high O 2 group (pO 2 > 10% PAL).  Tyson and Pearson (1991), and further elaborated on by Cole et al. (2020). Anoxic/Suboxic refers to vanishingly low [O 2 ], the presence of anaerobic metabolisms, heterotrophic aerobic bacteria. Dysoxic refers to a lack of anaerobic respiration but oxygen levels low enough to impact benthic ecology. Oxic is similar to modern surface waters.
suggesting that these effects would impact the majority of the marine environment. In contrast, at atmospheric oxygen levels above 60% PAL, both [O 2 ] and [PO4 3− ] are consistent with the range of concentrations in deep waters of the modern oceans, with enough phosphate available to support roughly modern levels of primary productivity ( Figure 5). However, non-trivial albeit increasingly localized portions of the marine environment will still likely have oxygen levels low enough to substantially impact the extent of habitable space for larger, more complex organisms. This is true even up to 100% PAL, as is observed with recent deoxygenation and impacts on the world's fisheries (e.g., Pauly and Cheung, 2018;Stortini et al., 2017), although three-dimensional ocean biogeochemistry models with higher spatial resolution are required to explore these effects.

Nutrient Availability and Productivity Dynamics
Our model results suggest global export production is strongly muted at atmospheric oxygen levels below 10% PAL relative to the modern (median value ∼7% of modern estimates) in agreement with previous work (Laakso & Schrag, 2019;Ozaki et al., 2019; Figure 4, Table 2). As atmospheric oxygen levels increase, median productivity remains depressed below ∼50% of modern estimates until pO 2 > 40% PAL. These limited levels of primary productivity can be tied directly to a lack of available nutrients (specifically PO4 3− ) in the marine system (Figure 4). The flux of riverine P delivered to the ocean is included in our stochastic analysis, however, this is not the primary control on the marine P reservoir size across four orders of magnitude in atmospheric pO 2 . Instead, marine PO4 3− availability is dictated by the efficacy of removal from the marine system, which in our model is parameterized based on the hypothesis that pervasively anoxic and Fe-rich oceans will give rise to a deep sea P   Table 2 Global Data Summary trap (e.g., Bjerrum & Canfield, 2002;Derry, 2015;Laakso & Schrag, 2014;Reinhard et al., 2017). Our results further highlight the development of a better understanding of the mechanistic links between pO 2 , ocean Fe inventory, and marine P availability as an important topic of future work.
While it is widely thought that P was the ultimate limiting nutrient for primary productivity globally and on geologic timescales through the mid-Proterozoic (Derry, 2015;Laakso & Schrag, 2014, 2018Ozaki et al., 2019;Reinhard et al., 2017), there are instances where N or other trace nutrients such as Fe, Mo, or V may play a key role in capping marine productivity. In our model results, we examine the response of the N cycle to increasing pO 2 and find that below 10% PAL median N fixation is <10% of modern estimates (median of 0.92 ± 3.5 (1σ) Tmol N yr −1 ), indicating that trace nutrient limitation is unlikely under these conditions, given that denitrification would be limited in a largely anoxic ocean ( Figure S2 in Supporting Information S1). Above 10% PAL, N fixation rates increase most dramatically, reaching within error of modern estimates; from 10% to 40% PAL, median N fixation is 9.24 ± 6.66 (1σ) Tmol N yr −1 , or ∼95% of modern estimates ( Figure S2 in Supporting Information S1). As atmospheric oxygen levels increase further, so too do median rates of N fixation. However, for a given level of export productivity, we observe that N fixation rates are highest at moderate levels of pO 2 (∼10-40% PAL) providing more robust support for the relationship observed by Reinhard et al. (2017). Reinhard et al. (2017) suggested that this may be the result of relatively higher productivity due to a weaker P scavenging sink, combined with still-limited water column O 2 (relative to the modern), leading to bioavailable N loss via denitrification and anaerobic ammonium oxidation. Our results support this hypothesis, with increased organic C burial and export productivity at pO 2 > 10% while water column [O 2 ] remains suppressed.
While the oxygenation of the marine system would lead to a significantly increased supply of dissolved and bioavailable trace nutrients critical for nitrogenase (Mo and V), this would lead to a decrease in dissolved and bioavailable Fe-also a necessary cofactor for nitrogenase (Falkowski, 1997;Raven, 1988). It is important to stress that within the current CANOPS model it is assumed N fixation will always keep pace with productivity, which is defined by P availability. However, further development of explicit Fe and trace metal cycling represents a promising avenue for future work. In any case, we find that some scenarios at higher atmospheric oxygen levels can result in extremely high rates of N fixation (see Section 2). In these cases, Fe availability could limit the capacity of N fixation which in turn may have played a key role in preventing the Earth system from moving into an extremely high productivity/low O 2 (eutrophic) regime. Under such conditions there is potential for a further stabilizing negative feedback whereby rapidly increasing productivity and N fixation could lead to higher demand for trace nutrients alongside water column deoxygenation resulting in an increased drawdown of necessary trace nutrients Mo and V, though this would likely be more localized. A nitrogen fixation cap on productivity tied to trace metal limitation has been explored in the context of Archean conditions (Zerkle et al., 2006). However, we only have a qualitative sense of this type of productivity cap in a high pO 2 /high productivity regime, highlighting a promising avenue for future work. This work would be aided by implementation of a complete iron cycle in this or similar biogeochemical models (e.g., van de Velde et al., 2020), as well as incorporation of trace metal cycling, and would have potentially broad applications for both Earth system evolution and predictions regarding Earthlike exoplanetary productivity regimes.

An Evolving Biosphere
While oxygen availability is traditionally the focus for thinking about habitability of benthic environments and the evolution and diversification of animals, considering the availability of food is also critically important Reinhard et al., 2020;Sperling and Stockey, 2018). At low pO 2 (<10% PAL), we find that the PO4 3− reservoir is limited to ∼7% of modern estimates. This extremely nutrient-limited regime-possibly characteristic of much of the mid-Proterozoic-is expected to strongly favor the dominance of small phytoplankton cells, limiting the expansion of eukaryotic ecosystems . In contrast, increased nutrient availability inherently tied to ocean ventilation may have acted as an important "bottom up" driver, providing a strong control on dominant cell size (e.g., cyanobacteria vs. photosynthetic eukaryotes) as well as the extent of the size spectrum, and expansion to higher trophic levels (Armstrong, 1999;Reinhard et al., 2020;Ward et al., 2014). Reinhard et al. (2020) suggested that this transition would require marine nutrient inventories above ∼10% of modern [PO 4 3− ], and in line with this prediction, we observe roughly an order of magnitude increase in median benthic carbon flux between the 0.01%-10% and 10%-40% PAL binned results ( Figure 6). As such, it is likely that algal expansion, increasing micro predation, and increasing heterotrophic complexity-mechanistically linked to the expansion of the nutrient reservoir-can be tied to atmospheric oxygen levels of at least ∼10% PAL ( Figure 5). Significantly, at pO 2 of 10% PAL, water column oxygen would have still been vanishingly low in the global average, and yet possibly still higher than estimated requirements of the earliest metazoans (Mills et al., 2014;Sperling et al., 2013). The resolvability of these nutrient-redox dynamics in our results provides further strong support for the notion that ventilation-associated changes to nutrient availability may have been more important for the evolutionary progress of early animals than simply just the oxygenation of the water column.
This expansion of the nutrient reservoir and subsequently of trophic structure and heterotrophic complexity would have been a necessary precursor to the arrival of more complex animals with relatively higher oxygen demands towards the late Neoproterozoic and Cambrian. This is especially relevant as an improved understanding of oxygen requirements for early animals and their communities is emerging (e.g., Boag et al., 2018;Levin, 2003;Sperling et al., 2016). As energy fluxes and water column [O 2 ] continue to increase as pO 2 rises, this dramatic shift may act as a positive feedback to support increasingly large and diverse fauna well into the Paleozoic.
Even in the modern ocean, the majority of the marine environment and especially the abyss is characterized by substantial food limitation (<0.167 mol C m −2 y −1 ; Lutz et al., 2007;Watling et al., 2013). Despite major variation in benthic C org flux across depth, latitude, seasonality, and sea-ice regimes, this is roughly double the median flux retrieved by our model ensemble for 200 m depth when pO 2 < 10% PAL. Broadly, this suggests that food limitation would have been a major factor in controlling benthic biodiversity throughout the ocean and strongly favored the microbial community (Sweetman et al., 2017). Recent stu7dies have explored the effects of declining POC flux as a result of climate change and suggest that a 3-fold reduction in POC may reduce nematode and benthic microbial biomass by ∼50%, macrofaunal biomass by 80%, and lead to major decreases in bioturbation, benthic respiration, and sediment mixed-layer depth (D. O. Jones et al., 2014;Laws et al., 2000;Levin et al., 2009;Smith et al., 2008). While our results lack spatial resolvability, it is clear that with an increase to higher pO 2 , especially above ∼50% PAL, benthic food supply would increase substantially and the area of severely food-limited regions of the ocean would shrink, potentially leading to increased ecological space for non-microbial life.

Conclusions and Key Findings
Across four orders of magnitude of pO 2 , ranging from levels representative of the Proterozoic Eon to the modern Earth, our analysis identifies viable biogeochemical regimes that are potentially capable of sustaining prescribed atmospheric conditions. We find that the most dramatic transformation of the marine biogeochemical system occurs between 10% and 40% PAL, with the water column not reaching modern-like oxygen and nutrient levels until atmospheric oxygen exceeds 60% PAL. The onset of ocean ventilation also marks a critical increase in energy flux to the benthic environment as export productivity increases in response to a muted Fe-P trap. Heightened food supply across this transition likely helped to drive the expansion of increasingly complex life across the Neoproterozoic and well into the Paleozoic. We stress that this work does not explain the evolution of atmospheric oxygen through time since we dictate pO 2 as a boundary condition. However, future efforts should aim to allow the secular evolution of pO 2 in order to further refine our solution space. In sum, our results provide a new quantitative framework for linking insights about marine redox-productivity dynamics gleaned from the rock record with atmospheric composition across an interval that marks one of Earth's most transformative biogeochemical transitions.

Data Availability Statement
Data and code archiving is in compliance with the FAIR data guidelines. All relevant data generated for this study and the code used to generate and analyze the data can be found at https://doi.org/10.5281/zenodo.4716158.