Module 6: The Great Oxygenation Event

Around 2.4 billion years ago, Earth's atmosphere underwent its most dramatic chemical transition: the partial pressure of O₂ rose from ~10^-5 of modern levels to ~1% in roughly 100 Myr. This Great Oxygenation Event (GOE) was a planetary-scale biogeochemical revolution driven by cyanobacterial oxygenic photosynthesis. It locked vast amounts of iron into banded iron formations, triggered the Huronian global glaciation, and set the stage — a billion years later — for the second oxygenation and the rise of complex life. Cross-links: Ecological Biochemistry, Earth Sciences.

6.1 The Pre-Oxygen Archean World

Before oxygenic photosynthesis, Earth's biosphere was entirely anaerobic. Primary production relied on chemolithotrophy (H₂, H₂S, Fe²⁺ as electron donors) and anoxygenic photosynthesis using bacteriochlorophylls:

\[ 2\,\text{H}_2\text{S} + \text{CO}_2 \xrightarrow{h\nu} \text{CH}_2\text{O} + \text{H}_2\text{O} + 2\,\text{S} \]

\[ 4\,\text{Fe}^{2+} + \text{CO}_2 + 11\,\text{H}_2\text{O} \xrightarrow{h\nu} \text{CH}_2\text{O} + 4\,\text{Fe(OH)}_3 + 8\,\text{H}^+ \]

These substrates are locally abundant but globally scarce — their flux is set by volcanic outgassing and hydrothermal circulation. Primary productivity was limited to roughly \(\sim 0.1\%\) of modern levels (Canfield 2014 estimate).

6.2 PSII & the Water-Splitting Complex

The origin of oxygenic photosynthesis required linking two photosystems in series to achieve enough redox potential to oxidise water. The standard Z-scheme:

\[ \text{H}_2\text{O} \xrightarrow{\text{PSII}} 1/2\,\text{O}_2 + 2\,e^- + 2\,\text{H}^+ \]

The PSII reaction centre reaches an oxidising potential of +1.25 V, enough to drive water oxidation (\(E^\circ = +0.82\) V at pH 7) through a four-electron process on the Mn₄CaO₅ oxygen-evolving cluster (OEC). The S-state cycle of Kok:

\[ S_0 \xrightarrow{h\nu} S_1 \xrightarrow{h\nu} S_2 \xrightarrow{h\nu} S_3 \xrightarrow{h\nu} S_4 \to S_0 + \text{O}_2 \]

Four sequential photon absorptions advance the oxidation state of the Mn cluster; on the fifth the cluster releases O₂ and returns to S₀. Crystal structures at 1.9 Å (Umena 2011) show the cubane Mn₄CaO₅ geometry in striking detail. The evolutionary assembly of PSII remains a major open question: Cardona (2019) argues it may have evolved as early as 3.5 Gya, long before the GOE, with initial oxygen being consumed locally before atmospheric accumulation became significant.

6.3 Banded Iron Formations

Banded iron formations (BIFs) are alternating millimetre- to metre-scale layers of iron oxides (hematite, magnetite) and chert or silica. They represent the largest geological reservoir of iron on Earth — the Hamersley Basin in Australia contains ~20×10¹² tonnes of iron ore. BIF deposition peaked between 2.6 and 2.4 Gya.

The basic reaction is oxidation of dissolved ferrous iron by dissolved O₂:

\[ 4\,\text{Fe}^{2+} + \text{O}_2 + 10\,\text{H}_2\text{O} \longrightarrow 4\,\text{Fe(OH)}_3 + 8\,\text{H}^+ \]

Rust precipitates on the ocean floor. As long as ferrous iron flux from hydrothermal vents exceeded oxygen production, all O₂ was consumed in BIF deposition and none accumulated in the atmosphere. When the two fluxes balanced (around 2.4 Gya), free O₂began to escape to the atmosphere. The BIF record ends (for ~1 Gyr) as dissolved Fe²⁺was exhausted from the oceans. The brief revival of BIFs at 0.75-0.6 Gya is associated with the Neoproterozoic snowball Earth events, during which sea-ice cover isolated the oceans and allowed Fe²⁺ to reaccumulate.

Alternative mechanism: anoxygenic photoferrotrophy

Kappler et al. (2005) showed that anoxygenic photoferrotrophic bacteria can deposit BIFs directly without oxygen, using sunlight and Fe²⁺ to fix CO₂:

