# Mercury isotope evidence for marine photic zone euxinia across the end-Permian mass extinction

Ruoyu Sun, Yi Liu, Jeroen E. Sonke, Zhang Feifei, Yaqiu Zhao, Yonggen Zhang, Jiubin Chen, Cong-Qiang Liu, Shuzhong Shen, Ariel D. Anbar, et al.

► To cite this version:

Ruoyu Sun, Yi Liu, Jeroen E. Sonke, Zhang Feifei, Yaqiu Zhao, et al. Mercury isotope evidence for marine photic zone euxinia across the end-Permian mass extinction. Communications Earth & Environment, 2023, 4 (159), 10.1038/s43247-023-00821-6. insu-04831614

HAL Id: insu-04831614

https://insu.hal.science/insu-04831614v1

Submitted on 12 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Mercury isotope evidence for marine photic zone euxinia across the end-Permian mass extinction

Ruoyu Sun<sup>1</sup>, Yi Liu<sup>1</sup>, Jeroen E. Sonke<sup>2</sup>, Zhang Feifei<sup>3</sup>, Yaqiu Zhao<sup>1</sup>, Yonggen Zhang<sup>1</sup>, Jiubin Chen<sup>1</sup>, Cong-Qiang Liu<sup>1</sup>, Shuzhong Shen<sup>3</sup>, Ariel D. Anbar<sup>4,5</sup> & Wang Zheng<sup>1✉</sup>

The driving forces, kill and recovery mechanisms for the end-Permian mass extinction (EPME), the largest Phanerozoic biological crisis, are under debate. Sedimentary records of mercury enrichment and mercury isotopes have suggested the impact of volcanism on the EPME, yet the causes of mercury enrichment and isotope variations remain controversial. Here, we model mercury isotope variations across the EPME to quantitatively assess the effects of volcanism, terrestrial erosion and photic zone euxinia (PZE, toxic, sulfide-rich conditions). Our numerical model shows that while large-scale volcanism remains the main driver of widespread mercury enrichment, the negative shifts of Δ¹⁹⁹Hg isotope signature across the EPME cannot be fully explained by volcanism or terrestrial erosion as proposed before, but require additional fractionation by marine mercury photoreduction under enhanced PZE conditions. Thus our model provides further evidence for widespread and prolonged PZE as a key kill mechanism for both the EPME and the impeded recovery afterward.

<footnote n="1">¹Institute of Surface-Earth System Science, School of Earth System Science, Tianjin University, Tianjin 300072, China. ²Laboratoire Géosciences Environnement Toulouse, CNRS/IRD/Université Paul Sabatier-Toulouse III, 31400 Toulouse, France. ³State Key Laboratory for Mineral Deposits Research, School of Earth Sciences and Engineering and Frontiers Science Center for Critical Earth Material Cycling, Nanjing University, 163 Xianlin Avenue, Nanjing 210023, China. ⁴School of Earth and Space Exploration, Arizona State University, Tempe, AZ 85287, USA. ⁵School of Molecular Sciences, Arizona State University, Tempe, AZ 85287, USA. ✉email: zhengw3@tju.edu.cn</footnote>
The end-Permian mass extinction (EPME, ~252 Ma) is the largest known Phanerozoic extinction, with a loss of ~81% of species in the ocean and ~89% of species on land¹. The causes, controls, and chronostratigraphic framework for the biotic extinction and recovery dynamics across the EPME have been extensively studied (see the recent review by Dal Corso et al.²). The main trigger of the EPME has been proposed to be the emplacement of large-scale volcanism, particularly the Siberian Traps large igneous province (STLIP)²,³. The long-term (~0.9 Myr) and multiphase LIP emplacement (2–7 × 10⁶ km³) released large amounts of volcanic gases including carbon dioxide, sulfur dioxide, mercury (Hg) and halogens, leading to abrupt climate and environmental changes that altered the hydrological cycle and the fluxes of bioessential elements (e.g., carbon, sulfur, nitrogen and phosphorous)⁴⁻⁶. This biogeochemical cascade eventually led to adverse marine conditions detrimental to organisms. However, the direct kill mechanism for marine biota during the EPME is still under debate.

Oceanic anoxia or euxinia ($H_2S$ rich conditions) stimulated by terrestrial nutrient influxes has been invoked as one of the key kill mechanisms for the EPME$^{7–10}$. Multiple lines of evidence point to enhanced terrestrial weathering and soil erosion and, consequently, increased inputs of organic matter (OM), sulfate, and nutrients to oceans$^{11,12}$. These processes may have increased

primary productivity and hence enhanced OM degradation via microbial sulfate reduction, producing widespread and prolonged anoxic and even euxinic conditions$^{7–9}$. A small number of studies reported evidence for photic zone euxinia (PZE), the occurrence of euxinic water in the photic zone, based on lipid biomarkers of anoxygenic photoautotrophs that require $H_2S$ to live (such as green and purple sulfur bacteria)$^{11,13,14}$. PZE typically occurs over a nutrient-rich continental shelf$^9$, which is the critical zone that supports marine primary productivity and hosts the majority of marine life$^{15}$. Therefore, the occurrence of PZE has been considered as a potent kill mechanism during almost all “The Big Five” Phanerozoic mass extinctions$^{16}$. However, the interpretations of biomarkers are sometimes ambiguous due to their susceptibility to contamination by anthropogenic or natural hydrocarbons and degradation during OM maturation$^{17–19}$. Furthermore, the mechanism that drove PZE as well as its relationship to volcanism, terrestrial erosion, and the EPME is not well understood. Therefore, the role of PZE during the EPME warrants further investigations.

Anomalous Hg enrichment in sedimentary strata has been widely used as key evidence linking large-scale volcanism (e.g., LIPs) and mass extinctions because volcanism is the dominant natural Hg source to the Earth's surface and Hg is dispersed globally via the atmosphere$^{20,21}$. In addition to Hg concentration, Hg stable isotope composition, especially the mass-independent fractionation of odd Hg isotopes (odd-MIF, represented by $\Delta^{199}\text{Hg}$)$^{22,23}$, is emerging as a useful tracer of Hg sources (volcanic vs. non-volcanic) and pathways (atmospheric deposition vs. terrestrial erosion) in paleoenvironments$^{24–26}$. Many marine sedimentary strata across the Permo–Triassic boundary (PTB) retrieved from different depositional settings show a short-lived (~several kyr) but pronounced spike (enrichment factor >10 in some strata) of Hg concentrations and Hg/TOC (Hg concentration normalized to a common host phase, total organic carbon) immediately before or at the onset of the main EPME interval, and a long-term (~400 kyr) but lower extent of Hg enrichment (enrichment factors = 2–8) afterward (see compilations in Supplementary Fig. 1, and by Shen et al.$^{24}$, Grasby et al.$^{20}$, and Dal Corso et al.$^{12}$). The short-lived spike of Hg and Hg/TOC was coincident with a sharp $\delta^{13}\text{C}$ minimum$^{27,28}$. Many sections also show a negative shift of $\Delta^{199}\text{Hg}$ by ~0.1‰ at or post-EPME relative to the Late Permian background levels before the extinction$^{24–26,29–33}$ (Supplementary Figs. 1–4; Supplementary Table 1). Notably, the most well-studied shallow marine (<100 m) Meishan section shows a short-lived, negative $\Delta^{199}\text{Hg}$ spike (~0.2‰ lower than the background values) within the main EPME interval, and this negative shift of $\Delta^{199}\text{Hg}$ postdated the Hg/TOC maximum and $\delta^{13}\text{C}$ minimum by ~40 kyr (Fig. 1a–c). Other sections (e.g., Daxiakou, Shangsi, Chaohu and Guryul River) with limited Hg isotope data also appear to show Hg enrichment prior to the negative shift of $\Delta^{199}\text{Hg}$ (Supplementary Fig. 4).

Previous studies often assumed the enhanced Hg enrichment and the negative shift of $\Delta^{199}\text{Hg}$ were synchronous across the EPME and interpreted them as the result of increased atmospheric Hg deposition from volcanic emission, and/or of massive terrestrial Hg input to oceans from biomass burning and soil erosion, which are often also linked to volcanism$^{24–26,29–33}$. Such interpretations follow because mantle-sourced volcanic Hg typically has circum-zero $\Delta^{199}\text{Hg}$ values$^{34,35}$, and modern terrestrial vegetation and soil typically have negative $\Delta^{199}\text{Hg}$ values$^{23}$, both of which are lower than background values of pre-EPME marine sediments (0.1–0.2‰) (Supplementary Table 1). However, these interpretations lack a full consideration on the impact of ocean redox conditions on Hg enrichment and Hg isotope variations, and they cannot account for the timing lag between the peaks of

