Mechanistic links between intense volcanism and the transient oxygenation of the Archean atmosphere

Introduction

The evolution of the Earth’s atmosphere plays a crucial role in maintaining long-term habitability and shaping the coupled evolution of life and its environments. One of the most important aspects of the evolution of Earth’s atmosphere was the first substantial increase in the partial pressure of atmospheric oxygen (pO₂) at ~2.5–2.2 Ga, known as the Great Oxidation Event (GOE) (Fig. 1a)^1,2,3. Before the GOE, the atmospheric pO₂ was much lower than the present atmospheric level (<10^–6 PAL), as constrained by geological records of mass-independent isotope fractionation of sulfur (S-MIF)^1,2,4,5. Across the GOE, atmospheric pO₂ rose to ~10^–4–10^–1 PAL (Fig. 1a)^1,6,7,8,9,10.

Mechanistic links between intense volcanism and the transient oxygenation of the Archean atmosphere — **Fig. 1: Earth system evolution over the Archean-Paleoproterozoic.**

Despite the protracted reducing conditions before the GOE, redox-sensitive element, and isotope records indicate the occurrence of transient oxygenation events in the ocean–atmosphere system, known as whiffs of oxygen events^{11,12,13,14,15,16,17,18,19,20,21,22,23} (Fig. 1). The most well-studied whiff is recorded in Mt. McRae Shale (Hamersley Basin, Australia) deposited at ~2.5 Ga (Fig. S1)¹¹. Enrichment of molybdenum (Mo) and rhenium (Re) in the shale indicates a transient increase in atmospheric pO₂ because Mo and Re are sourced from the oxidative weathering of continental sulfide minerals^11,21. This view is further supported by the enrichment of selenium (Se), which also indicates the oxidative weathering of continental crusts¹⁹. During the whiff intervals, the marine sulfate (SO₄^2–) concentrations likely increased owing to the oxidative weathering of sulfide minerals, promoting the development of local euxinia on continental margins^23,24. Additionally, the osmium isotopic records in Mt. McRae Shale indicate a transient increase in oxidative weathering rate of continental sediments²¹. This whiff event would have persisted for several millions years to ~11 Myr (Fig. S1)^11,23. During this event, a positive shift in δ¹⁵N suggests that the ocean was at least locally oxygenated^18,23,25. Although the occurrence of whiffs is subject to vigorous debate^26,27,28, the whiff event recorded in Mt. McRae Shale was also coeval with redox-sensitive elements enrichment in the Klein Naute Formation (South Africa)^23,29, reinforcing the idea of the whiff of oxygen. Furthermore, recent phylogenic estimates place the emergence of oxygen-utilizing enzymes long before the GOE, supporting the presence of O₂ during the late Archean³⁰.

Large igneous provinces (LIPs), voluminous volcanic rocks formed by massive volcanic events, are a potential driver of global transient changes in the biogeochemistry of the ocean–atmosphere system³¹. The activity of LIPs is linked to upwelling mantle plumes that cause the extensive emplacement of magma³², leading to global warming through the release of vast amounts of greenhouse gases (e.g., CO₂) into the atmosphere. It has been discussed that global warming associated with LIP eruptions during the Phanerozoic caused oceanic eutrophication by elevating the riverine input rate of nutrients. This, in turn, would have promoted primary productivity and O₂ production in the euphotic zone, leading to excess organic matter production. The subsequent decomposition of organic matter would have caused oceanic deoxygenation due to the elevated oxygen demand^{33,34,35,36,37,38,39,40}. Finally, the excess organic matter that was not decomposed in the water column would be buried in marine sediments, which is equivalent to the net oxygen production in the geological timescale, possibly causing the oxygenation of the atmosphere. This sequence of events is thought to have occurred repeatedly throughout the Phanerozoic (~0.54–0 Ga), sometimes causing biotic crises⁴¹.

Given the increase in biological productivity and O₂ production in the euphotic zone during oceanic eutrophication, repeated LIP eruptions during the late Archean³² may have contributed to the occurrence of whiffs (Fig. 1b). The timing of the period of intense subaerial volcanism and the interval of the whiff^22,42,43, together with the abundance of mercury (Hg) and its isotopic signatures, suggest that subaerial volcanism preceded the interval of the whiff²², implying close association between these phenomena. However, the biogeochemical consequences of LIP volcanism on the reducing environment on Archean Earth and its mechanistic link to the whiffs occurrence remain ambiguous due to the lack of a quantitative framework that allows for a comprehensive investigation of the coupled dynamics of the ocean–atmosphere–biosphere–climate during the Archean. Additionally, the impact of the long-term evolution of the Archean environment (e.g., continental growth) on the frequency and magnitude of the whiffs remains elusive.

Here, we developed a theoretical model of the global C–P–S–Fe–O₂–Ca biogeochemical cycles (Fig. 2; see Materials and Methods) to support the hypothesis that whiffs could be reasonably explained as the result of LIP eruptions. This model was based on the model used in the previous study⁴⁴, in which we introduced the dynamical S cycles and C and S isotopes (δ¹³C and δ³⁴S). The ocean–atmosphere system is represented by an atmospheric box, four ocean boxes, and a pore-space water box in oceanic crust (Fig. 2), with the marine S cycle represented by a single ocean box. The model simulates variations in atmospheric pCO₂, pCH₄, and pO₂. In the ocean boxes, oceanic levels of dissolved inorganic C (DIC), carbonate alkalinity (Alk_c), [PO₄^3–], [SO₄^2–], [Fe²⁺], [Ca²⁺], dissolved O₂, and dissolved CH₄. The behavior of DIC, Alk_c, and [Ca²⁺] are also considered in the pore-space water box. Using this model, we simulate a transient response of the system to the injection of CO₂ and reducing gases into the atmosphere assuming intense LIP-induced volcanism. Although the box modeling approach does not resolve spatial heterogeneity within the ocean for elements with short residence times, it is sufficient to simulate the long-term response of the biogeochemical cycle to the LIP-induced volcanism. This model was also used to quantify the conditions suitable for whiff occurrence on Archean Earth. The results demonstrate that the susceptibility of atmospheric O₂ levels to LIP eruptions is closely related to continental growth.

**Fig. 2: Schematic illustration of the biogeochemical model.**

Results

Eruption of a LIP causes a whiff of oxygen

