Ice core evidence for atmospheric oxygen decline since the Mid-Pleistocene transition

The 1.5-million-year-old Antarctic ice indicates a balanced atmospheric oxygen budget before the Mid-Pleistocene Transition.

locked while smaller molecules like Ar, O2, Ne, and He can still move within the firn.
However, even their motion would eventually be inhibited when ice lattice size is further reduced. A critical size of the ice lattice at which close-off fractionation of O2/N2 and Ar/N2 begins is 0.36 nm, concluded by multiple studies (69,71,72). The magnitude of this bubble close-off fractionation appears to be strongly size-dependent, evidenced by the observation that Ne is fractionated 34 ± 2 times than O2 in South Pole firn air (69).
It is worth noting that before the air motions inside the firn were constrained by detailed observations, ice core δO2/N2 was empirically found to have a high coherency with local summertime insolation (24). The subtle influence on ice grain properties by insolation at the time of snow deposition is thought to modulate the magnitude of bubble close-off fractionation (73). Since then, this empirical relationship between δO2/N2 and local summer insolation makes ice core δO2/N2 a useful tool of orbital tuning to constrain ice core chronologies (24,57). δAr/N2 also covaries with insolation and have shared variance with  (22) offers a quantitative estimate of the magnitude of the bubble close-off fractionations. This model shows that microbubbles control the magnitude of δAr/N2 and δO2/N2 fractionation within the firn with a δAr/N2-δO2/N2 slope of 0.55, consistent with the range of the observed slopes reported in Bender et al (13). The covariance between δO2/N2 and δAr/N2 originating from insolation modulation thus emerges as a potential way of deducing the true atmospheric O2/N2 ratio from the ice core δO2/N2 and δAr/N2, where δAr/N2 can be seen as a proxy for insolation.
However, the use of δAr/N2 to correct δO2/N2 for bubble close-off fractionations is complicated by another gas-fractionating process: post-coring gas loss. Gas losses further lower the ice core δO2/N2 and δAr/N2 and enriches δ 18 O of O2. This effect becomes apparent during ice storage at > −50 °C (57,74), when the ice has many fractures in it (23), and if the ice is subject to prolonged pumping [e.g. Allan Hills blue ice δO2/N2 data reported in Higgins et al (53)]. The temperature-dependency is consistent with a physical model that attributes gas loss to the outward diffusion of gases trapped inside the ice (75,76). Intriguingly, gas loss fractionation is still associated with an apparent ~0.5 slope of δAr/N2 against δO2/N2 (23), which is not very different from the slope caused by bubble close-off fractionations.
Two modes of gas loss fractionation-one mass-dependent and the other size-dependentare invoked to explain the observed slope as well as clearly mass-dependent isotopic fractionation (22,23). The slow leakage of gas molecules through the ice lattice is sizedependent and has little fractionation to the isotope ratios, while rapid gas loss associated with core cracking in response to rapid depressurization (and similarly to pumping) is mass dependent (23). The varying degree of each proposed process could cause scatters in the paired δO2/N2-δAr/N2 data but we cannot quantify their relative contribution. However, if gas loss fractionates δO2/N2 and δAr/N2 quantitatively similar to what bubble close-off does, δAr/N2 may still be capable of correcting δO2/N2 in a lump-sum way.
Below, we evaluate the effect of gas losses on O2/N2 and Ar/N2 ratios preserved in the Allan Hills ice (batch 2 and 3 only) by plotting the difference between δO2/N2 replicates (ΔδO2/N2) against the difference between δAr/N2 replicates (ΔδAr/N2) measured on the same depth (Supplementary Data Table 2; Fig. S8). This difference should reflect the varying degree of gas losses experienced by each replicate plus the analytical uncertainties associated with gas handling and mass spec analyses. In the case where more than 2 replicates are available, the pair with the largest δO2/N2 difference is selected. The replicate with higher δO2/N2 does not necessarily experience no gas loss. That sample simply has the least degree of gas losses. As a matter of fact, we have no way of identifying the δO2/N2 or δAr/N2 value that is truly "gas loss-free" for blue ice. At any rate, a slope of 1.47 ± 0.15 (2σ) is found between ΔδO2/N2 and ΔδAr/N2 in the Allan Hills blue ice (Fig. S8), close to the slope observed in the δO2/N2-δAr/N2 plot ( Table 1). The correction by δAr/N2 would therefore also compensate the influence of post-coring gas loss on δO2/N2, in addition to bubble close-off fractionations. When δAr/N2 is corrected to a value smaller than 0 (e.g. the mean value of all observed δAr/N2), the corrected δO2/N2 is expected to have the same degree of gas loss and a constant offset from the true atmospheric δO2/N2 (Fig. 1). When the δAr/N2 is normalized to 0, the effect of gas loss is fully removed. We opted to apply the first type of correction to remove the imprint of insolation in δO2/N2 while retaining a (presumed) constant offset. The offset was subsequently accounted for by the a priori knowledge that the extrapolated δO2/N2 should be 0 at present (Fig. 4) (8). The approach adopted here has smaller uncertainties than extrapolating δAr/N2 to 0 does.
In sum, despite a lack of definitive knowledge on how insolation modulates δO2/N2 and δAr/N2 and how ice physics dictates post-coring gas loss fractionation, the prevalent, consistent, and strong correlation between ice core δO2/N2 and δAr/N2 creates an opportunity to use δAr/N2 for empirical δO2/N2 corrections and infer atmospheric pO2 from the ice core δO2/N2. Unported License). The study area is indicated by the black square. Left, a magnified image of the drilling site from same WorldView03 satellite file. The imagery is processed (gammaadjusted) on ESRI® ArcGIS software package to enhance the color contrast within the blue ice. This color rendition causes blue ice to exhibit a green hue. The brown contours in the ice are exposed dust bands, providing a first-order tracer of surface ice stratigraphy. The locations of the cores reported in this work (see text) are marked with red circles.  19 and 172 m (shown as the dashed rectangle). We use δ 18 Oatm as a stratigraphic marker to better identify the transition. In ALHIC1503, the first appearance of MPT ice is accompanied by a rapid decline of δ 18 Oatm from 1.361‰ to 0.488‰. In ALHIC1502, this transition occurs between 146.40 and 147.99 m (marked by the black arrow). We therefore choose the average depth of 147.20 m as the divide between post-MPT and MPT ice. The deepest 40 Aratm datum marked in red was excluded from age classifications. (25). Six age bins were selected as (A)-(F), named after the average age of all samples in each bin with the full age range shown in brackets. Linear regression results (intercepts and slopes with 95% CI) are also shown. Statistically significant correlations were marked with asterisks after the correlation coefficients. The poor correlation in (D) and (F) likely results from the small range of the observed δAr/N2, leading to large uncertainties in the estimated δO2/N2 at 426 ka and 664 ka, as shown in Fig. 2. However, despite the large errors, the corrected δO2/N2 in (D) and (F) still falls on the trend line of the composite δO2/N2 record within the 95% CI. . This trend is reproduced by corrected δO2/N2 data (orange squares with error bars representing the 95% confidence interval). (B) Paired δO2/N2-δAr/N2 data color-coded according to age. Older samples generally have higher δO2/N2 values compared to younger samples and qualitatively reflect the pO2 decline in the Late-Pleistocene. The red and blue lines are regressions lines of δO2/N2 against δAr/N2 for samples older and younger than 300 ka, respectively, the results of which are shown in the box. The inferred rate of change solely based on these two regression lines is −26.00 ± 12.37 ‰/Myr (95% CI).