\[ 4\,\text{Fe}^{2+} + \text{CO}_2 + 11\,\text{H}_2\text{O} \xrightarrow{h\nu} \text{CH}_2\text{O} + 4\,\text{Fe(OH)}_3 + 8\,\text{H}^+ \]

If this mechanism was important, pre-GOE BIFs may not be direct evidence of oxygen at all, but of chemoautotrophic anoxygenic photosynthesis.

6.4 Oxygen Mass Balance

The atmospheric oxygen inventory \(N_\text{O}_2\) obeys:

\[ \frac{dN_{O_2}}{dt} = P_\text{photo} - S_\text{org} - S_\text{Fe} - S_\text{volc} \]

where \(P_\text{photo}\) is photosynthetic production and the three sinks are organic-carbon re-oxidation (respiration + burial), Fe²⁺ oxidation (BIFs), and reduction by volcanic gases (H₂, SO₂, H₂S).

At steady state, O₂ accumulates only if organic carbon is buried:

\[ \Delta N_{O_2} = B_\text{org} - S_\text{volc} - S_\text{Fe} \]

where \(B_\text{org}\) is the burial rate of organic matter. Without burial, respiration converts every O₂ back to CO₂. Current estimates for the GOE: \(B_\text{org} \approx 0.4\) times the modern rate, sustained for \(10^8\) yr, yielded the observed ~1% PAL atmospheric O₂.

6.5 Snowball Earth & Albedo Feedback

Two episodes of global glaciation have been linked to oxygenation: the Huronian (2.4-2.1 Gya, coinciding with the GOE) and the Neoproterozoic Cryogenian (0.72-0.635 Gya, including the Sturtian and Marinoan glaciations). The Huronian has been attributed to methane collapse: as O₂ rose, it destroyed atmospheric CH₄ (a stronger greenhouse gas than CO₂), removing a warming gas and triggering global cooling.

Budyko-Sellers energy balance

The positive ice-albedo feedback (ice reflects sunlight, cooling the Earth, causing more ice) can be captured in a simple 0-D energy balance:

\[ C\frac{dT}{dt} = \frac{S_\odot}{4}\bigl(1 - \alpha(T)\bigr) - \sigma T^4 + F_\text{GHG} \]

where \(\alpha(T)\) is the temperature-dependent albedo, increasing sharply below \(\sim 270\) K as sea ice forms. For certain values of the solar constant, the system is bistable: a warm state (ice-free equator) and a cold state (global snowball) coexist. Once Earth tips into the snowball state, escaping requires accumulated volcanic CO₂ (since the carbonate-silicate thermostat is frozen off) to push the greenhouse forcing past the hysteresis threshold.

6.6 The Neoproterozoic Second Oxygenation

For almost 1 Gyr after the GOE, atmospheric oxygen remained stuck at around 1-10% of modern levels. This “Boring Billion” (1.8-0.8 Gya) saw relatively little macroevolutionary novelty. Then, starting around 800 Mya and accelerating through the Ediacaran, a second oxygenation event raised atmospheric O₂ to near modern levels. Proposed triggers:

Weathering of orogenic belts (Transgondwana, Pan-African) increasing phosphorus flux and primary productivity
Emergence of eukaryotic algae with larger cells, faster sinking rates, and higher carbon burial efficiency
Nitrogen-fixation expansion
Feedbacks with Sturtian/Marinoan snowball Earth cycles — glaciations grind rock and release nutrients, fuelling biological bloom during deglaciation

Whatever the trigger, the biological consequences were enormous: sufficient atmospheric O₂ enabled energetically expensive large body plans, setting the stage for the Cambrian explosion 540 Mya (Module 7).

6.7 Biological Consequences of Oxygen

Rising O₂ was simultaneously a toxicity crisis for most existing life and a metabolic windfall for those that adapted. Two key innovations defined the post-GOE biosphere:

Aerobic respiration: couples glucose oxidation to O₂ reduction, yielding \(\sim 32\) ATP per glucose vs only 2 ATP per glucose from anaerobic glycolysis — a 16-fold improvement.
Antioxidant defences: superoxide dismutase, catalase, glutathione, thioredoxin, and ascorbate protect cells from reactive oxygen species (superoxide, H₂O₂, hydroxyl radical).

The combination of high ATP yield and atmospheric O₂ made energetically expensive multicellular body plans viable — a prerequisite for the Cambrian revolution.

Python: Atmospheric O₂ History & GOE Balance