To assess the biogeochemical impacts of LIP eruptions on the Archean Earth system, we ran the model under LIP eruptions with different total C influxes (Fig. S2). The maximum total C influx was set based on the case of intense volcanism that formed the Ontong Java Plateau, which is the largest Phanerozoic LIP known to have erupted in the Pacific Ocean during the Cretaceous. The eruption of the Ontong Java Plateau is known to have caused oceanic eutrophication and anoxia, and thus such a large-scale LIP eruption would potentially cause a whiff. We adopted a maximum estimate of the total CO₂ influx from the Ontong Java Plateau⁴⁵ as the maximum total C influx. We cannot rule out the possibility that the total C influxes from Archean LIPs were higher than this value. Nevertheless, we consider it unlikely that higher influx scenarios would dramatically alter our overarching conclusions. The intense volcanism is assumed to have continued for a period of 400 kyr in the standard case, which is also based on the estimated value for the case of the Ontong Java Plateau³⁶, while the sensitivity of the results to different values of the duration of the volcanism (100–1000 kyr) is also tested. We also assumed a subaerial LIP eruption (i.e., CO₂ and reducing gasses are injected directly into the atmosphere). Although a submarine eruption would have been more likely to occur during the Archean⁴³, this assumption has only a minor impact on the results because the timescale for reaching equilibrium of gas exchange between the ocean and the atmosphere would be at most the timescale of the oceanic circulation (<~1000 years), which is sufficiently short.

The result demonstrates that during LIP eruptions (yellow-hatched region in Fig. 3d), atmospheric pO₂ levels slightly decrease from their initial values (pO₂ ~ 10^–8 PAL) owing to the influx of reducing gases into the atmosphere. After the cessation of LIP volcanism, atmospheric pO₂ begins to increase as oceanic primary productivity increases. With a sufficient influx of total C, atmospheric pO₂ exceeds the threshold value necessary for Mo supply from oxidative weathering (10^–6.9 PAL), which corresponds to whiffs⁴⁶. Specifically, when the total C influx is larger than ~20% of the maximum value, atmospheric pO₂ exceeds the S-MIF limit and reaches ~10^–4–10^–2 PAL, which persists for several to ten million years (Fig. 3d). The atmospheric pO₂ increases rapidly to this level owing to positive feedback caused by the self-shielding effect of O₂ by its byproduct O₃⁴⁷. This atmospheric pO₂ is also beyond the threshold value necessary for oxidative weathering (10^–6.9 PAL), which would sufficiently explain the enrichment of redox-sensitive elements that are hosted in sulfides during the whiffs interval^12,24,46. The primary productivity required to trigger the whiffs (pO₂ > 10^–6.9 PAL) is ~270 Tmol C yr^–1 (Fig. 3c), which is comparable to estimates required for triggering the GOE^44,48. Increases in primary productivity are driven by elevated marine P concentration caused by enhanced supply of riverine P due to intensified continental weathering (Fig. 3b).

**Fig. 3: Response of the atmospheric and oceanic biogeochemical cycles caused by an intense eruption of a large igneous province (LIP) (*Exp_std* in Table S5).**

For the case with maximum C influx (red line in Fig. 3), atmospheric pCO₂ level increases to ~230 PAL during intense volcanism. Following the cessation of volcanism, pCO₂ gradually decreases over a period of >10 Myr. This timescale is relatively long compared with the perturbations caused by LIP eruptions during the Phanerozoic, which typically cease within several million years^{36,49,50,51,52,53,54,55,56,57,58,59,60}. This difference is caused by the low terrestrial weatherability during the Archean, i.e., the absence of land plants and small continental areas diminished the continental weathering efficiency^61,62. In the Archean, weathering of seafloor basalt was likely the primary driver of climatic adjustment via carbonate–silicate geochemical cycles^62,63,64,65. Despite the slow recovery of atmospheric pCO₂, the surface temperature declines soon after the cessation of the LIP volcanism owing to the marked reduction in atmospheric pCH₄ during the interval of whiffs. Following the cessation of LIP volcanism, marine P concentrations also decrease sharply, because the whiff of oxygen and oceanic oxygenation suppress the recycling of P from marine sediments (Fig. 3b). Marine P concentration then gradually declines following its removal via deposition of organic matter. Once marine P concentration falls below a critical value required for atmospheric oxygenation, atmospheric pO₂ level decreases rapidly, returning to the Archean-like levels, consistent with predictions from a photochemical model⁶⁶.

The duration of whiffs triggered by the LIP eruption typically spans several to ten million years (Fig. 3d). This timescale is primarily controlled by the timescale of the removal of atmospheric pCO₂ via weatherings of continental and oceanic crusts. Because of the low terrestrial weatherability during the Archean, the residence time of CO₂ was long, and CO₂ injected from the intense volcanism was removed from the atmosphere via weatherings over a period of 1–10 Myr. Consequently, the rate of riverine P supply, marine P concentrations, and oxygen production rate to the atmosphere remain elevated over this timescale. The timescale for oxygen loss in the atmosphere during the whiffs would be shorter than the timescale of the excess oxygen production because it would be much shorter than the present-day residence time of ~1 Myr. Therefore, the timescale of the whiffs follows that of the P accumulation in the ocean after the LIP eruption.

Conditions for a whiff of oxygen

Our results indicate that whiffs can be triggered by the accumulation of P in the ocean following LIP volcanism, and therefore factors influencing the rate of riverine P supply would impact the conditions necessary for a whiff. The most critical factor is the areal fraction of the continent, as in the conditions of the GOE^62,64, because it is a primary factor that affects the terrestrial weatherability in the absence of land plants. Results obtained using various continental areas and total C influxes are shown in Fig. 4. With a small continental area relative to that of the present (f_a), the steady-state atmospheric pCO₂ is ~290 PAL (Fig. S3a). With high total C input, atmospheric pCO₂ increases to ~500 PAL after LIP eruptions. Despite this increase, marine P concentration changes only slightly due to low continental weathering rates, preventing whiff occurrence (Fig. 4a). As f_a increases, the total C influx from the LIP eruptions required to trigger the whiffs becomes lower, making whiffs more likely than with the standard f_a. Both the maximum pO₂ achieved following LIP eruptions and the duration of the whiff increase with increasing f_a. When f_a is sufficiently large, background atmospheric pO₂ level exceeds the S-MIF limit, corresponding to the GOE. The conditions for the GOE are affected by multiple factors, including f_a and the influx of reducing power into the system^{44,48,62,64,67} (Fig. 5a).