<figure><img src="image_3.png" /><figcaption>**Fig. 1** Published geochemical data across the EPME. $\delta^{13}\text{C}$ of carbonates (a), Hg/TOC (b) and $\Delta^{199}\text{Hg}$ (c) of whole rocks, and reconstructed seawater temperature (d, based on $\delta^{18}\text{O}$ of conodont apatite) and seawater $\delta^7\text{Li}$ (e, based on $\delta^7\text{Li}$ of sedimentary rocks) from the Meishan section. f $\delta^{238}\text{U}$ of carbonate from global oceans (including South China, reflecting a global picture given the long residence time of seawater U). All the data in (a–f) have been chronologically calibrated. The vertical blue shading and dashed line: the main EPME interval and PTB; the range of grey dashed lines: approximate time range of STLIP sill intrusion phase. Note two samples with high Hg/TOC values (1825 and 3149 ppb/‰) at ~251.941 Ma (onset of the EPME) are excluded due to low TOC values (<0.2%). The 2σ analytic uncertainty (error bar) of $\Delta^{199}\text{Hg}$ is taken from the cited references (see data availability). Refer to Supplementary Fig. 1 for a global dataset including the Meishan section.</figcaption></figure>
$\Delta^{199}$Hg negative shift and Hg enrichment observed in the Meishan section. Thus, the mechanism of Hg enrichment and Hg isotopic variations across the EPME requires further studies.

To further elucidate the relationships between volcanic activity, enhanced weathering, and ocean anoxia during the EPME, here we develop a global box model of Hg isotopes (GBM-Hg-ISO-EPME) to investigate the mechanism of Hg enrichment and Hg isotopic variations across the EPME. This numerical model is adapted from previous models and comprises 13 fully-coupled Hg reservoirs (Supplementary Fig. 5, see Methods). It adopts a similar framework for the timing of Hg enrichment across the EPME as the Hg-C cycling model of Dal Corso et al.¹². It includes an elevated Hg enrichment phase of ~400 kyr (251.96–251.56 Ma) during the STLIP sill intrusion stage⁴,²⁸, and a peak Hg enrichment phase of ~1 kyr (251.950–251.949 Ma) immediately before the marine EPME starting at ~251.941 Ma. We mainly evaluate the model output against the published data of the Meishan section²⁴⁻²⁶,³¹,³² where the EPME is best studied and chronologically calibrated²⁸ (Supplementary Fig. 2). Other P-T sections deposited in global diverse marine settings were also evaluated by grouping their Hg isotope data in different time bins: pre-EPME, EPME, and post-EPME (Supplementary Figs. 3 and 4). We only focus on the $\Delta^{199}$Hg signature, because $\delta^{202}$Hg (representing mass-dependent fractionation) can be overprinted by multiple complex processes²³, and $\Delta^{200}$Hg (representing MIF of even Hg isotopes) variability is unresolvable relative to its analytical uncertainty in most EPME strata (Supplementary Fig. 2).

Our modeling results show that while large-scale volcanism remains a plausible driver of widespread Hg anomalies in sedimentary rocks across the EPME, perturbations due to volcanic (Run #1) or terrestrial Hg inputs (Run #2) alone cannot fully account for the extent of the negative shift of $\Delta^{199}$Hg. Instead, MIF in response to the development of marine PZE across the EPME can best explain both the extent of $\Delta^{199}$Hg negative shift (Run #3) and the timing difference between the peaks of Hg enrichment and $\Delta^{199}$Hg shift. Therefore, Hg isotopes provide new evidence for PZE as a kill mechanism for the EPME.

## Results and discussion

### Effect of subaerial volcanism on Hg enrichment and $\Delta^{199}$Hg.
Model Run #1 simulates how changes in Hg fluxes from subaerial volcanism affect Hg enrichment and $\Delta^{199}$Hg of various Earth reservoirs. To reproduce the long-term (~400 kyr) Hg enrichment (enrichment factors = 2–8) and the short-lived (~1 kyr) Hg spikes (enrichment factor >10) in global marine sediments across the EPME, the model requires an increase in subaerial volcanic Hg flux of 2–21 times and >21–27 (coastal-deep oceans) times its background value (~300 Mg a⁻¹), respectively (Supplementary Fig. 6a). The estimated increases of subaerial volcanic Hg flux are well within the previously estimated range (i.e., ~1,000–14,000 Mg a⁻¹, which is translated from the estimated total Hg emissions of 0.4 × 10⁹ to 5.6 × 10⁹ Mg for the ~400 kyr duration of STLIP sill intrusion stage, and is 2–45 times higher than background values)⁶,²⁶,³⁶. As shown in Fig. 2a, a long-term increase in subaerial Hg fluxes of 2–15 times with a short, intensive increase of ~50 times could explain the Hg enrichment pattern around the PTB in the Meishan section. However, Run #1 shows that the increase of subaerial volcanic Hg fluxes would increase $\Delta^{199}$Hg in sediments relative to the background values (Supplementary Fig. 6b), which is contrary to the negative shift of $\Delta^{199}$Hg observed across the EPME (Supplementary Fig. 3; Supplementary Table 1). This is because Hg emitted from subaerial volcanism would undergo photochemical transformations during atmospheric transport, which typically result in net positive $\Delta^{199}$Hg for Hg species that deposit to the ocean²²,²³ (see more discussion in

Supplementary Note 1). The model shows that an increase of submarine Hg fluxes could reproduce the observed marine Hg enrichment as well. However, because submarine Hg is poorly constrained during the EPME and its influence is thought to be limited to the vicinity of its source³⁷, we only focus on subaerial volcanic Hg in the model runs.

We additionally assess the potential contribution of Hg liberated from sediment degassing during sill intrusion to the long-term negative shift of $\Delta^{199}$Hg. As suggested by previous studies, the STLIP sill intrusion stage likely released massive amounts of Hg, CH₄, and CO₂ by heating organic-rich sediments in the Siberian Tunguska Basin⁵,³⁸,³⁹. The terrestrial organic-rich sediments (and sedimentary rocks like coals) in the Paleozoic⁴⁰ generally have more negative $\Delta^{199}$Hg values than mantle-sourced Hg (~0‰)³⁵. Thus, Hg emitted from sediment degassing rather than direct magmatic degassing could be a potential cause of the negative shift of $\Delta^{199}$Hg. However, assuming a negative $\Delta^{199}$Hg (~0.10‰ to ~0.30‰)⁴¹,⁴² for Hg from sediment degassing (Supplementary Fig. 7), our simulations show that the amount of degassing-sourced Hg required to reproduce the observed Hg enrichment with a ~0.1‰ shift of $\Delta^{199}$Hg is at least on the order of 10⁸ Mg, whereas the estimated sediment Hg release from STLIP sill intrusion is only ~2 × 10⁵ Mg⁶. Thus, Hg liberated from sediment degassing during sill intrusion cannot account for the long-term negative shift of $\Delta^{199}$Hg.

We further explore the potential effect of coal intrusion⁴³⁻⁴⁶, as coal deposits typically have negative $\Delta^{199}$Hg values (~0.11 ± 0.18‰, 1σ)⁴⁰. The initial emplacement pulse of intrusive sills which are hypothesized to drive the largest $\delta^{13}$C excursion and the onset of the EPME (Fig. 1) has been proposed to intrude into the West Siberian Coal Basin³⁸. It is likely that the coal deposits were rapidly metamorphosed and/or combusted, and massive amounts of Hg with negative $\Delta^{199}$Hg and isotopically light C were emitted into the atmosphere. Even assuming a lower-end $\Delta^{199}$Hg value of ~0.2‰ for the Siberian coal deposits⁴⁰, the model shows that > 6–7 × 10⁶ Mg of coal-derived Hg must be emitted in a duration of ~1 kyr to reproduce the short-lived Hg enrichment (enrichment factor >10) and $\Delta^{199}$Hg shift (ca. ~0.2‰ relative to the background value). However, a rough estimation assuming a mean coal Hg concentration of 100 ng g⁻¹⁴⁷ shows that > 6–7 × 10¹³ Mg of coal would have to have been combusted during the 1 kyr interval. This is more than one order of magnitude larger than the previously estimated values for coal reserves during that time (~0.3 × 10¹³ Mg)⁵,³⁹, and several times higher than the world’s total coal reserve at present (~1 × 10¹³ Mg). Thus, it seems unlikely that coal combustion caused the spikes of Hg enrichment and $\Delta^{199}$Hg.

