Molecular oxygen (O 2 ) is, and has been, a primary driver of biological evolution and shapes the contemporary landscape of Earth’s biogeochemical cycles. Although “whiffs” of oxygen have been documented in the Archean atmosphere, substantial O 2 did not accumulate irreversibly until the Early Paleoproterozoic, during what has been termed the Great Oxygenation Event (GOE). The timing of the GOE and the rate at which this oxygenation took place have been poorly constrained until now. We report the transition (that is, from being mass-independent to becoming mass-dependent) in multiple sulfur isotope signals of diagenetic pyrite in a continuous sedimentary sequence in three coeval drill cores in the Transvaal Supergroup, South Africa. These data precisely constrain the GOE to 2.33 billion years ago. The new data suggest that the oxygenation occurred rapidly—within 1 to 10 million years—and was followed by a slower rise in the ocean sulfate inventory. Our data indicate that a climate perturbation predated the GOE, whereas the relationships among GOE, “Snowball Earth” glaciation, and biogeochemical cycling will require further stratigraphic correlation supported with precise chronologies and paleolatitude reconstructions.

Keywords

Because the oldest presently known S-MDF signals are found around the Rooihoogte–Timeball Hill boundary, we measured, at high stratigraphic resolutions, multiple sulfur isotopic ratios in diagenetic pyrite in a continuous sequence from the lower Timeball Hill down to the upper Rooihoogte formations in South Africa (fig. S1). The samples in this study were collected from three diamond drill cores (KEA-4, EBA-2, and EBA-4) that penetrate the sequence in the western Transvaal Basin (see the Supplementary Materials). The Rooihoogte–Timeball Hill boundary has been dated to be 2316 ± 7 My old on the basis of a Re-Os isochron in spatially close synsedimentary to early diagenetic pyrite having the same initial 187 Os/ 188 Os ratios in one of our cores (EBA-2) ( 40 ). This is supported by a tuff zircon U-Pb age of 2309 ± 9 My for the lowermost Timeball Hill Formation in the nearby UD49 core ( 41 ) and the youngest detrital zircon U-Pb age of 2324 ± 17 My for the Timeball Hill Formation ( 42 ) ( Fig. 1 ). Therefore, these new data fill the previous gap in data and pinpoint the transition from S-MIF to S-MDF or the GOE as being present in a sedimentary rock sequence with compelling indications for continuous deposition.

One underlying problem is that there has been no consensus about the definition of GOE despite the term having been used in the literature for decades. Here, we define the GOE as the time interval during which the atmosphere changed from an anoxic state to one with sufficient pO 2 to prevent S-MIF signal preservation in sedimentary rocks. Therefore, the key questions we aimed to address in this study relate to when the atmosphere first irreversibly passed the pO 2 > 10 −5 PAL threshold (the conclusion of the GOE) and the rate at which the transition occurred (the duration of the GOE).

All isotope data are reported relative to VCDT (Vienna Cañon Diablo Troilite). The dashed lines represent three-point moving averages, whereas the shaded regions represent 1 SEM on the three points. The transitional interval highlighted in pink is defined on the basis of variation of Δ 33 S records. The stars adjacent to the KEA-4 lithologic column correspond to samples for which thin-section photomicrographs are shown in fig. S2. The Re-Os date of 2316 ± 7 Ma on syngenetic pyrite from the boundary between the Timeball Hill and Rooihoogte formations in the EBA-2 core is from Hannah et al. ( 40 ).

Stratigraphic correlation of the Rooihoogte and Duitschland Formation between the Carltonville and Duitschland area in the Transvaal Basin, South Africa. Note that the diagram is not drawn according to the real stratal thicknesses. The coarse green wavy line represents the sequence boundary in the middle Rooihoogte/Duitschland Formation that was previously proposed as the time interval in which the GOE occurred ( 34 ). Previously, the transition from S-MIF to S-MDF was obscured by an apparent hiatus in the Duitschland Formation. The pink interval represents the GOE identified in this study. The dashed vertical lines in the S-MIF box represent small Δ 33 S values (<2‰), whereas the solid lines represent the presence of large Δ 33 S values (>2‰). *From Hannah et al. ( 40 ). †From Rasmussen et al. ( 41 ).