**Fig. 4: The conditions of continental area and total C input from a LIP required for the occurrence of the whiffs of oxygen.**

**Fig. 5: The conditions of continental area and the hydrothermal input rate of Fe(II) required for the occurrence of the whiffs of oxygen.**

The conditions for whiffs are also affected by other factors controlling the global redox budget. The required influxes of hydrothermal Fe(II) and reducing gases from volcanoes for the whiffs are summarized in Fig. 5d and Fig. S4, respectively. This result also demonstrates that as f_a increases, the range of the influx of reducing materials that allow whiff occurrence extends (Fig. 5a, b). When influxes of reducing materials decrease sufficiently, the GOE occurs. Therefore, whiffs likely indicate a gradual shift toward more oxidizing conditions in the ocean-atmosphere system during the Archean. Influxes of reducing materials from LIPs, however, do not strongly affect the conditions for whiffs (Figs. S5 and S6) because the timescale for the loss of the reducing gases from LIPs is shorter than the duration of the LIP eruptions (Fig. 3d). The condition for whiffs occurrence is also dependent on the duration of LIP eruptions (Fig. S7), but the dependence is considerably weak compared with that of the total C input because the timescale for the responses of oceanic eutrophication and the whiff of oxygen are longer than that of the LIP eruption. In summary, our results demonstrate that the total C influx from LIPs, continental area (f_a), and the influx of reducing power (i.e., Fe(II) and reducing gases) are the critical factors for triggering whiffs.

Discussion

In this study, we assessed the impact of LIP volcanism on Archean biogeochemical cycles. The results demonstrated that intense volcanism associated with LIP emplacement could have caused transient oxygenation of the atmosphere lasting ~1–10 Myr. This timescale is consistent with the value of the Mt. McRae shale (<~11 Myr)^11,23, which was estimated assuming a typical shale accumulation rate¹¹, as well as the estimated maximum value of ~50 Myr for the whiff in the Jeerinah Formation, Australia, at ~2.66 Ga¹⁶. Our model also indicated that higher f_a increases the susceptibility of the system to a whiff by increasing the P supply rate during the period with elevated atmospheric pCO₂ caused by CO₂ accumulation after LIP volcanism. During the mid-late Archean (3.5–2.5 Ga), the continental crust mass would have increased (Fig. 1d)^{68,69,70,71,72,73,74,75,76,77,78}. This might have provided the backdrop to the transient oxygenation events in the late Archean, which would have ultimately contributed to the GOE^{44,62,64,70,79,80,81}. The evolution of the continental crust was likely a critical step for the evolution of the atmosphere and biosphere during the late Archean and early Paleoproterozoic.

Even if the atmospheric pO₂ exceeds the S-MIF limit transiently, it does not necessarily mean the disappearance of the S-MIF in geological records because S compounds with S-MIF signals would be continuously transported from the continental crust until it is completely lost, which is called the crustal memory effect⁸². When a whiff occurs, our model estimates that atmospheric pO₂ level increases and could readily exceed ~10^–4 PAL owing to positive feedback among atmospheric O₂, O₃, and CH₄^47,83. In the Mt. McRae shale, on the other hand, negative S-MIF signals are observed during the interval of whiffs⁸⁴, which seems not to be consistent with the high atmospheric pO₂ simulated in our model. However, this inconsistency would be explained by the crustal memory effect because both physical and chemical weathering would have been very intense to transport the crustal S-MIF signals to the ocean under the conditions with elevated atmospheric pCO₂ that followed the LIP volcanism. The impact of the crustal memory effect following LIP volcanism should be investigated to better constrain the maximum level of atmospheric pO₂ during the whiffs. Alternatively, the S-MIF signals in the Mt. McRae shale may actually indicate that the atmospheric pO₂ level might not have exceeded the S-MIF limit. This interpretation might not necessarily be inconsistent because some records suggest an increase in the rate of oxidative weathering, even with a level of atmospheric pO₂ below the S-MIF limit^28,46. This can still be explained by LIP eruptions with small total CO₂ influxes, because the time required for returning to the background level of atmospheric pO₂ following the LIP eruptions is very long (>10 Myr) even if a maximum level of pO₂ is lower than the S-MIF limit owing to the low terrestrial weatherability as discussed above.

Although there is as yet no clear evidence linking LIP eruptions to whiffs of oxygen, the age distribution of known LIPs suggests that such eruptions were common during the late Archean (e.g., BLIP-1 at 2.505 Ga, Crystal Springs at 2.510 Ga, and Mistassini-Ptarmigan-Irsuaq at 2.51–2.50 Ga)³² (Fig. 1b). Additionally, Hg isotope signatures in the lower Mt. McRae Shale, which precedes the whiff in the upper Mt. McRae Shale, provide a temporal link between the whiff and the LIP volcanism²². The stratigraphic separation of these two signals is at least 20 m, most of which is covered by silicified ferroan carbonate²². The average sedimentation rate of the carbonate-bearing rocks in the Hamersley group is assumed to be 12 m Myr^–1 ^22,85. With this assumption, the time lag between these two signals is ~1.7 Myr. This time lag could be explained if the timescale for the eruption of a LIP is comparably long because the supply of reducing gases from a LIP prohibits the occurrence of the whiff until the end of the eruption of a LIP.

The timing of both the onset of the GOE and the permanent oxygenation of the atmosphere has been actively discussed based on the presence or absence of S-MIF signals in geological records^82,86,87,88. Indeed, there were LIP eruptions during the period of the GOE (Fig. 1b). Therefore, our result indicates that LIP eruptions could deviate the timing of the onset and the end of the GOE by causing whiffs of oxygen^82,87,88,89. In other words, the occurrence of whiffs during the late Archean indicates that the redox budget of the ocean–atmosphere system was close to the tipping point for the GOE, beyond which the abrupt rise in atmospheric pO₂ would take place.

As the system approaches the GOE tipping point, it becomes more vulnerable to environmental perturbations and variations. Consequently, transient oxidation may also result from natural fluctuations unrelated to LIP eruptions. For example, the global biogeochemical cycles can be fluctuated by variations of the orbital parameters as discussed in various ages of the Phanerozoic^{90,91,92,93,94} and/or a bolide impact event^95,96. Although these potential perturbations may have contributed to causing the whiffs, these processes typically operate on timescales shorter than 1 Myr. In this light, the LIP eruption events have an advantage in explaining the long timescale of the whiffs inferred from the geological records.