### Effect of soil erosion on Hg enrichment and $\Delta^{199}$Hg.
In model Run #2, we test a second scenario: the input of terrestrial (soil) Hg to the oceans, which is generally characterized by negative $\Delta^{199}$Hg (~0.3‰ to ~0.1‰) in modern environments²²,²³. Across the EPME and the subsequent biotic recovery, soil erosion has been suggested to have increased greatly due to terrestrial ecosystem destruction resulting from the volcanism⁴⁸,⁴⁹. We test if a long-term increase in soil erosion rates could drive obvious Hg enrichment and $\Delta^{199}$Hg shift. Sensitivity analysis shows that an increase in soil erosion rates by up to 40 times could increase the marine Hg enrichment by 1–3 times and shift marine $\Delta^{199}$Hg values down to ~0.13‰ (Fig. 3). An increase in subaerial volcanic Hg flux is still needed to reproduce the observed marine Hg enrichment (enrichment factors = 2–8) around the PTB. To reproduce the long-term marine $\Delta^{199}$Hg shift, the model needs an increase of soil erosion rates by 3–5 times for the sedimentary sections (Meishan and Guryul Ravine) at the continental shelf in
<figure><img src="image_4.png" /></figure>

**Fig. 2** Simulated Hg enrichment and Δ¹⁹⁹Hg shift in model runs #1 and #3. Simulated Hg enrichment factor (EF) and Δ¹⁹⁹Hg shift in the coastal sediments in response to an increase in volcanism only (a), combined volcanism and soil erosion (b), and volcanism and coastal PZE (c). The amounts of soil erosion are expressed as times background soil erosion rate coefficients ($k_{se\_bg}$); the intensities of PZE are expressed as increasing isotope enrichment factors of odd-MIF associated with photoreduction of Hg(II) in the coastal ocean ($E_{occ(II-0),pho}$). The thick and thin solid lines represent simulations under high and low volcanic Hg emissions, respectively; the thick and thin solid lines of Δ¹⁹⁹Hg shift are largely overlapped due to insensitive response of Δ¹⁹⁹Hg to volcanic Hg emissions. The vertical shading and dashed lines are the same as Fig. 1. The simulated EF and Δ¹⁹⁹Hg shifts (lines) are calculated relative to their background steady-state values of the model. The measured EF and Δ¹⁹⁹Hg shifts (data symbols) are calculated relative to their background Hg/TOC (40 ppb/%) and Δ¹⁹⁹Hg values (0.09‰) of the Late Permian. Note two samples with high EF values (45 and 78) at -251.941 Ma (onset of the EPME) are excluded due to low TOC values (<0.2%). $F_{subaerial \ Hg\_bg}$: background subaerial volcanic Hg flux.

coastal oceans but an increase by >10 times for most sedimentary sections in open oceans (Fig. 3). Although soil erosion was enhanced at the PTB, long-term enhanced soil erosion rates that are >10 times higher than natural background rates seem unlikely.

Furthermore, the terrestrial crisis has been suggested to start earlier than the marine extinction², and the terrestrial ecosystem collapse has been proposed to result in the short-lived spikes of Hg enrichment just before the EPME in the Meishan section¹². Thus, we further explore if a fast increase of soil erosion rates within a short duration (~1 kyr) could drive these changes. Assuming a fast increase of soil erosion rates by 100 times, we find that soil erosion alone can only cause a transient and small increase of the marine Hg enrichment, followed by a large drop, and the negative shift of Δ¹⁹⁹Hg can at most reach -0.15‰ (Supplementary Fig. 8). These disproportionate changes in Hg enrichment and Δ¹⁹⁹Hg relative to the increase of soil erosion rate are because the soil Hg pool would be rapidly mobilized at high soil erosion rates, preventing continued terrestrial

contribution to marine Hg enrichment and Δ¹⁹⁹Hg shift. It is also important to note that the above simulation is actually unrealistic because it only considers the Hg input via soil erosion under background volcanic emission, but additional atmospheric Hg input via enhanced volcanism was also present, which should drive Δ¹⁹⁹Hg to more positive values. As shown in Fig. 2b, while an increase of soil erosion accompanied by less prominent volcanic Hg fluxes can explain the Hg enrichment pattern and the long-term Δ¹⁹⁹Hg shift around the PTB in the Meishan section, neither soil erosion alone nor soil erosion superimposed on enhanced volcanic Hg emission can reproduce the short-lived strong negative Δ¹⁹⁹Hg shift (~ -0.20‰).

Our conclusion based on the fully-coupled Hg isotope model is different from the previous Hg-C isotope model developed by Dal Corso et al.¹² who suggested that a millennial-scale increase of soil erosion by 100 times (ca. $4 \times 10^4$ Mg Hg a⁻¹ river flux, cumulative $4 \times 10^7$ Mg Hg) could generate the short-lived spikes of Hg enrichment and nadirs of δ²⁰²Hg values in the Meishan section. This was possible because Dal Corso et al.'s model¹²
<figure><img src="image_5.png" /><figcaption>Fig. 3 Simulated long-term Hg enrichment and Δ¹⁹⁹Hg shift in model run #2. Simulated enrichment factor (EF) and Δ¹⁹⁹Hg shift in coastal (a), surface (b), intermediate (c), and deep (d) marine sediments in response to increased soil erosion. The amounts of soil erosion are expressed as 1 (black line), 2 (orange line), 5 (sky blue line), 10 (bluish green line), 20 (yellow line) and 40 (blue line) times background soil erosion rate coefficients ($k_{se\_bg}$) under varying intensities of additional subaerial Hg fluxes. The simulated EF and Δ¹⁹⁹Hg shifts are calculated relative to their background steady-state values of the model. The measured Δ¹⁹⁹Hg shifts (data symbols) are calculated for the combined EPME and post-EPME relative to their pre-EPME Δ¹⁹⁹Hg values (Supplementary Table 1). For visual display, the Δ¹⁹⁹Hg shifts of marine sections (mean ± 1 SD) are anchored to arbitrary x-axis values given the insensitive responses of Δ¹⁹⁹Hg shifts to changes in volcanic Hg emission flux. $F_{subaerial \ Hg\_bg}$: background subaerial volcanic Hg flux.</figcaption></figure>

assumed invariant soil erosion fluxes that do not scale with the
changes of soil Hg reservoir size through time.

PZE drove the negative shift of marine Δ¹⁹⁹Hg. As shown in model Runs #1 and #2, neither volcanic emission nor soil erosion can fully account for the large negative shifts of Δ¹⁹⁹Hg across the EPME. Thus, in model Run #3, we add a third scenario: Hg isotope fractionation in the ocean caused by PZE. A pilot study on Mesoproterozoic black shale showed that sediments deposited under PZE tend to develop significantly more negative Δ¹⁹⁹Hg values (down to -0.20‰) than those deposited under oxic water columns⁵⁰. Two mechanisms have been proposed to explain the effect of PZE on sedimentary Δ¹⁹⁹Hg values: (1) photoreduction of Hg(II) in a sulfide-rich photic zone, and (2) enhanced oxidation of atmospheric Hg(0) in sulfidic surface water. Although the exact mechanism is still under investigation, experimental evidence⁵¹⁻⁵⁵ suggests that both mechanisms would shift sedimentary Δ¹⁹⁹Hg to more negative values. Multiple experiments have demonstrated that the sign of Δ¹⁹⁹Hg during aqueous photoreduction of Hg(II) depends on the type of ligands bonding with Hg. Photoreduction of oxygen-bound Hg(II) produces positive Δ¹⁹⁹Hg in the residual Hg(II), whereas photoreduction of sulfur-bound Hg(II) tends to produce negative Δ¹⁹⁹Hg in the residual Hg(II)⁵¹⁻⁵³ (Supplementary Fig. 9). Furthermore, sulfidic water is found to promote aqueous Hg(0) oxidation⁵⁴,⁵⁵. Thus, PZE may promote the uptake of atmospheric Hg(0), which typically has negative Δ¹⁹⁹Hg based on modern samples²³, via enhanced Hg(0) oxidation in the surface ocean and its sequestration to sediments. This oxidation process via a non-photochemical abiotic pathway has been shown to produce negative Δ¹⁹⁹Hg in the oxidized Hg(II) pool⁵⁶, which is also consistent with the negative shift of Δ¹⁹⁹Hg across the EPME. However, since the rate coefficients of dark Hg(0) oxidation are orders of magnitude smaller than those of Hg(II) photoreduction (Supplementary Table 2), the contribution of photoreduction to

global Hg flux is estimated to be much higher than that of dark
oxidation (Supplementary Fig. 5). Thus we focus on the odd-MIF
associated with photoreduction in the following model runs.