On the basis of the transition from S-MIF to S-MDF (mass-dependent fractionation of sulfur isotopes), dating the GOE stems mainly from work carried out on sedimentary rocks from South Africa ( 18 , 34 , 35 ). Although data from the Huronian Supergroup in Canada and the Hamersley and Turee Creek groups in Western Australia also provide information on the timing of the GOE ( 36 , 37 ), the age constraints from these sedimentary sections are less precise, roughly between ~2.2 and ~2.4 Ga. However, there have been large gaps in sedimentary sulfur isotope records in South Africa ( Fig. 1 ). The original work by Bekker et al. ( 18 ) was based on S-MDF signatures in pyrite from the uppermost Rooihoogte and the lowermost Timeball Hill formations (~4-m interval straddling the Rooihoogte–Timeball Hill boundary) in the EBA-2 core ( Figs. 1 and 2 ) and suggested that the GOE should be older than 2316 ± 7 million years ago (Ma). There is a gap of ~150 million years (My) between these S-MDF signatures and the youngest S-MIF signals found in pyrite in the Brockman Iron Formation in Western Australia ( 1 , 38 ) and a gap of ~100 My between S-MIF in the Koegas Formation in the Griqualand West Basin in South Africa ( 39 ). Subsequently, Guo et al. ( 34 ) reduced this time and information gap from sulfur isotope measurements of carbonate-associated sulfate (CAS) and pyrite in the Duitschland Formation ( Fig. 1 ). However, there remained a large uncertainty about the timing and exact nature of the GOE as a major sequence boundary and a gap of about 300 m for the shale succession of which no unweathered materials are available in outcrop, which separates the S-MIF signals obtained in the lower Duitschland Formation from the S-MDF signatures obtained in the uppermost part of the upper Duitschland Formation ( Fig. 1 ). Thus, there is a pressing need and an opportunity to gather continuous S-MIF records to provide a more complete picture of the timing, tempo, and consequences of the GOE.

Mass-independent fractionation of sulfur isotopes (S-MIF) is thought to originate from the photolysis and/or photoexcitation of volcanogenic SO 2 by ultraviolet light in an anoxic atmosphere ( 1 , 20 – 23 ). Although the specific fractionation mechanisms are still under investigation ( 20 – 22 , 24 ), the atmospheric origin of S-MIF signals, such as the Δ 33 S proxy (see definition in Materials and Methods), is supported by the widespread occurrence of anomalous sulfur isotope fractionation in sulfate and sulfide minerals in Archean and early Paleoproterozoic sedimentary rocks and their absence from younger strata ( 25 – 30 ). Furthermore, results of one-dimensional (1D) photochemical modeling suggested that atmospheric pO 2 must have been less than 10 −5 times the present atmospheric level (PAL) for the production and preservation of S-MIF signals ( 31 ). Although a previous study by Zahnle et al. ( 32 ) suggested a slightly lower pO 2 (~10 −7 PAL), the link between S-MIF and the absence of atmospheric O 2 remains strong. By accepting the original value of ~10 −5 PAL as the nominal threshold above which S-MIF is no longer preserved, we apply it here as a robust tool for tracing the first rise in pO 2 to >10 −5 PAL because of its potential for high stratigraphic resolution ( 1 , 31 , 33 ).