When a whiff of oxygen occurs following LIP volcanism, our model expects a negative carbon isotope excursion (CIE), followed by a longer positive CIE due to oceanic eutrophication and extensive organic C burial in the ocean. This behavior is consistent with CIEs observed after the Phanerozoic LIP eruptions^36,97. When the duration of LIP volcanism is relatively short (e.g., 400 kyr), the negative CIE caused by C accumulation during the LIP eruption would be more pronounced than the positive CIE (Fig. S7). A gradual transition to higher δ¹³C values after the whiff is recorded in the Mt. McRae shale, consistent with the model, although no clear negative CIE is observed. This can be explained by a prolonged LIP eruption, which would diminish the negative CIE while maintaining a similar positive CIE magnitude (Fig. S7). Additionally, carbonate and TOC records in the Mt. McRae shale exhibit characteristics similar to responses to Phanerozoic LIPs, suggesting global C cycle variations similar to those associated with Phanerozoic LIPs (Fig. S8). Our model further demonstrated that marine SO₄^2– concentration increases when a whiff occurs after the LIP eruption because the oxidative weathering rate increased (Fig. S9c). Concurrently, [ΣH₂S] also increases due to the enhanced microbial sulfate reduction (MSR) (Fig. S9d). Although our box model does not resolve local euxinic environments near continental margins, this is consistent with geological observations²³. The model also demonstrated that S isotope excursions would be unlikely to be observed during the interval of a whiff caused by LIP volcanism because marine SO₄^2– concentration during the whiff remains too low for the burial of gypsum (Fig. S9e). Our model also expects that increases in organic C and reactive P abundances would also be recorded during and after LIP eruptions. Further examination of these abundances during the whiffs would provide better verification of the link between LIP eruptions and whiffs occurrences expected based on this study.

The timing of the emergence of oxygenic photoautotrophs (OP), which would be a prerequisite for atmospheric oxygenation, is a subject of active debate. Some studies suggest that the emergence of OPs would have been immediately before the GOE^26,98,99, inferring that the records of redox-sensitive elements during the whiff intervals may not reflect true atmospheric oxygenation. Although the present study does not rule out this possibility, the theoretical link between LIP eruptions and whiffs of oxygen supports transient atmospheric oxygenation^{11,12,13,14,15,16,17,18,19} and the emergence of OPs more than several hundred million years earlier than the GOE³⁰. The transient atmospheric oxygenation might also have stimulated the evolution of life, including O₂ detoxification mechanisms, which might have been a critical step in adapting to the increase in the atmospheric pO₂ level across the GOE. Further investigation into the coevolution of the atmosphere, biosphere, and solid Earth during the late Archean and Paleoproterozoic would be fruitful for phylogenic, geological, and theoretical research.

Conclusion

Under a reasonable parameter set, our biogeochemical model indicates that a whiff of oxygen—a transient oxygenation event of the atmosphere during the late Archean—could have been triggered by intense volcanism caused by LIP emplacement. Our model explains the observed increases in atmospheric pO₂ levels during the whiff intervals, as documented in geological records. The duration of a whiff is estimated at ~1–10 Myr, primarily influenced by the total amount of C supplied from LIPs and the extent of the continental area. The evolution of the continental crust during the Archean likely enhanced P supply following LIP eruptions, consistent with the occurrence of the whiffs during the late Archean. This suggests that the occurrence of the whiffs indicates that the late Archean Earth system was approaching the tipping point for a permanent rise in atmospheric pO₂. The coupled dynamics of the atmosphere, biosphere, and continents, as recorded in the whiffs, would have been a fundamental step for the onset of the oxygenated world of early Earth.

Methods

Biogeochemical box model

We developed a biogeochemical box model that considers the biogeochemical cycles of C, P, S, Fe, O₂, and Ca in the ocean–atmosphere system, which is an improved version of the previous study⁴⁴. More specifically, the original model was improved by including the dynamical S cycles and C and S isotopes (δ¹³C and δ³⁴S). The new model considers the dynamical responses of atmospheric CO₂, CH₄, and O₂ concentrations, oceanic levels of DIC, carbonate alkalinity (Alk_c), [PO₄^3–], [SO₄^2–], [Fe²⁺], [Ca²⁺], dissolved O₂, and dissolved CH₄, and DIC, Alk_c, and [Ca²⁺] in the pore spaces of oceanic crust. Atmospheric budgets of O₂, CH₄, and CO₂ via photochemical reactions are calculated using a parameterization created in the previous study⁴⁷. The surface temperature is calculated using the temperature parameterization of Watanabe and Tajika⁶⁴, which accounts for the greenhouse effects of CO₂, H₂O, CH₄, and C₂H₆ based on the results of the radiative-convective climate models^100,101. The atmospheric CO₂ is removed via the weathering of the continental crust and the gas exchange with the ocean⁴⁴.

The ocean is divided into a low-latitude surface ocean box (<100 m depth and 90% of the oceanic area), a high-latitude surface ocean box (<400 m depth and 10% of the oceanic area), an intermediate ocean box (100–1500 m depth below the low-latitude surface ocean box), and a deep ocean box⁴⁴. The oceanic DIC is removed via the deposition of organic matter in the sediments and calcium carbonate in sediments and in the pore spaces of the oceanic crust. Marine phosphate is supplied via rivers and removed from the ocean via the deposition of organic matter and adsorption to Fe-bearing minerals and Ca-bound P burial. Marine Fe(II) is supplied via rivers and hydrothermal systems and removed from the ocean via the deposition of Fe(III) hydroxide and pyrite in sediments. The primary production in the ocean is driven by OP and Fe-(II)-based anoxygenic photoautotrophs (AP). O₂ is produced in surface oceans by OPs and is partly consumed by methanotrophs and the oxidation of Fe(II), with the remaining O₂ being outgassed into the atmosphere. CH₄ is generated by the decomposition of organic matter through methanogenesis and is partly consumed by methanotrophs, with the remainder outgassed into the atmosphere. In this model, the activities of H₂-based photoautotrophs and/or H₂-/CO-consuming chemoautotrophs are not considered. Their activities would increase the primary productivity of APs under low atmospheric pO₂ conditions^48,102,103. If their activities are considered, the atmospheric pCH₄ and the surface temperature under Archean-like low atmospheric pO₂ conditions would likely increase. However, this would not strongly alter the fundamental behavior of the atmosphere–biosphere interactions after the eruption of LIPs and whiffs.

