Destabilization of deep oxidized mantle drove the Great Oxidation Event.

The rise of Earth's atmospheric O2 levels at ~2.4 Ga was driven by a shift between increasing sources and declining sinks of oxygen. Here, we compile recent evidence that the mantle shows a significant increase in oxidation state leading to the Great Oxidation Event (GOE), linked to sluggish upward mixing of a deep primordial oxidized layer. We simulate this scenario by implementing a new rheological model for this oxidized, bridgmanite-enriched viscous material and demonstrate slow mantle mixing in simulations of early Earth's mantle. The eventual homogenization of this layer may take ~2 Ga, in line with the timing of the observed mantle redox shift, and would result in the increase in upper mantle oxidation of >1 log(fO2) unit. Such a shift would alter the redox state of volcanic degassing products to more oxidized species, removing a major sink of atmospheric O2 and allowing oxygen levels to rise at ~2.4 Ga.

INTRODUCTION

The rise of atmospheric O2 was one of the greatest transitions in Earth’s history, changing the surface mineralogy, overturning the biota to more oxygen-resistant species, and facilitating the diversification of eukaryotes (, ). However, what ultimately caused the largest increase in O2 2.4 billion years (Ga) ago—known as the Great Oxidation Event (GOE)—remains contentious. A plethora of mechanisms have been suggested (–). Particularly puzzling is the lag between the appearance, based on geochemical proxies, of low levels (“whiffs”) of environmental O2 during at least two Archaean Oxidation Events and of biological O2 producers at 3.5 Ga or earlier, and the GOE at ca. 2.4 Ga (). During this time, geological evidence suggests a generally anoxic environment, with detrital pyrite and surficial uraninite common, and banded iron formations indicative of high concentrations of dissolved Fe(II) in the oceans (, ). Part of the resolution to this paradox may lie in changes in the available sinks for atmospheric O2, particularly changes in the redox state of volcanic gases that are a significant sink of produced O2 (, ).

The redox state of volcanic gases in mantle-derived melts is set by that of the mantle source. In contrast to earlier findings [e.g., (, )], recent work focusing on carefully vetted, well-characterized, and geographically widely distributed sample suites, spanning ages from ca. 3.5 Ga to recent, has documented resolvable changes in fO2 over the GOE, with an average increase of ~1.3 fO2 log units (Fig. 1). This finding, although debated, is based on two different sample types (komatiites and eclogite xenoliths) and two different approaches [partitioning of the redox-sensitive element V between olivine and komatiite melts (, ) and forward modeling of V/Sc in picritic melts (, )]. The data (including outliers) exhibit a statistically significant shift over the 2.4 Ga boundary: A two-sample Kolmogorov-Smirnov test performed on the data populations >2.4 and <2.4 Ga returns a P statistic of 1.750 × 10−5, much lower than P = 0.01, and therefore, the two populations are statistically different at a confidence of 99%. The magnitude of the ƒO2 change is sufficient to explain the GOE (), and its tempo has been demonstrated to also account for the timing of the GOE (, ). However, the underlying causes for the postulated secular change of mantle ƒO2 and its timing remain unclear.

Fig. 1.

Documented changes in mantle fO2 over time, shown in log units relative to the fayalite-magnetite-quartz buffer (ΔFMQ).

Data from Aulbach and Stagno () and Aulbach et al. () are for melts derived by ambient convecting mantle, as sampled by mantle eclogite xenoliths, orogenic eclogites, and ophiolites. Data from Nicklas et al. (, ) are for komatiitic and picritic samples through time. Both datasets are normalized to the ƒO2 of modern ambient mantle [mid-ocean ridge basalt (MORB)], as derived by the respective method (whole-rock V/Sc versus crystal-komatiite melt partitioning of V). The timing of the GOE is shown as a vertical bar. The dashed line represents an exponential mixing trend fit to the data (fit based on fluid mixing experiments, early oxidized outliers not included; see section S6 for detail) ± 1 SD.

The observed change in mantle redox state may be driven by either subduction of oxidized surface material or the by the mixing of deep oxidized mantle (–). A subduction mechanism requires the recycling of oxidized surface components, necessitating an existing oxygenated atmosphere due to, for example, burial of organic carbon or loss of H to space (, ). However, the increase in mantle fO2 predates surface oxygenation, and recycling of oxidized surface material would take too long (on the order of billions of years) to effect changes in the mantle () by the GOE. The second mechanism implies a deep oxidized mantle layer. Such a layer is believed to have formed during Earth’s magma ocean stage, due to the disproportionation of FeO during bridgmanite formation, and loss of Fe0 to the core (). This mechanism is favored for the mantle redox shift (, ). The challenge for the deep oxidized mantle mechanism is one of timing; the reservoir forms in an early Earth magma ocean and yet does not manifest as an atmospheric redox change until 2 Ga later. Archaean convection may have been slowly mixed due to the bridgmanite phase transition (), delayed basal magma ocean cooling (), or slow tectonic rates (). This timing is the issue that we quantitatively address here, implementing a new, composition-dependent rheological model for the lower mantle.