The initial episode of demonstrable pO 2 increase is generally termed the Great Oxygenation Event (GOE) and was first identified by the disappearance of oxygen-sensitive detrital minerals such as pyrite, uraninite, and siderite from sedimentary rocks together with the widespread appearance of Fe-rich paleosols in geological sequences ( 3 , 13 ). A paucity of unambiguously continuous sedimentary records around this time meant that the GOE was only roughly constrained to the interval of 2.45 to 2.2 Ga ( 3 , 14 ). Geochemical data, including redox-sensitive trace element concentrations and isotope compositions, as well as the isotope composition of sulfide minerals that may have higher stratigraphic resolutions, help unravel the history of atmospheric pO 2 and the timing of the GOE ( 1 , 7 , 8 , 15 – 17 ). However, caution is warranted because these geochemical proxies likely have different sensitivities to changes in pO 2 and may not respond uniformly ( 1 , 7 , 8 , 18 , 19 ).

Earth’s present atmosphere is notable for its remarkably high concentration of ~21% by volume of molecular oxygen (O 2 ) upon which all complex life depends, whereas geologic and geochemical constraints suggest that before ~2.45 billion years ago (Ga), the atmosphere was essentially devoid of this gas ( 1 – 4 ). The rise of the atmospheric O 2 partial pressure (pO 2 ) is proposed to have occurred in two main episodes, one at the beginning and the other at the end of the Proterozoic Eon ( 3 , 5 , 6 ). More recent studies, most of which were driven by the emergence of redox-sensitive trace element proxies, suggest a more complex and dynamic evolution scenario for pO 2 ( 7 – 11 ). Therefore, continuous records and precise geochronological constraints are absolutely critical to defining and understanding these patterns in relation to the myriad of underlying biological and geological drivers ( 12 ).

Both the pyrite in bulk rock and that drilled from layered or granular pyrite show similar sulfur isotopic compositions (see Materials and Methods, Fig. 2 , and table S1). On the basis of the time series variation trends of Δ 33 S and the relationships between Δ 33 S and Δ 36 S values (see Materials and Methods for the definition of Δ 36 S), we have divided the sequence in the three cores into three intervals, termed S-MIF, transitional, and S-MDF ( Fig. 2 ). The pyrite in the S-MIF interval (below the pink band in Fig. 2 ) carries distinct S-MIF signals, with Δ 33 S reaching as high as +8‰ and the Δ 36 S/Δ 33 S ratios centering around −0.98, which is characteristic of the Neoarchean array ( Fig. 3 ) ( 30 , 43 ). In contrast, in the S-MDF interval (above the pink band in Fig. 2 ), all pyrite is characterized by near-zero Δ 33 S values (|Δ 33 S| < 0.3‰) and Δ 36 S/Δ 33 S ratios close to the theoretical MDF value (approximately –6.7) ( 44 ). The intervening transitional interval (the pink band in Fig. 2 ) contains pyrite with both small (|Δ 33 S| < 0.5‰ in the lower part) and large (0.5 to 2.0‰ in the upper part) Δ 33 S values and a Δ 36 S/Δ 33 S slope that is higher (less negative) than both the underlying S-MIF interval and the overlying S-MDF interval. Thus, the evolution from the S-MIF to the transitional interval is characterized by a sharp negative shift in Δ 33 S and a significant increase in the slope of Δ 36 S/Δ 33 S, whereas the changes from the transitional to the S-MDF interval are characterized by the disappearance of S-MIF signals and a decrease in the slope of Δ 36 S/Δ 33 S.

Multiple lines of evidence show that the pyrite grains analyzed in this study are authigenic. First, the sizes of the pyrite grains are much larger than coexisting terrestrial particles, mainly comprising quartz and clay minerals (fig. S2, A to I). A detrital origin for the pyrite is not feasible on the basis of hydraulic equivalence because the density of pyrite is higher than that of other minerals. Second, the size range of the pyrite grains does not change according to that of the matrix (fig. S2F). Third, some pyrite grains distinctly disturb sedimentary laminae (fig. S2, B to E and G to I), suggesting that the crystals grew in place before compaction and, therefore, are early diagenetic in origin. Furthermore, pyrite morphology and size do not change across the rock sequence studied here (fig. S2, A to C versus D to I).

DISCUSSION