Carbon isotopes

The budget of ¹³C of the atmospheric CO₂ and CH₄ are represented as follows:

$$mu frac{d}{{dt}}({p}^{13}C{O}_{2}) =, {f}_{13,{DICs}}{F}_{{oa},{co}2,s,uparrow }-{f}_{13,{co}2a}{F}_{{oa},{co}2,s,downarrow }+{f}_{13,{DICh}}{F}_{{oa},{co}2,h,uparrow }\ -{f}_{13,{co}2a}{F}_{{oa},{co}2,h,downarrow }+{f}_{13,{vc}}{F}_{{vc}}+{f}_{13,{wog}}{F}_{{wog}}+{f}_{13,{ch}4}{F}_{{oxi},{co}2}\ +{f}_{13,{ch}4}{F}_{{esc}}-{f}_{13,{co}2a}{F}_{{ws}}-{f}_{13,{co}2a}{F}_{{wc}}$$

(1)

$$mu frac{d}{{dt}}left({p}^{13}C{H}_{4}right) =, {f}_{13,{ch}4s}{F}_{{oa},{ch}4,s,uparrow }-{f}_{13,{ch}4a}{F}_{{oa},{ch}4,s,downarrow }+{f}_{13,{ch}4h}{F}_{{oa},{ch}4,h,uparrow }\ -{f}_{13,{ch}4a}{F}_{{oa},{ch}4,h,downarrow }-{f}_{13,{ch}4a}{F}_{{oxi},{ch}4}-{f}_{13,{ch}4a}{F}_{{esc}}$$

(2)

where μ is the number of moles in the present atmosphere (μ = 1.773 × 10²⁰ mol); f₁₃ represents the abundance of ¹³C relative to total C; f_13,vc is f₁₃ of the CO₂ outgassing from the interior of the Earth; f_13,wog is f₁₃ of the sedimentary organic matter; the subscripts ↑ and ↓ represent the upward and downward counterpart, respectively; F_oa,co2 is the gas exchange rate between the ocean and atmosphere; F_vc is the volcanic outgassing rate of CO₂; F_wog is the oxidative weathering rate of sedimentary organic matters; F_oxi is the oxidation rate of atmospheric CH₄ to form CO₂; F_esc is the consumption rate of atmospheric CH₄ via H₂ escape; and F_ws and F_wc are the terrestrial continental weathering rate of silicates and carbonates. The budgets of ¹³C of the oceanic and pore-space DIC are represented as follows:

$${V}_{s}frac{d}{{dt}}({DI}{{C}^{13}}_{s}) = – {f}_{13,{borg},s}{F}_{{po},s}-{f}_{13,{methano},s}{F}_{{dgch}4,s}-{f}_{13,{pc},s}{F}_{{pc},s}-{f}_{13,{DICs}}{F}_{{oa},{co}2,s,uparrow }\ +{f}_{13,{co}2a}{F}_{{oa},{co}2,s,downarrow }+{f}_{13,{co}2a},{f}_{{areas}}({F}_{{ws}}+{F}_{{wc}})+{f}_{13,{sedc}},{f}_{{areas}}{F}_{{wc}}\ +{f}_{13,{ch}4s}({F}_{{met},{co}2,s}+{F}_{{AOM},s})+({,,f}_{13,{DICm}}{F}_{{dif},{ms},uparrow }-{f}_{13,{DICs}}{F}_{{dif},{ms},downarrow })\ -({,,f}_{13,{DICs}}{F}_{{dif},{sh},to }-{f}_{13,{DICh}}{F}_{{dif},{sh},leftarrow })+({,,f}_{13,{DICm}}{F}_{{adv},{ms},uparrow }-{f}_{13,{DICs}}{F}_{{adv},{sh},to })$$

(3)

$${V}_{h}frac{d}{{dt}}({DI}{{C}^{13}}_{h}) = -{f}_{13,{borg},h}{F}_{{po},h}-{f}_{13,{methano},h}{F}_{{dgch}4,h}-{f}_{13,{pc},h}{F}_{{pc},h}\ -{f}_{13,{DICh}}{F}_{{oa},{co}2,h,uparrow }+{f}_{13,{co}2a}{F}_{{oa},{co}2,h,downarrow }+({,,f}_{13,{DICs}}{F}_{{dif},{sh},to }\ -{f}_{13,{DICh}}{F}_{{dif},{sh},leftarrow })-({,,f}_{13,{DICh}}{F}_{{dif},{hd},downarrow }-{f}_{13,{DICd}}{F}_{{dif},{hd},uparrow })\ +({,,f}_{13,{DICs}}{F}_{{dif},{sh},to }-{f}_{13,{DICh}}{F}_{{dif},{hd},downarrow })+{f}_{13,{co}2a},{f}_{{areah}}({F}_{{ws}}\ +{F}_{{wc}})+{f}_{13,{sedc}},{f}_{{areah}}{F}_{{wc}}+{f}_{13,{ch}4h}({F}_{{met},{co}2,h}+{F}_{{AOM},h})$$

(4)

$${V}_{m}frac{d}{{dt}}({DI}{{C}^{13}}_{m}) =, {f}_{13,{borg},s}({F}_{{po},s}-{F}_{{bo},m})-{f}_{13,{methano},m}{F}_{{dgch}4,m}\ -({,,f}_{13,{pc},m}{F}_{{pc},m}-{f}_{13,{pc},s}{F}_{{dc},m})+{f}_{13,{ch}4m}({F}_{{met},{co}2,m}+{F}_{{AOM},m})\ + ({,,f}_{13,{DICd}}{F}_{{dif},{dm},uparrow }-{f}_{13,{DICm}}{F}_{{dif},{dm},downarrow })-({,,f}_{13,{DICm}}{F}_{{dif},{ms},uparrow }\ -{f}_{13,{DICs}}{F}_{{dif},{ms},downarrow })+({,,f}_{13,{DICd}}{F}_{{adv},{dm},uparrow }-{f}_{13,{DICm}}{F}_{{adv},{ms},uparrow })$$

(5)