The disproportionation of FeO to create an oxidized layer implicitly increases the amount of bridgmanite in the layer; following Williams et al. (), the simplified reaction proceeds asBridgmanite [(Mg,Fe)SiO3] is the major mineral in the present-day lower mantle (~79%), followed by ferropericlase/magnesiowüstite [(Fe,Mg)O]. The disproportionation reaction shown in Eq. 1 will deplete the mantle of FeO and enhance the concentration of Fe3+-bearing bridgmanite, with free Fe0 lost to the core in the magma ocean stage.

Variations in the mineral assemblage affect density and viscosity and thus lower mantle flow. Gu et al. () suggested that a reduced lower mantle material may contain Al2O3 and be up to 1.5% denser than oxidized mantle, resulting in—problematically—rapid early mixing. However, bridgmanite is roughly three orders of magnitude stronger than ferropericlase (), and varying the proportion of these minerals will have a large effect on mantle viscosity. Ballmer et al. () used simple imposed viscosity variations to argue for long-lived bridgmanite enriched domains in the lower mantle. Here, we develop this idea by constructing a self-consistent lower mantle viscosity model, based on the constrained mineral physics of bridgmanite and ferropericlase, and integrate this with a mixing model to calculate effective viscosity.

RESULTS

Adopting the approach of Yamazaki and Karato (), we use refined solidus curves and updated mineral physics constraints, to construct diffusion profiles for bridgmanite and ferropericlase, and use a homologous scaling assumption to calculate viscosity (see section S3 for details). We then apply laboratory-constrained mixing laws to calculate the effective viscosity of this composite with depth. Our results are shown in Fig. 2.

Fig. 2.

Calculated homologous viscosity profiles for bridgmanite [(Mg,Fe)SiO3] and ferropericlase/magnesiowuestite [(Mg,Fe)O].

The variation in viscosity for these minerals is over five orders of magnitude in the lowermost mantle. We show the effective viscosity profiles (dashed lines) for mixtures of (right to left) 5, 10, 15, 20, 25, and 30% MgO, respectively. The kink in the profiles is due to a documented iron spin transition in ferropericlase/magnesiowuestite admixture (where MgO is 75% and FeO is 25%). The observed viscosity curves of Forte and Mitrovica () and Rudolph et al. () (the latter model was for Vs to density conversion factor R = 0.4, and spherical harmonic degrees l = 2 to 3), based on glacial rebound and geoid constraints, which document similar features, are also shown.

We implemented this viscosity model into the community convection code Aspect () and explored the variation in mixing rates of a bridgmanite-enriched oxidized layer under early Earth conditions (Fig. 3). Mixing rates on the early Earth have been suggested to be sensitive to convective vigor, or Rayleigh number, of the system, as well as tectonic state. To encompass these sensitivities, we developed equilibrated initial conditions at different times in Earth’s early history (4.5, 3.5, and 2.5 Ga), under a variety of tectonic scenarios (average surface velocities of 0, 1, and 10 times present plate motions). We also vary the thickness of the bridgmanite-enriched layer (0 to 500 km), its bridgmanite content (70 to 100%), the sharpness of the compositional transition to background mantle, and the mixing law exponent (which governs the smoothness of the viscosity changes with % bridgmanite; details for these models are provided in the Supplementary Materials). This approach allows us to develop controlled experiments to determine mixing rates under conditions relevant to the early Earth, with the caveats that we are not developing specific Earth history models, and mixing rates are conservative based on the mature initial conditions. The results for calculated mixing rates are shown in Fig. 4.

Fig. 3.

Evolution of mantle composition through time for two mantle convection models.