The extents of PZE and its impact on marine sediment Hg isotopic variations are likely different in different depositional settings. For example, coastal oceans (e.g., the Meishan section) typically develop stronger or more persistent PZE than the open ocean (e.g., the Shangsi and Chaohu sections) because of the abundant nutrient influx to coastal environments⁷,¹¹. Thus, in Run #3, we consider two scenarios: (1) PZE only occurred in coastal oceans (Supplementary Fig. 10); (2) PZE occurred in both coastal and open surface oceans (Fig. 4). In scenario 1, the extent of the Δ¹⁹⁹Hg shift in coastal oceans is highly sensitive to the isotope enrichment factor of odd-MIF associated with photoreduction of Hg(II) in the coastal ocean (E¹⁹⁹Hg$_{occ(II-0),pho}$). When E¹⁹⁹Hg$_{occ(II-0),pho}$ changes from -1.0‰ (representing photoreduction of oxygen-bound Hg(II)) to 1.5‰ (representing photoreduction of sulfur-bound Hg(II)), the coastal ocean Δ¹⁹⁹Hg values could shift significantly by -0.20‰ (Supplementary Fig. 10). Values of 0.25‰ and 1.50‰ for E¹⁹⁹Hg$_{occ(II-0),pho}$ could well reproduce the long-term and short-lived Δ¹⁹⁹Hg negative shifts, respectively, observed in the Meishan section (Fig. 2c). These fitted E¹⁹⁹Hg$_{occ(II-0),pho}$ values are compatible with experimental values determined during photoreduction of Hg(II) binding with reduced sulfur ligands (around 1.02‰)⁵¹ and photomicrobial reduction of Hg(II) in marine microalga (0.08–1.03‰) which are thought to contain abundant thiols⁵³. Furthermore, preliminary results from our new experiments investigating photoreduction of Hg(II) bound to inorganic S²⁻ also show similarly negative MIF in the remaining Hg(II) (Supplementary Note 2 and Fig. 9). It is reasonable to assume different E¹⁹⁹Hg$_{occ(II-0),pho}$ for the long-term and short-lived Δ¹⁹⁹Hg negative shifts, because the model-fitted enrichment factors represent the combined (net) effect of MIF generated by photoreduction of Hg(II) bound to different ligands, which is a
<figure><img src="image_6.png" /><figcaption>**Fig. 4** Simulated marine Δ¹⁹⁹Hg shift in response to isotope fractionation during open surface ocean Hg(II) photoreduction. Δ¹⁹⁹Hg shifts (representing the relative change of PZE intensity) in coastal (a) and open marine sediments of different depths (b surface; c intermediate; d deep) are calculated relative to their background steady-state values of the model. The isotope enrichment factors of odd-MIF associated with photoreduction of Hg(II) in the open surface ocean (E<sub>ocs(II-0),pho</sub>) are assumed to vary from −0.10‰ (black line) to 0.30‰ (yellow line), while the isotope enrichment factor of odd-MIF associated with photoreduction of Hg(II) in the coastal ocean (E<sub>occ(II-0),pho</sub>) is set at 0.50‰. The Δ¹⁹⁹Hg shifts (black squares) of marine sections deposited in different environmental settings are calculated as the changes of Δ¹⁹⁹Hg values for the EPME and post-EPME relative to the pre-EPME (Supplementary Table 1). For visual display, the Δ¹⁹⁹Hg shifts of marine sections (mean ± 1 SD) are anchored to arbitrary x-axis values given the insensitive responses of Δ¹⁹⁹Hg shifts to changes in volcanic Hg emission flux. F<sub>subaerial Hg_bg</sub>: background subaerial volcanic Hg flux.</figcaption></figure>

more realistic scenario in both modern natural aqueous environments<sup>52,57,58</sup> and under PZE conditions across the EPME. PZE was likely most expansive or persistent during the main extinction interval of the Meishan section<sup>13</sup> when the short-lived Δ¹⁹⁹Hg negative shift occurred. More expansive and persistent PZE would likely favor the photoreduction of sulfur-bound Hg(II) and hence the manifestation of more negative Δ¹⁹⁹Hg in the residual marine Hg(II)<sup>50</sup>. In contrast to the coastal ocean, all open ocean settings exhibit little Δ¹⁹⁹Hg changes (<0.01‰) with varying E<sup>199</sup>Hg<sub>occ(II-0),pho</sub> according to our model (Supplementary Fig. 10), which is inconsistent with the observed negative shifts of Δ¹⁹⁹Hg in open marine sediment sections. This is because the Hg exchange between coastal and open oceans are rather limited and the simulated different Hg species in coastal seawater (i.e., Hg(0), Hg(II), and Hg(P)) have opposite Δ¹⁹⁹Hg values, which cancel each other out when these Hg species are transported to open oceans. Therefore, while this scenario (PZE only in coastal ocean) can explain the negative shifts of Δ¹⁹⁹Hg in continental shelf settings (e.g., Meishan), it cannot satisfactorily explain the shifts of Δ¹⁹⁹Hg in distal slope and abyssal plain settings (Supplementary Fig. 10).

The second scenario that assumes PZE in both coastal and open oceans can satisfactorily explain the negative shifts of Δ¹⁹⁹Hg in almost all sections (Fig. 4). This assumption is reasonable as biomarker evidence for PZE has been reported for both continental shelf (e.g., Meishan)<sup>10,13</sup> and slope settings (e.g., Shangshi and Chaohu)<sup>11</sup>. PZE in the surface open ocean could develop via an oceanward expansion of euxinic water from the coastal ocean (e.g., continental shelf)<sup>11</sup> or via shoaling of euxinic water from the deep ocean<sup>59</sup>. Thus, it is plausible that

photoreduction of Hg(II) in the open ocean could exhibit similar PZE-driven MIF as the coastal ocean. A slight increase in E¹⁹⁹Hg<sub>ocs(II-0),pho</sub> from −0.1‰ to 0.2‰ in the open ocean could well explain the observed negative Δ¹⁹⁹Hg shifts in sections from intermediate and deep open ocean settings (Fig. 4). If we assume the E¹⁹⁹Hg<sub>ocs(II-0),pho</sub> is proportional to the extent of PZE as explained in scenario 1, then the above modeling results suggest that even a relatively small increase in the extent of PZE could explain the MIF data of the open ocean. Note that the modeled Δ¹⁹⁹Hg shifts in the coastal ocean in scenario 2 are different from those in scenario 1 with the increase of E¹⁹⁹Hg<sub>ocs(II-0),pho</sub> values. The reason is that PZE in the open ocean would emit Hg(0) with positive Δ¹⁹⁹Hg to the atmosphere due to photoreduction of sulfur-bound Hg(II). Some of this Hg(0) would deposit to the coastal ocean and thus drives an increase in coastal Δ¹⁹⁹Hg. However, even with the additional deposition of Hg(0) with positive Δ¹⁹⁹Hg, PZE (with a E<sub>occ(II-0),pho</sub> of 0.50‰) in the coastal ocean is still able to reproduce the observed Δ¹⁹⁹Hg shifts in coastal sections (Fig. 4). Some Hg(0) emitted by photoreduction of sulfur-bound Hg(II) could also deposit to soils and thus in theory cause a positive Δ¹⁹⁹Hg shift in soil reservoirs. This is indeed observed in two terrestrial sections in the Southern Hemisphere that exhibit positive excursions of Δ¹⁹⁹Hg during EPME intervals relative to pre-EPME strata<sup>60</sup>. However, there are much less terrestrial Hg isotope data across the EPME as compared to oceanic records, and the terrestrial Hg reservoirs are more likely affected by local sources and processes (e.g., local volcanism and weathering rates), and thus may not always record the positive Δ¹⁹⁹Hg shift (Supplementary Note 3). More terrestrial records are needed in the future to verify this hypothesis.
Since PZE is also proposed to enhance oceanic uptake of atmospheric Hg(0) by providing an additional pathway of Hg(0) oxidation and retention via bonding to sulfide or thiols in the surface ocean⁵⁴,⁵⁵, we evaluate the effect of such uptake on marine Hg enrichment and Δ¹⁹⁹Hg shift across the EPME. As shown in Supplementary Fig. 11, Hg(0) uptake into coastal oceans only results in very minor changes of marine Hg enrichment and Δ¹⁹⁹Hg values, while Hg(0) uptake into open oceans could significantly elevate Hg concentrations and decreases Δ¹⁹⁹Hg values. This is because the atmospheric Hg(0) uptake flux to the coastal ocean (~60 Mg a⁻¹) is more than 50 times smaller than Hg(0) uptake to the open ocean (~3150 Mg a⁻¹) (Supplementary Fig. 5). Increasing atmospheric Hg(0) to ~90% of total atmospheric Hg deposition flux to open oceans can negatively shift Δ¹⁹⁹Hg by >0.1‰ and double marine Hg enrichment factor in open oceans. Thus, the enhanced uptake of atmospheric Hg(0) due to PZE may also contribute to the negative shifts of Δ¹⁹⁹Hg across the EPME.

In summary, our modeling suggests that the MIF changes due to the development of PZE in coastal and open oceans can satisfactorily explain the extent of the negative shift of Δ¹⁹⁹Hg in sedimentary strata across the EPME in different marine depositional settings, but increased Hg emissions from volcanism are still needed to reach the observed Hg enrichment. This result does not exclude the contribution of soil erosion to the Hg enrichment and negative shift of Δ¹⁹⁹Hg, especially for the shallow marine Meishan section. Instead, it suggests that soil erosion alone (as proposed by most previous studies²⁴⁻²⁶,²⁹⁻³³) cannot account for the full extent of the Δ¹⁹⁹Hg negative shift during the main EPME interval as shown in the model Run #2, and additional Hg MIF in seawater caused by PZE is required for this shift.