The Δ33S values of pyrite in the S-MIF interval (up to +8‰) are among the highest values in the entire geologic record [for example, Johnston (27)]. These data are significantly larger than published data in the interval between ~2.5 and ~2.3 Ga that are commonly located in the range of −1 to 3‰ (25, 27). This may partially reflect sample bias because only a few geologic formations that contain rock types suitable for sulfur isotope analysis are currently known. As a result, data are nearly exclusively from the Transvaal Supergroup in South Africa (34), the Hamersley and Turee Creek groups in Western Australia (36), and the Huronian Supergroup in Canada (37). The large Δ33S values might be due to some particular facet of the atmospheric chemistry at that time (22, 45). The high, spatially heterogeneous, and temporally variable Δ33S values and the large area of the Transvaal Basin where these rocks were deposited suggest that the photochemical products of SO 2 in an anoxic atmosphere were the primary origin of these signals. Recycling of continental sulfides with S-MIF signals, which has previously been used to interpret the small Δ33S values in this time interval (25, 46), cannot explain these new data (see the Supplementary Materials). This is also supported by the Δ36S/Δ33S slope of −0.9, which is widespread in Neoarchean sedimentary rocks and interpreted to be of atmospheric origin (9, 25, 30, 43). Accordingly, our data indicate that the atmosphere was effectively devoid of oxygen (nominally <10−5 PAL) when the lower part of the upper Rooihoogte Formation was being deposited (Figs. 2 and 3).

The near-zero Δ33S values and mass-dependent Δ36S/Δ33S signature slope (approximately −6.7) in the S-MDF interval indicate that mass-independent fractionation did not take place or was not preservable during the deposition of the uppermost part of the Rooihoogte and lower Timeball Hill formations. This is entirely consistent with previous research that only investigated a narrow stratigraphic interval straddling the Rooihoogte–Timeball Hill boundary and the top Duitschland Formation (18, 34). Thus, our data indicate that the atmosphere had a substantial pO 2 (>10−5 PAL) following the deposition of the uppermost part of the Rooihoogte Formation. Although mechanisms other than pO 2 , such as decrease of pCH 4 , have been proposed to explain the disappearance of mass-independent fractionation signals (32), geochemical cycles of CH 4 and O 2 are intimately linked, as we show below. A significant increase in marine sulfate levels in this interval (see below), widespread hematite-rich oolitic ironstone, and a high enrichment of chromium in the middle part of the Timeball Hill Formation (16, 47) all suggest that a substantial rise in pO 2 was at play.

Meanwhile, our data also indicate that before the substantial rise in pO 2 , the atmosphere underwent significant compositional changes indicated by the sharp decline in the magnitude of Δ33S values and dynamic fluctuations in the Δ36S/Δ33S slope observed in the pyrites of the transitional interval (Figs. 2 and 3). Because the Δ36S/Δ33S ratio in the transitional interval cannot be explained by the mixing of signals from MIF and MDF intervals (see the Supplementary Materials), we suggest that the decrease in Δ33S values in the transitional interval was caused by an increase of pO 2 , which weakened the preservation of Δ33S signals from photochemistry through mixing different atmospheric sulfur pools; however, this homogenization was incomplete (9, 25). Some samples in this interval in all three cores that have near-zero Δ33S values lend support to this hypothesis. This suggests that pO 2 started to increase slightly in the transitional interval. Decreases in the magnitude of Δ33S values have also been ascribed to the shielding effect of organic haze, such as the small Δ33S values present in the Mesoarchean (9, 37, 48). This effect would also decrease the Δ36S/Δ33S slope (9). For example, the Mesoarchean Δ36S/Δ33S slope (approximately −1.5) is lower (more negative) than that of the Neoarchean (approximately −0.9). The higher Δ36S/Δ33S slope in the transitional interval suggests that the effect of organic haze might not be the key factor for the decline in Δ33S values (Fig. 3). Therefore, we hypothesize that this change in atmospheric composition was the prelude to the substantial increase in pO 2 , representing the GOE as defined above.