$${V}_{d}frac{d}{{dt}}({DI}{{C}^{13}}_{d}) =, {f}_{13,{borg},s}(1-{alpha }_{m}){F}_{{bo},m}+{f}_{13,{borg},h}(1-{{alpha }_{{sd}}}^{2}){F}_{{po},h}\ -{f}_{13,{methano},d}{F}_{{dgch}4,d}+{f}_{13,{ch}4d}({F}_{{met},{co}2,d}+{F}_{{AOM},d})+{f}_{13,{pc},s}{F}_{{dc},s,d}+{f}_{13,{pc},h}{F}_{{dc},h,d}\ +{f}_{13,{pc},m}{F}_{{dc},m,d}+{f}_{13,{bpyrsdc}}{F}_{{bpyr},{sd},c}+({,,f}_{13,{DICh}}{F}_{{dif},{hd},downarrow }-{f}_{13,{DICd}}{F}_{{dif},{hd},uparrow })\ -({,,f}_{13,{DICd}}{F}_{{dif},{dm},uparrow }-{f}_{13,{DICm}}{F}_{{dif},{dm},downarrow })+({,,f}_{13,{DICh}}{F}_{{adv},{hd},downarrow }-{f}_{13,{DICd}}{F}_{{adv},{dm},uparrow })\ +({,,f}_{13,{DICp}}{F}_{{adv},{dp},uparrow }-{f}_{13,{DICd}}{F}_{{adv},{dp},downarrow })$$

(6)

$${V}_{p}frac{d}{{dt}}({DI}{{C}^{13}}_{p})=({,,f}_{13,{DICd}}{F}_{{adv},{dp},downarrow }-{f}_{13,{DICp}}{F}_{{adv},{dp},uparrow })-{f}_{13,{pc},p}{F}_{{pc},p}$$

(7)

where the subscripts → and ← represent the fluxes from the low-latitude surface ocean box to the high-latitude surface ocean box and vice versa, respectively; f_area is the areal fraction of the surface ocean box in total oceanic area; V is the volume of each ocean box; f_13,borg and f_13,pc are f₁₃ of the organic matter and carbonate depositing in the ocean, respectively; f_13,sedc is f₁₃ of the terrestrial carbonates; f_13,methano is f₁₃ of CH₄ produced via fermentation and methanogenesis; f_13,bpyrsdc is f₁₃ of the OC decomposed during the formation of pyrite in sediments; α_m and α_sd is the transfer efficiency of organic matter in the intermediate water box and deep ocean box below the low-latitude surface ocean box, respectively; F_dif and F_adv are the water exchange rate between ocean boxes via diffusion and advection; F_met is the oxidation rate of CH₄ by methanotrophs; F_dgc is the decomposition rate of organic matter in each ocean box; F_pc and F_dc are the formation and dissolution rates of CaCO₃ in each ocean box; F_po is the export production rate; F_bo is the burial rate of organic matter from each ocean box; F_AOM is the abiotic oxidation rate of methane by marine sulfate; F_bpyr,sd,c is the rate of OC decomposition during the formation of pyrite in sediments. The budgets of ¹³C of the oceanic CH₄ are represented as follows:

$${V}_{s}frac{d}{{dt}}([{{}^{13}}{CH}_{4}]_{s}) = -({,,f}_{13,{ch}4s}{F}_{{oa},{ch}4,s,uparrow }-{f}_{13,{ch}4a}{F}_{{oa},{ch}4,s,downarrow })+{f}_{13,{methano},s}{F}_{{dgch}4,s}\ -{f}_{13,{ch}4s}{F}_{{met},{ch}4,s}-{f}_{13,{ch}4s}{F}_{{AOM},s}+({,,f}_{13,{ch}4m}{F}_{{dif},{ms},uparrow }-{f}_{13,{ch}4s}{F}_{{dif},{ms},downarrow })\ -({,,f}_{13,{ch}4s}{F}_{{dif},{sh},to }-{f}_{13,{ch}4h}{F}_{{dif},{sh},leftarrow })+({,,f}_{13,{ch}4m}{F}_{{adv},{ms},uparrow }-{f}_{13,{ch}4s}{F}_{{adv},{sh},to })$$

(8)

$${V}_{h}frac{d}{{dt}}([{{}^{13}}{CH}_{4}]_{h}) = -({,,f}_{13,{ch}4h}{F}_{{oa},{ch}4,h,uparrow }-{f}_{13,{ch}4a}{F}_{{oa},{ch}4,h,downarrow })+{f}_{13,{methano},h}{F}_{{dgch}4,h}\ -{f}_{13,{ch}4h}{F}_{{met},{ch}4,h}-{f}_{13,{ch}4h}{F}_{{AOM},h}+({,,f}_{13,{ch}4s}{F}_{{dif},{sh},to }-{f}_{13,{ch}4h}{F}_{{dif},{sh},leftarrow })\ -({,,f}_{13,{ch}4h}{F}_{{dif},{hd},downarrow }-{f}_{13,{ch}4d}{F}_{{dif},{hd},uparrow })+({,,f}_{13,{ch}4s}{F}_{{dif},{sh},to }-{f}_{13,{ch}4h}{F}_{{dif},{hd},downarrow })$$

(9)

$${V}_{m}frac{d}{{dt}}([{{}^{13}}{{CH}}_{4}]_{m}) =, {f}_{13,{methano},m}{F}_{{dgch}4,m}-{f}_{13,{ch}4m}{F}_{{met},{ch}4,m}-{f}_{13,{ch}4m}{F}_{{AOM},m}\ +({,,f}_{13,{ch}4d}{F}_{{dif},{dm},uparrow }-{f}_{13,{ch}4m}{F}_{{dif},{dm},downarrow })-({,,f}_{13,{ch}4m}{F}_{{dif},{ms},uparrow }-{f}_{13,{ch}4s}{F}_{{dif},{ms},downarrow })\ +({,,f}_{13,{ch}4d}{F}_{{adv},{dm},uparrow }-{f}_{13,{ch}4m}{F}_{{adv},{ms},uparrow })$$

(10)