**PZE is compatible with the timing delay between Hg enrichment and Δ¹⁹⁹Hg shift.** The PZE scenario is also compatible with the observed 40 kyr time lag between the largest negative Δ¹⁹⁹Hg shift and the spike of Hg enrichment in the Meishan section (Fig. 2c), because in this scenario the factors driving Hg enrichment (volcanism) and the negative Δ¹⁹⁹Hg shift (PZE) are different. In contrast, the peak timings of modeled Hg enrichment and Δ¹⁹⁹Hg shift are nearly coincident in Runs #1 and #2, further supporting our argument that neither volcanism nor soil erosion can fully account for the changes in Hg enrichment and Δ¹⁹⁹Hg during the EPME. The trigger for PZE across the EPME has been proposed to be environmental perturbations associated with volcanism (e.g., global warming and increased nutrient input due to enhanced continental weathering, Fig. 5 and Supplementary Fig. 1)⁷⁻⁹,¹¹,¹⁵. For example, terrestrial-derived nutrients such as phosphorous are well known to enhance marine primary productivity, which has been proposed as the key mechanism of oceanic anoxia during the EPME⁷,⁸. Widespread anoxia caused by phosphorous influx has been found to lag behind the onset of enhanced terrestrial weathering by as much as ~0.3 Myr in some locations across the EPME because phosphorous can be sequestered from the water column under non-sulfidic conditions but remobilized from the sediments under euxinic conditions⁷. A previous modeling study found that PZE started to develop if riverine phosphate fluxes during the EPME reached ~10 times the modern value⁹, suggesting that PZE may lag behind the onset of the enhanced terrestrial nutrient influx. Furthermore, rapid global warming during the EPME (Supplementary Fig. 1)⁶¹ can also lead to expansion of PZE due to decreased solubility of oxygen in warmer surface waters⁶², increased microbial metabolic rates⁶³, and enhanced hydrologic cycle and terrestrial weathering rates¹⁵,⁴⁸. Previous studies on Meishan section suggest that

seawater temperature⁶¹ (Fig. 1d) and terrestrial weathering (Fig. 1e) (based on lithium isotopes etc.)⁶⁴ due to volcanism started to increase from beds 22–23 (~252.10–252.03 Ma), well before the main extinction interval. In contrast, organic biomarker data at the Meishan section shows that PZE reached its first peak between bed 24 and 28 (~251.98–251.88 Ma)¹⁰,¹³, consistent with the timing of the main extinction interval as well as the peak of the negative Δ¹⁹⁹Hg shift (Fig. 1c). It is noted that Δ¹⁹⁹Hg values of the Meishan section started to decline again after 251.75 Ma, which likely reflects the second expansion of PZE at beds 30–37 as shown by biomarkers¹⁰,¹³, although the resolution of Hg MIF data is relatively low at these beds and thus more data is needed to verify this hypothesis (Supplementary Note 4). Additionally, δ²³UED in carbonates from global oceans (including South China)⁶³ indicate that the extent of global seafloor anoxia was the largest within the EPME. Taken together, this suggests that PZE may lag behind the onset of volcanic eruptions and terrestrial erosion, further supporting PZE as the main cause of the strong negative shift of Δ¹⁹⁹Hg within the main extinction interval that lagged behind the peak of Hg enrichment.

**Implications for past mass extinctions and future climate changes.** Our model suggests that volcanism (e.g., STLIP) was the main cause of the marine Hg enrichment, especially the strong Hg concentration spike immediately before the EPME, whereas soil erosion and PZE were the main drivers of the negative shift of marine Δ¹⁹⁹Hg across the EPME. In particular, the strong short-lived negative Δ¹⁹⁹Hg shift within the main EPME interval is best explained by PZE. These three factors are inherently related, as PZE was likely stimulated by global warming and enhanced inputs of terrestrial nutrients triggered by volcanism (Fig. 5). In Supplementary Fig. 12, a new scenario that combines volcanism, soil erosion, and PZE is modeled for the Meishan section. This scenario illustrates an interplay between volcanism, soil erosion, and PZE across the EPME, and suggests that the Δ¹⁹⁹Hg shift due to PZE could be better quantified once we have better constraints on volcanism and soil erosion.

Our modeling results have important implications for the use of Hg concentration anomalies and Hg isotopes in sedimentary rocks for understanding ancient volcanic eruptions, soil erosion, and ocean redox conditions, and for evaluating their impacts on biotic extinction or climate change events. First, our model not only confirms the previous finding that large-scale volcanism can cause Hg concentration anomalies, but also represents an important step towards quantitatively constraining the relationship between the volcanic Hg emission fluxes and Hg enrichment as well as Δ¹⁹⁹Hg variation in marine sediments. This quantitative framework can be applied to assess the volcanic perturbation to Hg cycling during other mass extinction events. Second, our model suggests that the marine Δ¹⁹⁹Hg shift is highly sensitive to PZE, which has been widely reported across the EPME and other mass extinction events but rarely recognized as a driver of Δ¹⁹⁹Hg shift in previous studies. Our model reinforces the previous argument based on experimental and observational data that Hg isotopes can be a novel proxy of oceanic redox conditions⁵⁰, but as a step forward, our model also sets the stage for a new way of quantifying the extent and persistency of PZE based on the degree and duration of Δ¹⁹⁹Hg shift, although further experiments are still needed to refine the quantitative relationship between Hg MIF and PZE.

Our model is inevitably associated with limitations and uncertainties. The Hg isotope model adopts a mechanistically simple box structure and algorithm, and is not spatially resolved. Thus, our analysis only focuses on the global scenario of each simulated reservoir. However, the compiled P-T sections used for
<figure><img src="image_7.png" /><figcaption>b) extinction</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>c) post-extinction</figcaption></figure>

<figure><img src="image_9.png" /></figure>

model evaluation were deposited in diverse marine environments, which show a high heterogeneity in Hg isotope compositions between sites. The numerous site-specific Hg isotope data needs to be reconciled under the global scenario for a better understanding of the factors driving Hg cycling. Second, the isotope enrichment factors and rate coefficients associated Hg

transfer between reservoirs are either based on laboratory-controlled experiments or empirically optimized against observations. These parameters may not be representative of natural ecosystems, and could vary considerably during mass extinction events. Therefore, more work is still needed to provide more robust constraints on the environmental changes across the
**Fig. 5 A conceptual model for PZE and Hg MIF across the EPME.** Shown is the development of PZE and its relationship with volcanism and soil erosion, and its impact on mass extinction, as well as how sediment Hg MIF and Hg/TOC change under different redox conditions. **a** Pre-extinction: the photic zone was oxic and the continental shelf (coastal ocean) harbored a thriving ecosystem. The Hg in sediments was dominated by background atmospheric deposition with a small positive Δ¹⁹⁹H (-0.1‰). **b** During the EPME: emplacement of large-scale volcanism (e.g., STLIP) increased Hg emission and soil erosion, leading to Hg enrichment (represented by the peak of Hg/TOC). Volcanism also led to the development of strong PZE by enhancing terrestrial nutrient influx (see the chain of environmental deteriorations), which killed marine life by shrinking the habitats on the shelf. PZE caused a strong negative shift of Δ¹⁹⁹H (by ~0.2‰) relative to the pre-extinction value) due to photoreduction of Hg(II)-S compounds and enhanced uptake of atmospheric Hg(0) by sulfidic seawater. **c** Post-extinction: with a decrease in volcanic activity, sediment Hg enrichment decreased but was still higher than the pre-extinction level. PZE was weakened due to the decline of terrestrial nutrient inputs, but likely sustained via recycling of nutrients (e.g., P) deposited previously⁷. The weaker PZE alongside less soil erosion led to a smaller negative shift of Δ¹⁹⁹H (-0.1‰) relative to the pre-extinction value).

EPME. Overall, our model provides the basis for further understanding of the interplay among environmental triggers (volcanism, climate change), soil erosion, oceanic redox, and biotic evolution in the geological past (Fig. 5). Understanding this interplay is increasingly critical to understanding Earth's habitability in the future, as marine oxygen minimum zones have been expanding over many continental shelves due to climate and anthropogenically induced eutrophication⁶⁵.

## Methods

**Definition of Hg enrichment factor and isotope ratio.** Hg enrichment factor in sedimentary strata (EF<sub>strata</sub>) is defined as

$$EF_{\text{strata}} = \frac{(\text{Hg/TOC})_{\text{EPME}}}{(\text{Hg/TOC})_{\text{pre-EPME}}} \quad (1)$$