The high stratigraphic resolution and continuous record of the transition from S-MIF to S-MDF signals in these three cores speak to dynamic processes as underpinning the substantial rise in pO 2 and provide an opportunity to precisely constrain the time of pO 2 rise at the GOE to >10−5 PAL, as well as the duration of the transition. The package of sediment between where S-MIF signals cease to be discernible (that is, an oxidized atmosphere) and the Rooihoogte–Timeball Hill boundary is approximately 5 m thick in all three cores (Fig. 2). Sequence stratigraphic analysis indicates continuous deposition from the upper Rooihoogte Formation to the lower Timeball Hill Formation (49) and is supported by uniformly fine sedimentary rock (mudstone, including black shale and siltstone, and few carbonates) throughout (see the Supplementary Materials). Although with some uncertainty, a conservative sedimentation rate can be assumed on the basis of previously reported sedimentation rates in slightly older parts of the same basin—rates that are also similar to modern values (50). The decompacted sedimentation rate of modern siliciclastic sediments, except for deep-sea clay, is >5 m/My (51). Assuming 80% compaction for shale and 50% of time represented by deposition, the duration of the 5-m package is less than 10 My. Applying this estimate to the age of the Rooihoogte–Timeball Hill boundary (2316 ± 7 Ma), our data constrain the termination of the GOE to be within 2316 to 2326 Ma (±7 My).

The sedimentary package corresponding to the transitional interval, which represents the duration of the GOE, is thin. It spans about 5 m of strata in all three cores (Fig. 2), which, according to the conservative sedimentation rate assumed above, represents <10 My. If we assume a more rapid sedimentation rate (about 50 m/My) of siliciclastic sediments in a foreland basin as the Rooihoogte Formation deposited (52) and a lower compacted ratio of the siltstone/mudstone relative to black shale, the transitional interval would represent less than 1 My in time. In either case, the GOE took place on a geologically rapid time scale. Such a fast rise of pO 2 is consistent with atmospheric modeling results that predict a rapid and irreversible transition from stable anoxic to oxic atmospheres (53, 54) and may provide constraints on the controversial problem of whether the rise of O 2 during the GOE was primarily a consequence of an increase in the flux of oxygen sources or a decrease in the flux of oxygen sinks (34, 55–57).

The GOE not only ended the ultraviolet photoreaction of SO 2 in the troposphere but also had effects on exogenic sulfur cycling. A large magnitude (about 45‰) of negative shift in δ34S occurred in the lower part of the S-MDF interval (Fig. 2). The pyrite around the Rooihoogte–Timeball Hill boundary in these three cores is depleted in 34S, with δ34S values as low as −35‰. Although minor occurrences of low δ34S pyrite have recently been reported in late Archean strata by in situ analysis (58), samples studied here are the oldest samples that exhibit low δ34S values in bulk rock analyses (18, 59). The extremely 34S-depleted pyrite suggests that microbial sulfate reduction occurred under sulfate-unlimited conditions in the upper portion of the MDF interval (60). In contrast, high δ34S values of +2.4 to +24.4‰ further downcore in the S-MIF and transitional intervals are consistent with quantitative sulfate reduction under low oceanic sulfate levels (60) before the GOE. Therefore, these data suggest that marine sulfate levels increased considerably following the GOE. An increase in the marine sulfate reservoir at about 2.3 Ga is also supported by the large sulfur isotope separation between CAS and pyrite (attributed to microbial sulfate reduction) recorded in the topmost Duitschland Formation (34), as well as the recent finding of gypsum in sabkha-facies evaporites at ~2.31 Ga (41).