$${V}_{d}frac{{d}}{{dt}}([{{}^{13}}{CH}_{4}]_{d}) =, {f}_{13,{methano},d}{F}_{{dgch}4,d}-{f}_{13,{ch}4d}{F}_{{met},{ch}4,d}-{f}_{13,{ch}4d}{F}_{{AOM},d}\ +({,,f}_{13,{ch}4h}{F}_{{dif},{hd},downarrow }-{f}_{13,{ch}4d}{F}_{{dif},{hd},uparrow })-({,,f}_{13,{ch}4d}{F}_{{dif},{dm},uparrow }\ -{f}_{13,{ch}4m}{F}_{{dif},{dm},downarrow })+({,,f}_{13,{ch}4h}{F}_{{adv},{hd},downarrow }-{f}_{13,{ch}4d}{F}_{{adv},{dm},uparrow })$$

(11)

where F_dgch4 is the production rate of CH₄ by methanogenesis.

Global sulfur cycle

We introduced the dynamical S cycle to the biogeochemical box model. The budgets of marine sulfate (SO₄^2–), marine hydrogen sulfide (ΣH₂S), crustal pyrite (S_pyr), and crustal gypsum (S_gyp) are represented as follows:

$${V}_{{oc}}frac{d}{{dt}}([S{{O}_{4}}^{2-}])={F}_{{dpyr}}+{F}_{{dgyp}}+{F}_{{wpyr}}+{F}_{{wgyp}}-{F}_{{dg},S}-{F}_{{AOM}}+{F}_{{oxhs}}-{F}_{{bgyp}}-{F}_{{bpyr},{sd}}$$

(12)

$${V}_{{oc}}frac{d}{{dt}}([varSigma {H}_{2}S])={F}_{{dg},S}+{F}_{{AOM}}-{F}_{{oxhs}}-{F}_{{bpyr},{wc}}$$

(13)

$$frac{d}{{dt}}({S}_{{gyp}})={F}_{{bgyp}}-{F}_{{dgyp}}-{F}_{{wgyp}}$$

(14)

$$frac{d}{{dt}}({S}_{{pyr}})={F}_{{bpyr}}-{F}_{{dpyr}}-{F}_{{wpyr}}$$

(15)

where F_dpyr and F_dgyp are the outgassing rates of S from volcanos via metamorphic decompositions of pyrite and gypsum, respectively; F_wpyr is the oxidative weathering rate of sedimentary pyrite; F_wgyp is the weathering rate of sedimentary gypsum; F_dg,S is the decomposition rate of organic C by marine sulfate; F_AOM is the rate of abiotic oxidation of methane (AOM) by marine sulfate; F_oxhs is the oxidation rate of ΣH₂S; F_bgyp and F_bpyr are the gypsum and pyrite burial fluxes from the ocean, respectively. The budgets of ³⁴S are represented as follows:

$${V}_{{oc}}frac{d}{{dt}}([{{}^{34}}S{{O}_{4}}^{2-}]) =, {f}_{34,{pyr}}{F}_{{dpyr}}+{f}_{34,{gyp}}{F}_{{dgyp}}+{f}_{34,{pyr}}{F}_{{wpyr}}\ +{f}_{34,{gyp}}{F}_{{wgyp}}-{f}_{34,{so}4}{F}_{{dg},S}-{f}_{34,{so}4}{F}_{{AOM}}\ +{f}_{34,{hs}}{F}_{{oxhs}}-{f}_{34,{so}4}{F}_{{bgyp}}-{f}_{34,{bpyrsd}}{F}_{{bpyr},{sd}}$$

(16)

$${V}_{{oc}}frac{d}{{dt}}([varSigma {{H}_{2}}{{}^{34}}S])={f}_{34,{so}4}{F}_{{dg},S}+{f}_{34,{so}4}{F}_{{AOM}}-{f}_{34,{hs}}{F}_{{oxhs}}-{f}_{34,{bpyrwc}}{F}_{{bpyr},{wc}}$$

(17)

$$frac{d}{{dt}}({{}^{34}}{S}_{{gyp}})={f}_{34,{so}4}{F}_{{bgyp}}-{f}_{34,{gyp}}{F}_{{dgyp}}-{f}_{34,{gyp}}{F}_{{wgyp}}$$

(18)

$$frac{d}{{dt}}({{}^{34}}{{S}_{{pyr}}})={f}_{34,{bpyr}}{F}_{{bpyr}} – {f}_{34,{pyr}}{F}_{{dpyr}} – {f}_{34,{pyr}}{F}_{{wpyr}}$$

(19)

where f₃₄ represents the abundance of ³⁴S relative to total S; and f_34,bpyr represents f₃₄ of pyrite formed in the ocean, respectively.

The burial of gypsum from the water column is assumed to occur in the low-latitude surface ocean box, as follows^104,105,106:

$${F}_{{bgyp}}=frac{[S{{O}_{4}}^{2-}]_{s}}{[S{{O}_{4}}^{2-}]_{0}}cdot frac{[C{a}^{2+}]_{s}}{[C{a}^{2+}]_{0}}cdot {F}_{{bgyp},0}$$

(20)

where [SO₄^2–]₀ and [Ca²⁺]₀ are the present oceanic concentration of sulfate and calcium ions, respectively, and k_bgyp,0 is the rate constant for the burial of sedimentary gypsum. The buried gypsum enters the sedimentary gypsum reservoir. The sedimentary gypsum experiences metamorphism volcanism and provides S to the ocean–atmosphere system:

$${F}_{{dgyp}}={r}_{{sp}}cdot {k}_{{dgyp}}cdot {S}_{{gyp}}$$

(21)

where r_sp is the seafloor spreading rate relative to present, and k_dgyp is the rate constant for the pyrite degassing. The weathering rate of sedimentary gypsum is represented as follows:

$${F}_{{wgyp}}={F}_{{wgyp},0}cdot frac{{F}_{{ws}}}{{F}_{{ws},0}}cdot frac{{S}_{{gyp}}}{{S}_{{gyp},0}}$$

(22)

where F_wgyp,0 and F_ws,0 are the present sedimentary gypsum and continental silicate weathering fluxes, and S_gyp,0 is the present sedimentary gypsum reservoir size.

The burial of pyrite is assumed to occur in the water column (F_bpyr,wc) and sediment (F_bpyr,sd):

$${F}_{{bpyr}}={F}_{{bpyr},{wc}}+{F}_{{bpyr},{sd}}$$

(23)

The stoichiometric reaction for the pyrite formation in the water column is represented as follows:

$$2{{{rm{H}}}}_{2}{{rm{S}}}+{{{rm{Fe}}}}^{2+}to {{rm{Fe}}}{{{rm{S}}}}_{2}+{{{rm{H}}}}_{2}+{2{{rm{H}}}}^{+}$$