Panel 1 traces atmospheric O₂ partial pressure on a log scale from 4 Gya to the present, highlighting the GOE, the Lomagundi peak, the Boring Billion plateau, and the Neoproterozoic second rise. Panel 2 zooms in on the GOE itself: the crossover at ~2.4 Gya where biological O₂ production overtook the declining Fe²⁺ reductant flux, allowing free oxygen to accumulate rather than be locked up in BIFs.

Python

script.py99 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

# Atmospheric O2 accumulation model
# dO2/dt = P_photo(t) - S_sink(t) - S_oxid(t)
# Production rises after cyanobacteria ~ 2.7 Gya
# Sinks: BIF iron oxidation (dominant before GOE), volcanic reductants, sediment burial

fig, axes = plt.subplots(1, 2, figsize=(16, 7), facecolor='#03120a')
fig.subplots_adjust(wspace=0.35)

# ---- Panel 1: Atmospheric O2 history ----
ax1 = axes[0]
ax1.set_facecolor('#061a0d')

t_Gya = np.linspace(4.0, 0.0, 1000)
dt_from_start = 4.0 - t_Gya

# Stepwise model: before 2.7 Gya ~ 1e-8 PAL
# GOE 2.4 Gya: jumps to 1% PAL
# Lomagundi 2.2 Gya: peak ~5-20% PAL
# Boring Billion 1.8-0.8 Gya: ~1-10% PAL
# Neoproterozoic oxygenation: second rise to ~50%
# Phanerozoic: ~1 PAL +/- 50%

def log_o2_model(t):
    # piecewise in Gya
    y = np.ones_like(t) * -7  # log10 PAL
    y[t < 2.7] = -4
    y[t < 2.4] = -2
    y[(t < 2.2) & (t > 2.0)] = -1.5
    y[(t <= 2.0) & (t > 0.8)] = -1.8 + 0.05 * np.sin(5 * (2.0 - t[(t <= 2.0) & (t > 0.8)]))
    y[(t <= 0.8) & (t > 0.54)] = -0.5 + 0.1 * (0.8 - t[(t <= 0.8) & (t > 0.54)])
    y[t <= 0.54] = np.log10(np.clip(1.0 - 0.3 * np.sin(8 * (0.54 - t[t <= 0.54])), 0.3, 1.5))
    return y

logO2 = log_o2_model(t_Gya)

ax1.plot(t_Gya, logO2, color='#34d399', lw=2.8)
ax1.fill_between(t_Gya, -7, logO2, alpha=0.25, color='#34d399')

# Events
ax1.axvline(2.4, color='#fbbf24', ls='--', lw=1.5, alpha=0.9)
ax1.text(2.45, -3, 'GOE\n2.4 Gya', color='#fbbf24', fontsize=10, fontweight='bold')

ax1.axvline(0.58, color='#f472b6', ls='--', lw=1.5, alpha=0.9)
ax1.text(0.63, -3.5, '2nd rise\n~ 0.6 Gya', color='#f472b6', fontsize=10, fontweight='bold')

ax1.axvspan(1.8, 0.8, alpha=0.1, color='#64748b', label='"Boring Billion"')

ax1.invert_xaxis()
ax1.set_xlabel('Time (Gya)', color='#6ee7b7', fontsize=11)
ax1.set_ylabel('log10 [pO2 / PAL]  (Present Atmospheric Level = 1)', color='#6ee7b7', fontsize=11)
ax1.set_title('Atmospheric O2 History: 4 Gya to Present',
              color='#a7f3d0', fontsize=13, fontweight='bold')
ax1.legend(fontsize=10, labelcolor='#a7f3d0', facecolor='#061a0d', edgecolor='#10b981')
ax1.tick_params(colors='#6ee7b7')
ax1.grid(True, alpha=0.2, color='#10b981')
ax1.set_ylim(-7.5, 0.8)
for sp in ax1.spines.values(): sp.set_edgecolor('#10b981')

# ---- Panel 2: O2 vs iron sink balance before and through the GOE ----
ax2 = axes[1]
ax2.set_facecolor('#061a0d')

t2 = np.linspace(2.7, 2.0, 400)
# Iron flux declining
Fe_flux = 1.5 * np.exp((t2 - 2.0) / 0.2)
# O2 production rising
O2_prod = 0.2 + 1.2 * (1 - np.exp(-(2.7 - t2) / 0.15))