The negative shift in δ34S values occurred above the disappearance of S-MIF in all three cores by ~3 m of strata (Fig. 2), representing as much as ~6 My in time on the basis of the sedimentation rate assumed above. Thus, the increase in seawater sulfate lagged the rise in pO 2 by ~6 My. An elevated sulfate flux into the ocean through oxidative pyrite weathering is an expected consequence of higher pO 2 after the GOE and can partly contribute to the rise in the oceanic sulfate concentration. It has been suggested that the threshold of pO 2 required to induce high pyrite weathering rates is the same as or a little lower than that required for the preservation of S-MIF signals (7, 43, 61). Therefore, the rise in seawater sulfate levels should coincide with or predate the disappearance of MIF. The observed time lag and opposite sequence of events suggest that other factors, not only the changes in sources (for example, sulfide weathering), contributed to the expansion of the oceanic sulfate reservoir, possibly including variation in the sinks of the seawater sulfate.

Geochemical data and mass balance analysis suggest that sulfate sinks must have been relatively much stronger before the GOE than in the modern ocean. As shown in table S2, in the modern sulfur cycle, volcanic sulfur input can comprise as much as 10% of total sulfur in the ocean-atmosphere system, which is dominantly from pyrite and gypsum weathering. Very low pO 2 before the GOE would severely limit the amount of sulfur input from weathering processes. Therefore, sulfur inputs into the pre-GOE ocean-atmosphere system are dominated by volcanic fluxes, which can be one to three times higher than the modern value, considering the likely higher heat flux and the associated higher volcanic gas fluxes (53). Thus, assuming that the sulfate reduction sink is first order with respect to sulfate concentration and that residence times were similar, expected sulfate levels before the GOE would be at least 10% (and perhaps as much as 30%) of today’s values (2.8 to 8.4 mM). This is not consistent with models that suggest extremely low seawater sulfate levels before the GOE, possibly as low as 2.5 μM (34, 62). Therefore, a much higher sink of seawater sulfate is required to maintain the very low sulfate levels previously suggested (8, 34, 60). We describe below how pyrite burial through microbial sulfate reduction modulated by pO 2 can simultaneously account for the low seawater sulfate levels before the GOE and the delayed rise to high sulfate levels after the GOE.

To discern the sequence of carbon, sulfur, and oxygen cycle perturbations, we constructed a model of biogeochemical cycles. This model integrates the sulfur cycle model to a coupled atmospheric-biogeochemical model of CH 4 and O 2 (53, 54). In the sulfur cycle model, continental sulfide weathering and volcanic sulfur were the two main sources of seawater sulfate, whereas pyrite burial and sulfate burial were the two main sinks (see Materials and Methods and table S2). We assumed that sulfide weathering was coupled to pO 2 , and volcanic sulfur fluxes were scaled to the outgassing rate of reducing gases (63). Sulfate burial is proportional to seawater sulfate content. Pyrite burial is assumed to be governed by microbial sulfate reduction dependent on sulfate content and on pO 2 to account for the influence of quality and availability of organic matter to sulfate-reducing microbes.

Our model both reproduces the rapid rise of oxygen (53, 54) (consistent with the rapid disappearance of S-MIF in the cores we analyzed here) and predicts a time lag of ~10 My between the increase in pO 2 and the subsequent rise in sulfate from ~10 μM to >1 mM (Fig. 4); this is comparable to the lag in isotope signals estimated above. Although the elevated pO 2 enhances the oxidative weathering of sulfide minerals on land, thereby increasing sulfate flux into the oceans, the modeled increase in sulfide weathering flux alone cannot explain the increase in sulfate levels by approximately two orders of magnitude (Fig. 4). Rather, the required increase in seawater sulfate can be simulated if the flux of pyrite burial decreases in response to higher pO 2 (Fig. 4). This behavior can be accounted for by a reduced supply of labile organic matter due to enhanced aerobic remineralization of organic matter as a consequence of the rise of O 2 (see “Model description” in Materials and Methods and the Supplementary Materials).

Fig. 4 Model results showing the increase in pO 2 and seawater sulfate level across the GOE. The red line represents the atmospheric pO 2 . The blue solid and dashed lines represent seawater sulfate levels modeled with the pyrite burial flux coupled or decoupled with pO 2 , respectively.