(24)

The burial of pyrite from the water column (F_bpyr,wc) is estimated as follows¹⁰⁷:

$${F}_{{bpyr},{wc}}={k}_{{bpyr},{wc}}cdot {varSigma }_{i=s,h,m,d}{V}_{i}cdot [F{e}^{2+}]_{i}cdot [varSigma {H}_{2}S]$$

(25)

where k_bpyr,wc is the rate constant for the burial of pyrite from the water column. For pyrite burial in sediments, on the other hand, we assumed that part of H₂S originating from MSR in sediments deposits is pyrite. The stoichiometric reaction for pyrite formation in sediments is represented as follows:

$${8{{rm{S}}}{{{rm{O}}}}_{4}}^{2-}+{16{{rm{H}}}}^{+}+4{{rm{Fe}}}{({{rm{OH}}})}_{3}+17{{rm{C}}}{{{rm{H}}}}_{2}{{rm{O}}}to 4{{rm{Fe}}}{{{rm{S}}}}_{2}+4{{{rm{H}}}}_{2}+17{{{rm{CO}}}}_{2}+27{{{rm{H}}}}_{2}{{rm{O}}}$$

(26)

The burial of pyrite in sediments is represented as follows:

$${F}_{{bpyr},{sd}}=min left({k}_{{bpyr},{sd}}cdot frac{[S{O}_{4}]_{d}{2-atop}}{[S{O}_{4}]_{d}{2-atop}+{K}_{{MSR}}},frac{8}{17}{F}_{{bo},d}right)$$

(27)

where k_bpyr,sd is the rate constant for the burial of pyrite in sediments and K_MSR is the half-saturation constant for MSR. The net production rate of O₂ via deposition of pyrite in the water column and sediments via the production of H₂ can be represented in terms of the budget of O₂ as follows:

$${F}_{{bpyr},o2}({Tmol},,{O}_{2},,{y}{r}^{-1})=-frac{1}{4}{F}_{{bpyr},{wc}}({Tmol; S; y}{r}^{-1})-frac{1}{4}{F}_{{bpyr},{sd}}({Tmol; S; y}{r}^{-1})$$

(28)

where the first term of the right-hand side represents the net consumption of O₂ via H₂ produced by the production of pyrite in the water column via Eq. 24 and the second term represents the production of H₂ in sediments via Eq. 26.

The buried pyrite enters the sedimentary pyrite reservoir. The sedimentary pyrite experiences metamorphism volcanism and oxidative weathering, whose rates are given as follows, respectively¹⁰⁸:

$${F}_{{dpyr}}={r}_{{sp}}cdot {k}_{{dpyr}}cdot {S}_{{pyr}}$$

(29)

$${F}_{{wpyr}}={F}_{{wpyr},0}cdot frac{{F}_{{ws}}}{{F}_{{ws},0}}cdot {left(frac{p{O}_{2}}{p{O}_{2,0}}right)}^{!!0.5}cdot frac{{S}_{{pyr}}}{{S}_{{pyr},0}}$$

(30)

where k_dpyr is the rate constant for pyrite degassing, F_wpyr,0 is the present oxidative weathering flux of sedimentary pyrite, and S_pyr,0 is the present sedimentary pyrite reservoir size. The stoichiometric reaction for pyrite degassing and oxidative weathering of pyrite is represented as follows:

$$4{{rm{Fe}}}{{{rm{S}}}}_{2}+15{{{rm{O}}}}_{2}+14{{{rm{H}}}}_{2}{{rm{O}}}to {8{{rm{S}}}{{{rm{O}}}}_{4}}^{2-}+16{{{rm{H}}}}^{+}+4{{rm{Fe}}}{({{rm{OH}}})}_{3}$$

(31)

The net effect of the metamorphic volcanism and oxidative weathering of pyrite on the global redox budget can be represented in terms of the consumption of O₂, as follows:

$${F}_{{wdpyr},o2}=-frac{15}{8}({F}_{{wpyr}}+{F}_{{dpyr}})$$

(32)

LIP eruption scenario

For the standard scenario of volcanism after an eruption of LIPs, we assumed an eruption scenario of the Ontong Java Plateau that caused the deoxygenation of the ocean, the so-called ocean anoxic event 1a at ~120 Ma. For this standard scenario, we assumed a period of intense volcanisms of 400 kyr³⁶ and a total C influx of 10²¹ g CO₂ (~2.7 × 10⁵ GtC) which is an estimated maximum value for the eruption of the Ontong Java Plateau⁴⁵. The volcanic gases are supplied to the atmosphere at a constant rate over the period of the eruption of LIP, assuming the emplacement of LIP on continents. In reality, it is likely that the emplacement on the oceanic crust was more frequent than on continents considering the small mass of the continental crust. Nevertheless, this treatment would be sufficient for representing the response of the system at a timescale longer than the gas exchange between the atmosphere and the ocean. We assumed that the supply rate of reducing gas (represented by H₂) relative to C from LIP volcanism (f_red,lip) is constant throughout the LIP emplacement. Despite the uncertainty in f_red,lip^109,110, we set f_red,lip to 0.36, which is equal to the relative supply rate of C and H₂ in the background steady state of the model⁴⁴. However, it has been discussed that the different pressure conditions owing to the different locations of the LIP emplacement result in different values of f_red,lip^43,111,112. The impact of the choice of f_red,lip was investigated in Fig. S5 and S6. The supply of Fe(II) from the erupted silicate minerals on continents is not considered because we assumed a continental LIP eruption. Note that the supply of Fe(II) would affect the conditions for the occurrence of LIPs in the case of the marine LIP: however, it is beyond the scope of this study.

Using this model, we estimated the response of the global C–P–S–Fe–O₂–Ca biogeochemical cycles to the intense volcanism after the eruption of LIPs under the reducing conditions during the late Archean. We first obtain the Archean-like background steady-state using the set of the standard parameters listed in Table 1. Starting from this background state, we investigate the response of the system to the volcanism assuming the eruption of LIPs. We simulate the response with respect to different background hydrothermal Fe(II) supply rates, f_a, background reducing gas outgassing rate, reducing gas outgassing rate from LIPs, total C influx, and duration of LIP (Table S5).

Table 1 Standard parameter values

Full size table