The top row (A to C) shows 300 Ma of model evolution for a simulation with a passive mantle, i.e., the bridgmanite-enriched layer has the same viscosity as the rest of the mantle. The second row (D to F) shows a model with an initial, deep, and highly viscous bridgmanite layer, which mixes more slowly and has not reached equilibration by 300 Ma. The initial thermal state is imported from a previous simulation (4.5 Ga/v = 0 case; fig. S5), which equilibrated at these conditions. Equivalent bridgmanite concentrations (F) reach 95% in the enriched layer, which grades into the background mantle (Bdg = 79%; note that the initial condition colorbar now shows model domain). The upper mantle is modeled by an olivine rheology, but we track Fe3+/ΣFe as equivalent bridgmanite concentrations in the enriched domains. The compositional distinction is tracked throughout the whole mantle but only affects the lower mantle rheologically as a modal bridgmanite contrast.

Fig. 4.

Evolution of simulated mixing trends through time.

(A) Simulated mixing curves from two models (Fig. 3) showing initial layer overturns and long-term exponential mixing trends, with either no bridgmanite-viscosity increase (passive simulation 2) or an increase in viscosity with bridgmanite content (simulation 1). Average mixing rates are calculated using parametric curve fitting of exponential curves. (B) Calculated mixing rates for a variety of models with different thermal states and properties, surface velocity scales (V scale, the velocity multiplier on surface plate motions), bridgmanite enrichment in the oxidized layer, and mixing law exponent (see eq. S5). The light gray-shaded region encapsulates the uncertainty region of our models; the central dark dashed line is a curve through the kernel density estimate maxima at each time, ±1 SEM (dark fill). The magenta box is the inferred average mixing rate required by the GOE at 2.4 Ga. A suite of our models with high bridgmanite layer viscosity and moderate thermal conditions satisfy this mixing constraint. (C) Evolution of upper mantle Δlog(fO2), calculated from the gray bounding mixing curves in (B). The timing of the GOE (red) and the inferred mantle Δlog(fO2) shift from geological constraints (Fig. 1) are also plotted. Our best estimate mixing rates calculated from the modeling suite satisfy the observed Δlog(fO2) constraints and the delayed onset of the GOE.

For each simulation, we calculate the average mixing rate of the bridgmanite layer throughout its evolution. The geological record suggests an overall fO2 change of ~1.3 ΔFMQ (fayalite-magnetite-quartz) log units over the transition, more than enough to alter the composition of volcanic volatiles to drive the GOE (). To assess the time taken to raise the average upper mantle oxygen fugacity by ~1.0 to 1.5 log units in each simulation, we assess the mixing of bridgmanite-enriched material as a proxy for Fe3+/ΣFe. To do so, we have calculated Fe3+/ΣFe for our mantle composition using the approach of Williams et al. () (see section S8) and applied the method of Kress and Carmichael () to calculate the difference in fO2 from the estimated Fe3+/Fe2+. Considering that present-day bridgmanite, constituting ~79% of the lower mantle, has a Fe3+/ΣFe of ~0.37 (, ), this implies that a variation of 10% in bridgmanite relative abundance changes the log(fO2) by ~0.8. The total range of Δlog(fO2) over our primitive reservoirs (range of 74 to 100% bridgmanite) is ~2.034.

The mixing of a bridgmanite-enriched layer (Fig. 3) occurs in two stages: first, the catastrophic disruption of the intact layer, generally by mantle plumes, and second, the progressive mixing of this entrained material in the mantle. The concentration of the oxidized reservoir in the upper mantle in the progressive mixing phase generally follows a logarithmic relationship over time, and the mixing rates shown in Fig. 4 are an average of that over 2 Ga. Simulations performed at different times shown had different initial conditions, as did simulations with different surface velocities, part of the reason for moderate mixing rates in the fast surface velocity cases. The simulations also show a significant drop in the mixing rate with increasing bridgmanite content of the oxidized layer, which makes the layer more viscous and resistant to mixing (Fig. 2). Variation in the mixing exponent in our model, from −0.05 to −0.2, has a significant effect; more negative numbers give a greater variation in viscosity with bridgmanite content, and thus more stability, and are more appropriate for large-viscosity contrast systems ().

The most significant change in mixing rate is from the evolving thermal state of the mantle. Lower heat production and core temperatures in models of younger ages result in lower Rayleigh numbers and thus less vigorous mixing. From our calculations, we have developed mixing rate curves that encapsulate the range of simulated mixing responses and their decay through time and uncertainty, and we integrate these bounds to calculate the proportion of oxidized mantle and thus Fe3+/ΣFe and fO2. This is shown in Fig. 4C, which documents the shift in upper mantle log(fO2) calculated for these mixing curves. The lower bound mixing curves are able to replicate the delayed timing of the GOE and the change in oxidation state of the mantle across this boundary.

DISCUSSION