ax2.plot(t2, Fe_flux, color='#f87171', lw=2.5, label='Fe2+ reductant flux (BIF deposition)')
ax2.plot(t2, O2_prod, color='#34d399', lw=2.5, label='Biological O2 production')
idx_cross = np.argmin(np.abs(Fe_flux - O2_prod))
ax2.axvline(t2[idx_cross], color='#fbbf24', ls=':', lw=2, label=f'Crossover ~ {t2[idx_cross]:.2f} Gya')
ax2.fill_between(t2, 0, np.minimum(Fe_flux, O2_prod), alpha=0.2, color='#f87171',
                 label='Locked up as BIFs')
ax2.fill_between(t2, np.minimum(Fe_flux, O2_prod), np.maximum(Fe_flux, O2_prod),
                 where=(O2_prod > Fe_flux), alpha=0.3, color='#34d399',
                 label='Free O2 accumulates')

ax2.invert_xaxis()
ax2.set_xlabel('Time (Gya)', color='#6ee7b7', fontsize=11)
ax2.set_ylabel('Relative flux', color='#6ee7b7', fontsize=11)
ax2.set_title('GOE Balance: O2 Production vs Iron Sink',
              color='#a7f3d0', fontsize=13, fontweight='bold')
ax2.legend(fontsize=8, labelcolor='#a7f3d0', facecolor='#061a0d', edgecolor='#10b981')
ax2.tick_params(colors='#6ee7b7')
ax2.grid(True, alpha=0.2, color='#10b981')
for sp in ax2.spines.values(): sp.set_edgecolor('#10b981')

fig.suptitle('The Great Oxygenation Event (2.4 Gya)',
             color='#a7f3d0', fontsize=14, fontweight='bold')
plt.savefig('output.png', bbox_inches='tight', facecolor=fig.get_facecolor(), dpi=120)
plt.close()
print('GOE simulation complete.')

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Python: Snowball Earth Bistability & BIF Timeline

Panel 1 solves the Budyko-Sellers energy-balance model with ice-albedo feedback, showing the bistable hysteresis region around the Neoproterozoic solar constant — conditions under which a small perturbation can tip Earth into (and out of) a global snowball state. Panel 2 plots the relative BIF deposition rate through geological time, showing the Hamersley/Transvaal peak around the GOE and the Rapitan revival during the Cryogenian.

Python

script.py93 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

# Snowball Earth: Budyko-Sellers energy balance model with ice-albedo feedback
# dT/dt = (S*(1 - alpha(T)) - sigma*T^4) / C

fig, axes = plt.subplots(1, 2, figsize=(16, 7), facecolor='#03120a')
fig.subplots_adjust(wspace=0.35)

sigma = 5.67e-8
C = 1e8   # heat capacity
def albedo(T):
    # Smooth step: low T -> high albedo (icy)
    return 0.35 + 0.3 * (1 - 1/(1 + np.exp((270 - T) / 3)))

# ---- Panel 1: Equilibrium T vs solar constant ----
ax1 = axes[0]
ax1.set_facecolor('#061a0d')

S_range = np.linspace(1000, 1600, 200)
# Hysteresis: find stable T by cooling sweep and warming sweep
T_cool = []
T_warm = []
T_start_cool = 300
T_start_warm = 210
for S in S_range[::-1]:
    T = T_start_cool
    for _ in range(15000):
        dT = (S * (1 - albedo(T)) - 4 * sigma * T**4) / (4 * 1e7)
        T = T + dT
    T_cool.append(T)
    T_start_cool = T
T_cool = T_cool[::-1]

for S in S_range:
    T = T_start_warm
    for _ in range(15000):
        dT = (S * (1 - albedo(T)) - 4 * sigma * T**4) / (4 * 1e7)
        T = T + dT
    T_warm.append(T)
    T_start_warm = T

ax1.plot(S_range, T_cool, color='#60a5fa', lw=2.5, label='Cooling branch (warm initial)')
ax1.plot(S_range, T_warm, color='#f97316', lw=2.5, ls='--', label='Warming branch (icy initial)')
ax1.axvline(1361, color='#fbbf24', ls=':', lw=1.5, label='Modern solar constant')
ax1.axvspan(1220, 1280, alpha=0.15, color='#93c5fd', label='Neoproterozoic S')

ax1.axhline(273, color='#fde047', lw=0.7, alpha=0.5)
ax1.text(1500, 275, 'freezing', color='#fde047', fontsize=9)