Where Hg/TOC is the concentration ratio between Hg and total organic carbon (TOC) in samples across the EPME and before the EPME (i.e., the Late Permian background). Hg concentrations are commonly normalized to TOC to preclude Hg enrichment due to changes in organic burial rate²⁰. The enrichment factor in modeled reservoirs (EF<sub>model</sub>) is defined as the Hg mass ratio between individual reservoirs across the EPME and before the EPME. Theoretically, EF<sub>strata</sub> and EF<sub>model</sub> represent the same quantity.

Hg isotope composition is expressed as δ<sup>xxx</sup>Hg (‰, xxx = 199, 200, 201, 202, 204) by normalizing to a common NIST 3133 Hg standard:

$$\delta^{xxx}\text{Hg}(\‰) = [(^{xxx}\text{Hg}/^{198}\text{Hg})_{\text{sample}}/(^{xxx}\text{Hg}/^{198}\text{Hg})_{\text{NIST3133}} - 1] \times 1000 \quad (2)$$

MIF value is denoted as Δ<sup>xxx</sup>Hg (‰, xxx = 199, 200, 201, 204), representing the difference between the measured δ<sup>xxx</sup>Hg value and that predicted from δ<sup>202</sup>Hg using a kinetic MDF law⁶⁶:

$$\Delta^{xxx}\text{Hg}(\‰) = \delta^{xxx}\text{Hg} - \delta^{xxx}\beta \times \delta^{202}\text{Hg} \quad (3)$$

The mass-dependent scaling factor xxxβ is 0.2520 for ¹⁹⁹Hg, 0.5024 for ²⁰⁰Hg, 0.7520 for ²⁰¹Hg and 1.4930 for ²⁰⁴Hg.

**Model structure and algorithm, and time-emission framework.** The GBM-Hg-ISO-EPME model adopts the basic structure of the GBM-Hg-ISO-14box-v1 model of Sun et al.²² that was developed based on global Hg cycling box models of Harvard University⁶⁷. The reservoir sizes and fluxes used to derive the mass transfer coefficients are based on the well-known modern Hg cycle and taken from the recent GEOS-Chem Hg model (www.geos-chem.org) and the global 3-D simulation for Hg in the Massachusetts Institute of Technology ocean general circulation model (MITgcm)⁶⁸ (Supplementary Fig. 5). GBM-Hg-ISO-14box-v1 comprises 14 fully-coupled speciated Hg reservoirs (atmospheric Hg(0)atm/Hg(II)atm, fast/slow/armed terrestrial (soil) Hgft/Hgts/Hgla, open surface (0–100 m) oceanic Hg(0)ocs/Hg(II)ocs/Hg(P)ocs, open intermediate (100–1000 m) oceanic Hg(0)oci/Hg(II)oci/Hg(P)oci, open deep (>1000 m) oceanic Hg(0)ocd/Hg(II)ocd/Hg(P)ocd) with Hg inputs via volcanism and Hg removal via terrestrial and oceanic Hg burial. For GBM-Hg-ISO-EPME, we collapse the three terrestrial reservoirs into one terrestrial reservoir (Hg<sub>soil</sub>), and the three deep oceanic reservoirs into one deep oceanic reservoir (Hg<sub>ocd</sub>), because there is no process-related MIF among three terrestrial reservoirs and among three deep oceanic reservoirs. This also decreases the model's complexity and uncertainty. To evaluate the direct effect of terrestrial Hg inputs on coastal oceans (an average water depth ~90 m, including mainly continental shelf), we additionally add coastal Hg reservoir (Hg(0)occ/Hg(II)occ/Hg(P)occ) according to the box model of Liu et al.⁶⁹. Depending on the sedimentary settings and depths of the compiled P-T sections with reported Hg isotope compositions (Supplementary Table 1), we evaluate modeled oceanic reservoirs against different sections. In detail, we evaluate the modeled coastal ocean against the Meishan and Guryul Ravine sections which were located at shallow water depths (<100 m) of the continental shelf. The deeper (200–500 m) sections of Shangsi, Xiakou/Daxiakou, and Chaohu are evaluated against the simulated open intermediate ocean, and the bathyal to near-abyssal sections of Buchanan Lake and Gujo-Hachiman against the simulated open deep ocean.

Mass transfer of Hg between boxes (reservoirs) of the GBM-Hg-ISO-EPME model are controlled by first-order rate coefficients (k, a⁻¹), and are represented by a group of 13 coupled ordinary differential equations:

$$\frac{dM_i}{dt} = \sum_{j \neq i} (k_{j \to i} \times M_j) - \sum_{i \neq j} (k_{i \to j} \times M_i) + S \quad (4)$$

Where M<sub>i</sub> (Mg) is the Hg mass in box 'i', and M<sub>j</sub> (Mg) is the Hg mass in box 'j'. k<sub>j→i</sub> = F<sub>j→i</sub>/M<sub>j</sub> and k<sub>i→j</sub> = F<sub>i→j</sub>/M<sub>i</sub> are the first-order rate coefficients (calculated using reservoir sizes and fluxes, Supplementary Fig. 5 and Table 2), in which F<sub>j→i</sub> (Mg a⁻¹) is Hg flux from box 'j' to box 'i' and F<sub>i→j</sub> (Mg a⁻¹) is Hg flux from box 'i' to box 'j'. S (Mg a⁻¹) is external Hg input from the deep mineral reservoir by volcanic processes (0 for all reservoirs except for the atmosphere and deep ocean).

Odd-MIF between boxes of the GBM-Hg-ISO-EPME model is controlled by isotope enrichment factors of Δ¹⁹⁹Hg (E¹⁹⁹Hg, which relates to fractionation factor of odd-MIF (a¹⁹⁹Hg<sub>MIF</sub>) by E¹⁹⁹Hg (‰) = 1000 × (a¹⁹⁹Hg<sub>MIF-1</sub>)), and are also represented by a group of 13 coupled ordinary differential equations:

$$\frac{d(M_i \times \Delta^{199}\text{Hg}_i)}{dt} = \sum_{j \neq i} [k_{j \to i} \times M_j \times (\Delta^{199}\text{Hg}_j + \text{E}^{199}\text{Hg}_{(j \to i)})] - \sum_{i \neq j} [k_{i \to j} \times M_i \times (\Delta^{199}\text{Hg}_i + \text{E}^{199}\text{Hg}_{(i \to j)})] + S \times \Delta^{199}\text{Hg}_S \quad (5)$$

The background subaerial volcanic Hg emission flux in the Late Permian is assumed to be 300 Mg a⁻¹ according to Dal Corso et al.¹², which takes into account the active tectonic degassing at that time (e.g., pyroclastic eruptions and lava effusion of STLIP), and is also similar to modern fluxes⁷⁰,⁷¹. The same value is applied to background submarine volcanic Hg flux, which is the median of the ranges estimated in the recent UNEP report⁷¹. The Δ¹⁹⁹Hg values of subaerial volcanism are assumed to zero according to recent measurements of Siberian continental flood basals and other igneous rocks⁷²,⁷³. Before running the model, we first optimize the parameters at the natural steady-state (i.e., dM<sub>i</sub>/dt = 0; d(Δ¹⁹⁹Hg)/dt = 0) with background Δ¹⁹⁹Hg values of Earth's reservoirs in the latest Permian as constraints. We use Δ¹⁹⁹Hg values of Perrian coals (fossilized plants)⁴⁰ and oceanic sediments (0.10–0.20‰)²⁴–²⁶,²⁹ as the background constraints to represent Δ¹⁹⁹Hg values of coeval terrestrial and oceanic reservoirs, respectively (Supplementary Table 3; Supplementary Note 3). They broadly exhibit complementary Δ¹⁹⁹Hg signatures: open oceans receive important atmospheric Hg(II) deposition characterized by positive Δ¹⁹⁹Hg, and terrestrial soils mainly sequester atmospheric Hg(0) via the plant pump, and are characterized by negative Δ¹⁹⁹Hg. According to the optimization strategies of the GBM-Hg-ISO-14box-v1 model²² and the least squares method, the isotope enrichment factors of Δ¹⁹⁹Hg (E¹⁹⁹Hg) used in the model are either based on laboratory-controlled experiments or optimized against the observational background Δ¹⁹⁹Hg constraints for terrestrial and oceanic reservoirs (Supplementary Table 3).