Our results suggest that the mixing of a highly viscous, oxidized, and enriched bridgmanite layer may have driven the GOE. This was accompanied by a secular change in mantle ƒO2 of >1 log unit (–), sufficient to explain the magnitude and timing of atmospheric oxygenation (, ). The specific history of the transition was likely protracted (), and other factors may have contributed. For example, plate tectonics is far more efficient at long-term mixing of geochemical reservoirs, and its inception may have promoted the destabilization of a deep mantle reservoir (, ). A lull in global volcanic activity around 2.45 to 2.2 Ga (, ), accompanied by a decrease in volcanic degassing, could have resulted in the loss of a significant oxygen sink, permitting the rise of biogenic O2 during the GOE. The huge shift from primarily submarine volcanism before 2.5 Ga to predominantly subaerial volcanism, which is generally more oxidizing, in the Proterozoic () may have further reduced the volcanic gas sink for atmospheric O2, allowing it to rise. As discussed in (), rapid changes in atmospheric O2, such as the GOE, can ensue when tipping points are crossed even when the driving changes are gradual.

The GOE at 2.4 Ga appeared to have been preceded by a long rise in mantle fO2. This trend reflects the mixing of a deep oxidized mantle reservoir, which changed the composition of volcanic degassing products from reducing to more oxidizing, allowing oxygen levels to rise. Such a layer would have been enriched in bridgmanite, relative to the ambient mantle, and our rheological and convection models demonstrate long lag times in processing this material and raising upper mantle fO2. The scenario is potentially testable; long mixing times have been inferred previously from the persistence of anomalous signatures of short-lived radiogenic isotopes in the magmatic record that were created in the earliest Earth history (), and the melting of enstatite-enriched upper mantle has geochemical implications. The long lag in mixing an oxidized bridgmanite-rich layer may explain the long delay between low-level environmental oxidation starting at 2.8 and possibly as far back as 3.8 Ga () and the eventual rise of atmospheric oxygen at 2.4 Ga.

MATERIALS AND METHODS

The datasets used in this work are published in (–). The mantle convection simulations presented here were developed using the open-source community code Aspect (), which solves the convection equations for mass, momentum, and energy conservation, with the latter incorporating shear heating, adiabatic heating, and decaying radioactive heating sources. Our geometry is a two-dimensional axisymmetric plane, and we use adaptive mesh refinement of the computational mesh to resolve areas of large temperature gradients and near external boundaries (e.g., the surface), and our resolution in the thermal boundary layers is ~12 km.

Our rheology follows an Arrhenius form for olivine in the upper mantle, and we develop a homologous scaling approach for bridgmanite and magnesiuwuestite in the lower mantle. For the latter, the mantle viscosity η is determined from the effective diffusion Deff and by the relationshipHere, A and Ω are experimentally determined, R is the gas constant, G is the grain size, and T is the temperature. The effective diffusion itself is related to the mantle melting curve viaHere, D0 and g are, again, experimentally determined for individual minerals, and T(P)/T is the ratio of the mantle solidus to the actual mantle temperature, for different pressures (P). Our values for bridgmanite and magnesiowuestite in the lower mantle are provided in section S3, as are details of the solidus used.

The mixing of bridgmanite and magnesiowuestite is determined using the relationship of (). Here, the effective viscosity is given asHere, subscripts w and s refer to the weak and strong phases, respectively; ϕ is the volume fraction of the phase; and η is the viscosity. J is an empirical parameter, which, for a traditional Voigt model (uniform strain between phases), would be J = 1.0 or J = −1.0 for a Reuss model (uniform stress between phases). We adopt J = −0.15 as default.

We have coded our approach as a module for Aspect, and the source is available from github.org. Full details of the implementation, as well as initial and boundary conditions, are available in the Supplementary Materials.

To calculate fO2 changes due to bridgmanite composition, in a closed system constituted of only bridgmanite and ferropericlase-magnesiowüstite (Mg,FeO), we adopt the approach of Kress and Carmichael () to estimate relative differences in fO2, due to varying bridgmanite contents. Kress and Carmichael () present the relationshipThis equation assumes atmospheric pressure, and a, b, c, and d are experimentally calibrated. The summation is over the major oxide components (here, just FeO is considered). If we assume batch melting at constant temperatures, then we can calculate the difference in fO2 using changes in bulk Fe2O3 and FeO in our initial two-mineral assemblage. In this case, we can useKress and Carmichael () give a = 0.196 and dFeO = −1.828. Other parameters, benchmarking, and applications are outlined in section S8.