ax1.set_xlabel('Solar constant (W/m^2)', color='#6ee7b7', fontsize=11)
ax1.set_ylabel('Equilibrium surface T (K)', color='#6ee7b7', fontsize=11)
ax1.set_title('Budyko-Sellers: Ice-Albedo Bistability\nSnowball Earth hysteresis',
              color='#a7f3d0', fontsize=12, fontweight='bold')
ax1.legend(fontsize=9, labelcolor='#a7f3d0', facecolor='#061a0d', edgecolor='#10b981')
ax1.tick_params(colors='#6ee7b7')
ax1.grid(True, alpha=0.2, color='#10b981')
for sp in ax1.spines.values(): sp.set_edgecolor('#10b981')

# ---- Panel 2: BIF deposition rate over time ----
ax2 = axes[1]
ax2.set_facecolor('#061a0d')

t_Gya = np.linspace(4.0, 0.5, 400)
# Main BIF peak: 2.6-2.4 Gya (Hamersley, Transvaal, Superior)
# Rapaitan/Neoproterozoic revival: 0.75-0.6 Gya
peak1 = 12 * np.exp(-((t_Gya - 2.5)**2) / 0.07**2)
peak2 = 3.5 * np.exp(-((t_Gya - 0.65)**2) / 0.05**2)
trace = 1.5 * np.exp(-((t_Gya - 3.3)**2) / 0.15**2)
rate = peak1 + peak2 + trace + 0.05

ax2.plot(t_Gya, rate, color='#f87171', lw=2.6)
ax2.fill_between(t_Gya, 0, rate, alpha=0.3, color='#f87171')
ax2.axvline(2.4, color='#fbbf24', ls='--', lw=1.2, label='GOE')
ax2.axvline(0.72, color='#60a5fa', ls='--', lw=1.2, label='Sturtian snowball')
ax2.invert_xaxis()
ax2.set_xlabel('Time (Gya)', color='#6ee7b7', fontsize=11)
ax2.set_ylabel('Relative BIF deposition rate', color='#6ee7b7', fontsize=11)
ax2.set_title('Banded Iron Formations Through Time\n(Hamersley peak + Rapitan revival)',
              color='#a7f3d0', fontsize=12, fontweight='bold')
ax2.legend(fontsize=9, labelcolor='#a7f3d0', facecolor='#061a0d', edgecolor='#10b981')
ax2.tick_params(colors='#6ee7b7')
ax2.grid(True, alpha=0.2, color='#10b981')
for sp in ax2.spines.values(): sp.set_edgecolor('#10b981')

fig.suptitle('Snowball Earth & Banded Iron Formations',
             color='#a7f3d0', fontsize=14, fontweight='bold')
plt.savefig('output.png', bbox_inches='tight', facecolor=fig.get_facecolor(), dpi=120)
plt.close()
print('Snowball Earth and BIF simulation complete.')

Click Run to execute the Python code

Code will be executed with Python 3 on the server

References

Lyons, T.W., Reinhard, C.T. & Planavsky, N.J. (2014). The rise of oxygen in Earth's early ocean and atmosphere. Nature, 506, 307-315.
Holland, H.D. (2006). The oxygenation of the atmosphere and oceans. Philosophical Transactions B, 361, 903-915.
Umena, Y. et al. (2011). Crystal structure of oxygen-evolving photosystem II at 1.9 Å. Nature, 473, 55-60.
Kappler, A. et al. (2005). Deposition of banded iron formations by anoxygenic phototrophic Fe(II)-oxidizing bacteria. Geology, 33, 865-868.
Hoffman, P.F. & Schrag, D.P. (2002). The snowball Earth hypothesis. Terra Nova, 14, 129-155.
Budyko, M.I. (1969). The effect of solar radiation variations on the climate of the Earth. Tellus, 21, 611-619.
Sellers, W.D. (1969). A global climatic model based on the energy balance of the Earth-atmosphere system. Journal of Applied Meteorology, 8, 392-400.
Canfield, D.E. (2014). Oxygen: A Four Billion Year History. Princeton University Press.
Och, L.M. & Shields-Zhou, G.A. (2012). The Neoproterozoic oxygenation event. Earth-Science Reviews, 110, 26-57.
Cardona, T. (2019). Thinking twice about the evolution of photosynthesis. Open Biology, 9, 180246.
Bekker, A. et al. (2004). Dating the rise of atmospheric oxygen. Nature, 427, 117-120.
Kump, L.R. (2008). The rise of atmospheric oxygen. Nature, 451, 277-278.

←Previous: Panspermia Next: Major Evolutionary Transitions→