It is still uncertain how to build a precise temporal connection between volcanism and the EPME biotic crisis. According to high-precision radioisotope data, the EPME broadly coincides with the contemporaneous STLIP volcanism which is segmented into three distinct magmatic stages (S1: > 252.24 ± 0.1 to 251.907 ± 0.067 Ma, characterized by initial pyroclastic eruptions followed by lava effusion; S2: 251.907 ± 0.067 to 251.483 ± 0.088 Ma, characterized by cessation of extrusion and onset of widespread sill intrusion; S3: 251.483 ± 0.088 to <251.354 ± 0.088 Ma, characterized by resumed lava extrusion after a ~ 420 kyr hiatus and continued sill intrusion)⁴,²⁸. As proposed by Burgess et al.³⁸, the initial emplacement pulse of the sill intrusion stage mostly likely triggered the EPME where the main extinction interval occurred within 60 ± 48 kyr (251.941 ± 0.037 to 251.880 ± 0.031 Ma), followed by protracted biotic recovery²⁸. In the coupled Hg-C cycling model of Dal Corso et al.¹², they suggested that elevated Hg emission mainly occurred during the STLIP sill intrusion S2 stage with a 420 kyr duration (251.98–251.56 Ma), within which the peak Hg emission lasted for ~1 kyr (251.950–251.949 Ma) immediately before the marine EPME starting at ~251.941 Ma. Our model adopts the similar time-emission framework as Dal Corso et al.¹², but delays the starting time of elevated Hg emission by ~20 kyr to 251.96 Ma (now well within the geochronological uncertainty of the STLIP S2 starting time) after examining the time-series of Δ¹⁹⁹Hg values in the Meishan section, which declined since ~251.96 Ma (Fig. 2).

Starting from the reservoir Hg budgets and Δ¹⁹⁹Hg values at the natural steady-state, different model scenarios are tested under the above time-emission
framework to examine how volcanic emissions and associated geochemical perturbations (increased soil erosion, change of marine redox states) affect Hg budgets and $\Delta^{199}$Hg values in various Earth surface reservoirs across the EPME. Our aim is to converge plausible source and process parameters that could best reproduce the long-term (~400 kyr) and short-lived (~1 kyr) variations in Hg enrichment and $\Delta^{199}$Hg values across the EPME. We mainly test the model output of coastal ocean particulate-bound Hg (Hg(P)\_occ) against the published Hg/TOC and $\Delta^{199}$Hg values of the Meishan section.

## Data availability

The authors declare that the main data supporting the findings of this study are available within the article and its Supplementary Information file, and they are available at https://doi.org/10.5281/zenodo.7800694.

## Code availability

The model codes are available at https://zenodo.org/badge/latestdoi/444653796.

Received: 31 January 2023; Accepted: 26 April 2023;

Published online: 08 May 2023

## References

1. Fan, J-x et al. A high-resolution summary of Cambrian to Early Triassic marine invertebrate biodiversity. *Science* **367**, 272–277 (2020).

2. Dal Corso, J. et al. Environmental crises at the Permian–Triassic mass extinction. *Nat. Rev. Earth Environ.* **3**, 197–214 (2022).

3. Zhang, H. et al. Felsic volcanism as a factor driving the end-Permian mass extinction. *Sci. Adv.* **7**, eabh1390 (2021).

4. Burgess, S. D. & Bowring, S. A. High-precision geochronology confirms voluminous magmatism before, during, and after Earth's most severe extinction. *Sci. Adv.* **1**, e1500470 (2015).

5. Svensen, H. et al. Siberian gas venting and the end-Permian environmental crisis. *Earth Planet. Sci. Lett.* **277**, 490–500 (2009).

6. Sanei, H., Grasby, S. E. & Beauchamp, B. Latest Permian mercury anomalies. *Geology* **40**, 63–66 (2012).

7. Schobben, M. et al. A nutrient control on marine anoxia during the end-Permian mass extinction. *Nat. Geosci.* **13**, 640–646 (2020).

8. Schobben, M. et al. Flourishing ocean drives the end-Permian marine mass extinction. *Proc. Natl. Acad. Sci. USA* **112**, 10298–10303 (2015).

9. Meyer, K. M., Kump, L. R. & Ridgwell, A. Biogeochemical controls on photic-zone euxinia during the end-Permian mass extinction. *Geology* **36**, 747–750 (2008).

10. Cao, C. et al. Biogeochemical evidence for euxinic oceans and ecological disturbance presaging the end-Permian mass extinction event. *Earth Planet. Sci. Lett.* **281**, 188–201 (2009).

11. Zhou, W. et al. Expansion of photic-zone euxinia during the Permian–Triassic biotic crisis and its causes: Microbial biomarker records. *Palaeogeogr. Palaeoclimatol. Palaeoecol.* **474**, 140–151 (2017).

12. Dal Corso, J. et al. Permo–Triassic boundary carbon and mercury cycling linked to terrestrial ecosystem collapse. *Nat. Commun.* **11**, 2962 (2020).

13. Grice, K. et al. Photic zone euxinia during the Permian–Triassic superanoxic event. *Science* **307**, 706–709 (2005).

14. Hays, L. E., Beatty, T., Henderson, C. M., Love, G. D. & Summons, R. E. Evidence for photic zone euxinia through the end-Permian mass extinction in the Panthalassic Ocean (Peace River Basin, Western Canada). *Palaeoworld* **16**, 39–50 (2007).

15. Meyer, K. M. & Kump, L. R. Oceanic euxinia in earth history: causes and consequences. *Annu. Rev. Earth Planet. Sci.* **36**, 251–288 (2008).

16. Whiteside, J. H. & Grice, K. Biomarker records associated with mass extinction events. *Annu. Rev. Earth Planet. Sci.* **44**, 581–612 (2016).

17. French, K. L., Rocher, D., Zumberge, J. E. & Summons, R. E. Assessing the distribution of sedimentary C40 carotenoids through time. *Geobiology* **13**, 139–151 (2015).

18. French, K. L. et al. Reappraisal of hydrocarbon biomarkers in Archean rocks. *Proc. Natl. Acad. Sci. USA* **112**, 5915–5920 (2015).

19. Summons, R. E., Welander, P. V. & Gold, D. A. Lipid biomarkers: molecular tools for illuminating the history of microbial life. *Nat. Rev. Microbiol.* **20**, 174–185 (2021).

20. Grasby, S. E., Them, T. R., Chen, Z., Yin, R. & Ardakani, O. H. Mercury as a proxy for volcanic emissions in the geologic record. *Earth Sci. Rev.* **196**, 102880 (2019).

21. Pyle, D. M. & Mather, T. A. The importance of volcanic emissions for the global atmospheric mercury cycle. *Atmos. Environ.* **37**, 5115–5124 (2003).

22. Sun, R. et al. Modelling the mercury stable isotope distribution of Earth surface reservoirs: Implications for global Hg cycling. *Geochim. Cosmochim. Acta* **246**, 156–173 (2019).

23. Blum, J. D., Sherman, L. S. & Johnson, M. W. Mercury isotopes in earth and environmental sciences. *Annu. Rev. Earth Planet. Sci.* **42**, 249–269 (2014).

24. Shen, J. et al. Evidence for a prolonged Permian–Triassic extinction interval from global marine mercury records. *Nat. Commun.* **10**, 1563 (2019).

25. Grasby, S. E. et al. Isotopic signatures of mercury contamination in latest Permian oceans. *Geology* **45**, 55–58 (2017).

26. Wang, X. et al. Mercury anomalies across the end Permian mass extinction in South China from shallow and deep water depositional environments. *Earth Planet. Sci. Lett.* **496**, 159–167 (2018).

27. Xie, S. et al. Changes in the global carbon cycle occurred as two episodes during the Permian–Triassic crisis. *Geology* **35**, 1083–1086 (2007).

28. Burgess, S. D., Bowring, S. & Shen, S. Z. High-precision timeline for Earth's most severe extinction. *P. Natl. Acad. Sci. USA* **111**, 3316 (2014).

29. Wang, X. et al. Global mercury cycle during the end-Permian mass extinction and subsequent Early Triassic recovery. *Earth Planet. Sci. Lett.* **513**, 144–155 (2019).

30. Shen, J. et al. Mercury fluxes record regional volcanism in the South China craton prior to the end-Permian mass extinction. *Geology* **49**, 452–456 (2021).

31. Sial, A. N. et al. Hg Isotopes and enhanced Hg concentration in the Meishan and Guryul Ravine Successions: Proxies for volcanism across the Permian–Triassic boundary. *Front. Earth Sci.* **9**, 651224 (2021).

32. Sial, A. N. et al. Globally enhanced Hg deposition and Hg isotopes in sections straddling the Permian–Triassic boundary: Link to volcanism. *Palaeogeogr. Palaeoclimatol. Palaeoecol.* **540**, 109537 (2020).

33. Georgiev, S. V. et al. Late Permian–Early Triassic environmental changes recorded by multi-isotope (Re-Os-N-Hg) data and trace metal distribution from the Hovea-3 section, Western Australia. *Gondwana Res.* **88**, 353–372 (2020).

34. Sun, R. et al. Historical (1850–2010) mercury stable isotope inventory from anthropogenic sources to the atmosphere. *Elem. Sci. Anth.* **4**, 000091 (2016).

35. Moynier, F. et al. The mercury isotopic composition of Earth's mantle and the use of mass independently fractionated Hg to test for recycled crust. *Geophys. Res. Lett.* **48**, e2021GL094301 (2021).