Fig. S5. A simplified view on the geochemical cycling of oxygen, represented by electron transfers.
Assuming that O(-II), C(IV), Fe(III), and S(VI) are the most energetically favored redox states, deficit and surplus of electrons represent molecular O2 and reduced carbon, iron, and sulfur species, respectively. Arrows represent the flow of electrons (red: consuming O2; blue: producing O2). For example, the generation of 1 mol of O2 during photosynthesis means the transfer of 4 mol of electrons from O(-II) in water to another electron acceptor. Carbon, iron, and sulfur are common electron acceptors and the internal electron transfer between them (e.g. sulfate reduction) do not impact the oxygen budget. Numbers are estimates from Johnson and Bif (77) for photosynthesis and respiration rates, Galvez (78) for burial and weathering rates, and Walker (79) for reservoir sizes. Note that we do not differentiate igneous sink [e.g. Fe(II) upwelled through the mid-ocean ridges and volcanically outgassed H2S] from the general weather sink, because the igneous sink flux is ~20% of the sedimentary weathering flux and consequently could only affect pO2 on timescales beyond what ice core records can capture (78). These estimates all have their uncertainties, which are not discussed here because evaluating these estimates is beyond the scope of the present study.  Table S1) affects the estimates of slope and intercept. However, the estimated δO2/N2 is not very sensitive to the choice of regression methods when δAr/N2 is normalized to −7.1‰, where two regression lines intersect.   S8. Pair difference for measured δO2/N2 versus δAr/N2, calculated by subtracting the second-measured replicate from the first-measured replicate. For samples more than twice, we select two replicates with the largest δO2/N2 difference. The observed correlation is consistent with fractionation associated with post-coring gas losses in ice cores (23). Notably the slope of 1.474 ± 0.153 (2σ) is close to the observed slopes of δO2/N2 versus δAr/N2 in Allan Hills ice samples (Table 1).