36. Grasby, S. E. et al. Progressive environmental deterioration in northwestern Pangea leading to the latest Permian extinction. *GSA Bull.* **127**, 1331–1347 (2015).

37. Percival, L. M. E. et al. Determining the style and provenance of magmatic activity during the Early Aptian Oceanic Anoxic Event (OAE 1a). *Glob. Planet. Change* **200**, 103461 (2021).

38. Burgess, S. D., Muirhead, J. D. & Bowring, S. A. Initial pulse of Siberian Traps sills as the trigger of the end-Permian mass extinction. *Nat. Commun.* **8**, 164 (2017).

39. Retallack, GregoryJ. & Jahren, A. H. Methane Release from Igneous Intrusion of Coal during Late Permian Extinction Events. *J. Geol.* **116**, 1–20 (2008).

40. Sun, R., Sonke, J. E. & Liu, G. Biogeochemical controls on mercury stable isotope compositions of world coal deposits: A review. *Earth Sci. Rev.* **152**, 1–13 (2016).

41. Deng, C. et al. Mercury isotopic compositions of the Precambrian rocks and implications for tracing mercury cycling in Earth's interior. *Precambrian Res.* **373**, 106646 (2022).

42. Sun, R., Grasby, S. E., Shen, J., Xiao, J. & Yin, R. Climate/ocean dynamics and possible atmospheric mercury depletion events during the Late Sturtian deglaciation. *Chem. Geol.* **598**, 120830 (2022).

43. Grasby, S. E., Sanei, H. & Beauchamp, B. Catastrophic dispersion of coal fly ash into oceans during the latest Permian extinction. *Nat. Geosci.* **4**, 104–107 (2011).

44. Ogden, D. E. & Sleep, N. H. Explosive eruption of coal and basalt and the end-Permian mass extinction. *Proc. Natl. Acad. Sci. USA* **109**, 59–62 (2012).

45. McElwain, J. C., Wade-Murphy, J. & Hesselbo, S. P. Changes in carbon dioxide during an oceanic anoxic event linked to intrusion into Gondwana coals. *Nature* **435**, 479–482 (2005).

46. Elkins-Tanton, L. T. et al. Field evidence for coal combustion links the 252 Ma Siberian Traps with global carbon disruption. *Geology* **48**, 986–991 (2020).

47. Yudovich, Y. E. & Ketris, M. P. Mercury in coal: a review: Part I. *Geochemistry. Int. J. Coal Geol.* **62**, 107–134 (2005).

48. Sephton, M. A. et al. Catastrophic soil erosion during the end-Permian biotic crisis. *Geology* **33**, 941–944 (2005).

49. Shen, S. Z. et al. Calibrating the End-Permian Mass Extinction. *Science* **334**, 1367–1372 (2011).

50. Zheng, W., Gilleaudeau, G. J., Kah, L. C. & Anbar, A. D. Mercury isotope signatures record photic zone euxinia in the Mesoproterozoic ocean. *P. Natl. Acad. Sci. USA* **115**, 10594–10599 (2018).
51. Zheng, W. & Hintelmann, H. Isotope fractionation of mercury during its photochemical reduction by low-molecular-weight organic compounds. *J. Phys. Chem. A* **114**, 4246–4253 (2010).

52. Motta, L. C., Kritee, K., Blum, J. D., Tsz-Ki Tsui, M. & Reinfelder, J. R. Mercury isotope fractionation during the photochemical reduction of Hg(II) coordinated with organic ligands. *J. Phys. Chem. A* **124**, 2842–2853 (2020).

53. Kritee, K., Motta, L. C., Blum, J. D., Tsui, M. T.-K. & Reinfelder, J. R. Photomicrobial visible light-induced magnetic mass independent fractionation of mercury in a marine microalga. *ACS Earth Space Chem.* **2**, 432–440 (2018).

54. Zheng, W., Liang, L. & Gu, B. Mercury reduction and oxidation by reduced natural organic matter in anoxic environments. *Environ. Sci. Technol.* **46**, 292–299 (2012).

55. Zheng, W., Lin, H., Mann, B. F., Liang, L. & Gu, B. Oxidation of dissolved elemental mercury by thiol compounds under anoxic conditions. *Environ. Sci. Technol.* **47**, 12827–12834 (2013).

56. Zheng, W. et al. Mercury stable isotope fractionation during abiotic dark oxidation in the presence of thiols and natural organic matter. *Environ. Sci. Technol.* **53**, 1853–1862 (2019).

57. Bergquist, B. A. & Blum, J. D. Mass-dependent and -independent fractionation of Hg isotopes by photoreduction in aquatic systems. *Science* **318**, 417–420 (2007).

58. Zheng, W. & Hintelmann, H. Mercury isotope fractionation during photoreduction in natural water is controlled by its Hg/DOC ratio. *Geochim. Cosmochim. Acta* **73**, 6704–6715 (2009).

59. Shen, Y. et al. Multiple S-isotopic evidence for episodic shoaling of anoxic water during Late Permian mass extinction. *Nat. Commun.* **2**, 210 (2011).

60. Shen, J. et al. Mercury evidence from southern Pangea terrestrial sections for end-Permian global volcanic effects. *Nat. Commun.* **14**, 6 (2023).

61. Joachimski, M. M. et al. Climate warming in the latest Permian and the Permian–Triassic mass extinction. *Geology* **40**, 195–198 (2012).

62. Kennett, J. P. & Stott, L. D. Abrupt deep-sea warming, palaeoceanographic changes and benthic extinctions at the end of the Palaeocene. *Nature* **353**, 225–229 (1991).

63. Hülse, D. et al. End-Permian marine extinction due to temperature-driven nutrient recycling and euxinia. *Nat. Geosci.* **14**, 862–867 (2021).

64. Sun, H. et al. Rapid enhancement of chemical weathering recorded by extremely light seawater lithium isotopes at the Permian–Triassic boundary. *P. Natl. Acad. Sci. USA* **115**, 3782–3787 (2018).

65. Schobben, M. et al. Eutrophication, microbial-sulfate reduction and mass extinctions. *Commun. Integr. Biol.* **9**, e1115162 (2016).

66. Blum, J. D. & Bergquist, B. A. Reporting of variations in the natural isotopic composition of mercury. *Anal. Bioanal. Chem.* **388**, 353–359 (2007).

67. Amos, H. M. et al. Global Biogeochemical Implications of Mercury Discharges from Rivers and Sediment Burial. *Environ. Sci. Technol.* **48**, 9514–9522 (2014).

68. Zhang, Y. et al. Biogeochemical drivers of the fate of riverine mercury discharged to the global and Arctic oceans. *Global Biogeochem. Cycles* **29**, 854–864 (2015).

69. Liu, M. et al. Rivers as the largest source of mercury to coastal oceans worldwide. *Nat. Geosci.* **14**, 672–677 (2021).

70. Li, C. et al. Unequal anthropogenic enrichment of mercury in Earth's Northern and Southern Hemispheres. *ACS Earth Space Chem.* **4**, 2073–2081 (2020).

71. UN-Environment. Global Mercury Assessment 2018. (UN Environment Programme, Chemicals and Health Branch Geneva, Switzerland, 2019).

72. Deng, C. et al. Mercury isotopic composition of igneous rocks from an accretionary orogen: Implications for lithospheric recycling. *Geology* **50**, 1001–1006 (2022).

73. Yin, R. et al. Mantle Hg isotopic heterogeneity and evidence of oceanic Hg recycling into the mantle. *Nat. Commun.* **13**, 948 (2022).

## Acknowledgements

This study is supported by the National Key Research and Development Plan (2022YFF0800300), National Natural Science Foundation of China (42173011; 41973009; 41922019; 41625012; 41961144028; 42293280; 42261144668; 42073002) and National Science Foundation of United States (award No. EAR 1760203 to A.D.A). We thank Benjamin J. W. Mills and Jacopo Dal Corso for kindly providing us with the chronologically calibrated Hg/TOC and Hg isotope data of the Meishan section used in their Hg-C isotope model. Dominik Hülse is thanked for sharing with us their compiled data on δ²³U of carbonate and seawater temperature.

## Author contributions

R.S. and W.Z. conceived and designed this project. R.S. developed the box model code with the help of Y. Zhang. R.S. and W.Z. led manuscript writing with major supports from J.E.S. and A.D.A. Y.L., F.Z., Y. Zhao, J.C., C-Q.L., and S.S. contributed to the data interpretation.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s43247-023-00821-6.

**Correspondence** and requests for materials should be addressed to Wang Zheng.

**Peer review information** *Communications Earth & Environment* thanks the anonymous reviewers for their contribution to the peer review of this work. Primary Handling Editor: Joe Aslin. Peer reviewer reports are available.

**Reprints and permission information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_10.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2023