Permo-Triassic boundary carbon and mercury cycling linked to terrestrial ecosystem collapse.

Records suggest that the Permo-Triassic mass extinction (PTME) involved one of the most severe terrestrial ecosystem collapses of the Phanerozoic. However, it has proved difficult to constrain the extent of the primary productivity loss on land, hindering our understanding of the effects on global biogeochemistry. We build a new biogeochemical model that couples the global Hg and C cycles to evaluate the distinct terrestrial contribution to atmosphere-ocean biogeochemistry separated from coeval volcanic fluxes. We show that the large short-lived Hg spike, and nadirs in δ202Hg and δ13C values at the marine PTME are best explained by a sudden, massive pulse of terrestrial biomass oxidation, while volcanism remains an adequate explanation for the longer-term geochemical changes. Our modelling shows that a massive collapse of terrestrial ecosystems linked to volcanism-driven environmental change triggered significant biogeochemical changes, and cascaded organic matter, nutrients, Hg and other organically-bound species into the marine system.

Introduction

The Permo–Triassic mass extinction (PTME) is the largest known extinction in Earth′s history, with the loss of ~90% of species in the sea and ~70% of species on land[1-4]. The PTME has been causally linked to the emplacement of the Siberian Traps Large Igneous Province (LIP) and associated volcanic gas emissions (especially CO2, SO2 and halogens), via widespread environmental changes such as warming and oceanic anoxia[5-8]. The PTME also saw a crisis in terrestrial ecosystems, with loss of plant diversity, increased wildfire activity and consequent enhanced soil erosion[9-14]. Recent work has shown that the disruption of vegetation started before and culminated at the marine extinction level (Fig. 1), implying that the environmental disaster impacted terrestrial ecosystems first[12,14]. The cause of the terrestrial mass extinction is still unclear, and several kill mechanisms have been hypothesised. For example, a shift from a humid warm climate to an unstable highly seasonal climate and an associated increase in wildfires affected the equatorial Permo–Triassic peatlands, drastically reducing the abundance and diversity of the flora[14]; abnormal pollen and spores found in different localities around the world during the PTME interval suggest widespread mutagenesis possibly linked to an increase in UV-B radiation due to ozone depletion[15,16]; a terrestrial S-isotope record from the Karoo basin in South Africa could indicate volcanically driven acid rain at the P–T transition[17] that might have also severely impacted the flora. Whilst the taxonomic losses in terrestrial ecosystems are becoming clearer[12,18], and local, enhanced input of terrestrial material into marine environments has been recorded[9,11], the biogeochemical impacts and feedbacks on the exogenic C cycle are not known. It is possible that these impacts were severe: the PTME represents the largest, and maybe the only, known mass extinction of insects[19], suggesting that there may have been a substantial decrease in available food sources at the lowest levels of the food chain.

Fig. 1

Permian–Triassic events.

Existing records show the terrestrial extinction started before, and its collapse culminated at, the marine extinction interval, where large Hg and Hg/TOC spikes, negative δ202Hg, negative Δ199Hg and more negative δ13C values (EP. I[30]) are recorded. After this event, the δ13C record shows a second minimum (Ep. II), while Hg concentration decreases, but remains relatively higher than background values, and δ202Hg and Δ199Hg rebound towards more positive values. This pattern is seen in many P–T successions, in terrestrial and shallow- to deep-water settings. Meishan section: U-Pb ages from Burgess et al.[34]; TOC and C isotopes (δ13C) from Xie et al.[30]; Hg concentrations, δ202Hg and Δ199Hg data are from Grasby et al.[23] and Wang et al.[20]. Chinahe section: all data and correlation to Meishan section are from Chu et al.[14]. Note that to plot Hg (ppb) and Hg/TOC data, the same scale is used. MDF mass-dependent fractionation, MIF mass-independent fractionation, EP. episode.

The PTME is marked by an approximately two- to four-fold increase in marine sedimentary Hg concentration with respect to background levels during a ~400 kyr interval also characterised by negative δ13C values[20], which implies a relatively long-term injection of Hg and 13C-depleted CO2 into the atmosphere–land–ocean system during this time. Superimposed on this trend is a prominent, short-lived Hg spike, which is usually expressed as Hg/TOC, given the affinity of Hg with organic matter, and which is coincident with the collapse of the terrestrial ecosystems[14], the onset of the marine mass extinction interval and a sharp minimum in δ13C values (Fig. 1). While higher Hg/TOC values have been reported preceding the marine extinction in the deep-water settings of Japan[21], these are an artefact of normalisation to values of TOC which are below the analytical detection limit (<<0.1%), and do not track increases in Hg concentrations (see TOC and Hg data in the Supplementary Information of ref. [21]). The Hg record is interpreted as evidence of increased Hg input into the Earth's surface system from the Siberian Traps[13,21-24]. However, Hg is also stored in large quantities in terrestrial biomass and soils[25-27]. The mobilisation of these terrestrial reservoirs during the PTME could also result in the increased loading of Hg to aquatic environments, even without an elevated volcanic Hg flux to the atmosphere[23,27-29]. Any such changes in the soil and biomass carbon reservoir would also have direct implications for the release of C to the atmosphere and the sedimentary δ13C record.

Hg isotopes can be used to better understand how Hg has been transported into the sedimentary environment. Mass-independent fractionation (MIF—denoted by Δ) occurs due to aqueous photoreduction of Hg2+ to that takes place in surface waters and in clouds[29]. This results in positive Δ199Hg values in the remaining water Hg2+ pool, hence positive Δ199Hg values in sediments dominated by atmospheric Hg2+ deposition. Conversely, slightly negative Δ199Hg values characterise the terrestrial reservoir (soil and biomass), which primarily captures [29]. During Hg uptake by plants additional mass-dependent fractionation (MDF) and MIF occur, resulting in more negative δ202Hg and Δ199Hg values[29]. The Permian–Triassic boundary shallow-water record of Meishan shows a prominent negative δ202Hg excursion in correspondence to the Hg spike coupled to a small negative Δ199Hg shift, but deeper water successions show persistent positive Δ199Hg[20,21,23,24]. Therefore, isotope data appear to indicate that the Hg was transported to the deep-water settings mostly via the atmosphere, and to the shallow-water settings via continental runoff and via the atmosphere[20,21,23,24].

The δ13C records at the Permian–Triassic boundary show two minima[30,31], which are here called EP. I and II (EP. = episode) following ref. [30] (Fig. 1). The observed δ13C trends are similarly recorded in different depositional settings and by different substrates (carbonates, bulk organic matter, separate plant remains)[14,30-32], strongly indicating that they represent actual changes in the C-isotope composition of the reservoirs of the exogenic C cycle. The Hg concentration spike occurs in the same interval as the minimum in δ13C values associated with the PTME (EP. I, Fig. 1). At Meishan, where the chronostratigraphic framework is well established[33,34], the initial large Hg spike occurs about 60 Kyrs after the onset of the carbon-isotope perturbation (Fig. 1). Data from non-marine end-Permian successions confirms this diachrony[13,14] (Fig. 1) and show that the Hg spike is also coincident with a sudden decrease in Total Organic Carbon (TOC) values to almost zero[14] (Fig. 1). The combination of geochemical and palaeontological data from these sections shows that the terrestrial ecosystem disruption started with the onset of the carbon-isotope perturbation and climaxed at the very sharp δ13C minimum (EP. I), coincident with the Hg spike (Fig. 1), and the start of the main marine extinction interval. This mismatch in both timings and fluxes between the C and Hg cycles at the PTME, suggest that the 13C-depleted C and the Hg came from multiple sources.

Overall, the interplay between volcanism and terrestrial reservoir changes in controlling PTME biogeochemistry is not well known, and previous attempts to model the δ13C record have been fundamentally hampered because of the lack of an independent tracer of the C source. To overcome this, we use published records of δ13C and Hg systematics to jointly constrain a new coupled C and Hg biogeochemical model. Our model shows that the large, sudden geochemical shifts at the PTME are best explained by a massive pulse of terrestrial biomass oxidation, while Siberian Traps volcanism can explain the longer-term geochemical changes.

Results and discussion

A coupled Hg-C cycle model

Figure 2 shows the biogeochemical box model developed here. The full model derivation follows in the ‘Methods’ section. The model combines a multi-box sediment–ocean–atmosphere carbon-alkalinity cycle (based on previous work[35-37]), with the global mercury cycle[38,39]. The ocean is split into ‘surface’, ‘high-latitude’ and ‘deep’ boxes. It considers ocean circulation and carbonate speciation, and contains a simplified organic carbon cycle in which burial rates are prescribed. As well as computing the global C and Hg cycles it also computes δ13C of all C reservoirs and δ202Hg of the ocean reservoirs. A full atmosphere–ocean model of δ202Hg would require dynamic biosphere reservoirs, which would greatly increase model complexity. We therefore simplify the system to a mixing model for marine δ202Hg, in which atmospheric and riverine inputs have different isotopic signatures. The atmosphere is assumed to have ‰, and riverine input is assumed to have ‰, following ref. [40].

Fig. 2

Model schematic.

Concentrations of modelled species are tracked in boxes representing the atmosphere (a), surface ocean (s), high-latitude ocean (h) and deep ocean (d). Exchange between boxes via air–sea exchange, circulation and mixing, or evasion and deposition are shown as dashed arrows. Biogeochemical fluxes between the hydrosphere and continents/sediments are shown as solid arrows. See text, ‘Methods' and Supplementary Information for full details of fluxes. a C cycle. b Mercury cycle.

Concentrations of modelled species are tracked in boxes representing the atmosphere (a), surface ocean (s), high-latitude ocean (h) and pan class="Chemical">deep ocean (d). Exchange between boxes via air–sea exchange, circulation and mixing, or evasion and deposition are shown as dashed arrows. Biogeochemical fluxes between the hydrospn>here and continents/sediments are shown as solid arrows. See text, ‘Methods' and Supn>plementary Information for full details of fluxes. a C cycle. b n>an class="Chemical">Mercury cycle.

The model is set up for the late Permian by reducing the solar constant to that of 250 Ma, and increasing the background tectonic n>an class="Chemical">CO2 degassing rate to 1.5 times the present day (D = 1.5), in line with estimates for the Late Permian[41]. To obtain the observed pre-event ocean–atmosphere δ13C composition of ~3.5‰, we set the rate of land-derived organic carbon burial to 20% higher than present day, and adjust the composition of the weathered carbonate reservoir to 3‰. This is consistent with high terrestrial productivity and C burial in the Permian (e.g., coal forests and mires) and rapid recycling of more recently buried and 13C-enriched carbonate material.

In the following paragraphs, we test two model end-member scenarios: (I) the release of volcanic and thermogenic pan class="Chemical">Hg and C from Siberian Traps activity alone, and (II) with the additional release of n>an class="Chemical">Hg and C as a consequence of the collapse of the terrestrial ecosystems.

Volcanic and thermogenic degassing

We first model the release of volcanic/volcanogenic Hg and C from the Siberian Traps. Existing radioisotope data show that the extinction, the negative δ13C excursion and Hg spike might have all occurred during the intrusive phase of the Siberian Traps[20,34,42]. It is suggested that the emplacement of large sills caused the combustion or thermal decomposition of organic-rich sediments with the consequent release of thermogenic volatiles, such as C and Hg[20,22]. It has been proposed that over a ~400 Kyrs intrusive phase the Siberian Traps emitted ~7600–13,000 Mg yr−1 of volcanic Hg, which included both magmatic and coal-derived Hg[20,22,23]. Relating this Hg release to the background volcanic source is not straightforward because estimates of the background source vary, but taking the most likely present day range[25] (~90–360 Mg yr−1), and further constraining this by taking into account the need to balance overall sedimentary burial of Hg (~190 Mg yr−1), and the overall ~50% increase in tectonic degassing in the late Palaeozoic relative to today[43], we arrive at a best guess for the background late Permian Hg flux of ~300 Mg yr−1. This means that the Siberian Traps eruption increased the geogenic Hg input by a factor of ~25–43 over ~400 Kyrs.

To test this scenario, we model Siberian Traps intrusion by increasing the volcanic Hg source by 25–43-fold for 450 Kyrs, while also increasing the CO2 source in line with estimates[44,45] for Siberian Traps degassing based on magma volumes and sediment intrusion (by 4-8 × 1012 moles/year). The CO2 released by contact metamorphism at the PTME is assumed to have an average δ13C composition of −25‰[44]. Specifically, the input functions are:Here the first vector is time in millions of years and the second is the flux alteration at that time. Here is the additional CO2 release in mol yr−1, and is the relative Hg degassing rate increase. For the duration of these pulses, the thermohaline circulation is also assumed to collapse due to warming and freshwater input[46]. We reduce the circulation term to 1 Sv over this period, which allows more rapid change in the model surface ocean C isotopes and Hg loading. This is a large reduction, and also reflects the simple structure of the model in which the entire low-latitude surface ocean is represented by a single box, and so is well-mixed. Figure 3a–d shows that this magnitude and timing of release of C and Hg is capable of driving the longer-term decline in carbonate δ13C, and the coeval long-term approximately two- to four-fold enrichment in shallow sediment Hg/TOC that is observed in Meishan. However, the model scenario does not capture the spike in Hg concentration, or nadir in δ13C (EP. I[30] in Fig. 1) that are coincident with the final stage of the terrestrial extinction. It also does not produce any substantial change in marine δ202Hg isotopes (Fig. 1d), because the primary Hg source to the ocean is the atmosphere for the full model run.

Fig. 3

Model results.

a–d Volcanic C and Hg source only (thermogenic plus magmatic). e–h With addition of ‘vegetation loss′ Hg source and corresponding increase in organic carbon oxidation. C-isotope, Hg/TOC and δ202Hg data from ‘Meishan GSSP' section[20,23]. In both columns, top panel (a, e) shows input forcing in terms of relative rate of Hg input, next panel (b, f) shows δ13C of new shallow ocean carbonate, third panel (c, g) shows Hg/TOC data versus the model molar Hg/Corg burial rate, which has been scaled to match the background as the model does not calculate weights. Final panel (d, h) shows model δ202Hg against data, where solid line is shallow water (box s) and dashed line is deep water (box d) Blue dots = published data; Red and green lines = model Hg fluxes; Pink lines = model output.

a–d Volcanic C and Hg source only (thermogenic plus magmatic). e–h With addition of ‘vegetation loss′ n>an class="Chemical">Hg source and corresponding increase in organic carbon oxidation. C-isotope, Hg/TOC and δ202Hg data from ‘Meishan GSSP' section[20,23]. In both columns, top panel (a, e) shows input forcing in terms of relative rate of Hg input, next panel (b, f) shows δ13C of new shallow ocean carbonate, third panel (c, g) shows Hg/TOC data versus the model molar Hg/Corg burial rate, which has been scaled to match the background as the model does not calculate weights. Final panel (d, h) shows model δ202Hg against data, where solid line is shallow water (box s) and dashed line is deep water (box d) Blue dots = published data; Red and green lines = model Hg fluxes; Pink lines = model output.

Within the model, we have also explored a scenario wherein the large Hg pulse represents a further rapid pulse of LIP volcanism. We have attempted this scenario in Supplementary Note 1 (Scenario I–2), where a 1 Kyr volcanic pulse is assumed to raise the Hg and C input rates by a further factor of 5. While the Hg/TOC can indeed be explained by an additional short-lived pulse of Hg, we require the total release rate of Hg to be ~200 times greater than background levels, and even then, this scenario fails to reproduce any of the Hg isotope signature or the nadir in carbonate δ13C (Supplementary Fig. 1).

Terrestrial ecosystem collapse

For scenario II, we explore the additional effects of a geologically rapid (~1 Kyr) pulse of Hg and C as the result of the collapse of terrestrial ecosystems at the PTME. The magnitude of this Hg flux is again difficult to quantify precisely, and we explore an increase of 100-fold over background conditions. This level of increase represents the magnitude required to drive the sedimentary Hg signal that we observe, and is compatible with the available terrestrial biosphere Hg reservoir: total soil Hg is estimated to be on the order of ~106 Mg Hg when considering a soil depth of ~15 cm[47]. So, our model Hg delivery flux would require decimetre-scale soil organic matter oxidation over 1000 years, coincident with the PTME and the sharp EP. I negative δ13C shift[11]. The Hg pulse is delivered directly to the low-latitude surface ocean via runoff in the model, and is accompanied by a pulse of ‘soil oxidation’ C which we assume raises the global rate of oxidative weathering by a factor of 30—a number chosen to have the observed level of impact on the C-isotope record while being compatible with the Hg input change. We also assume a cessation of terrestrial organic C burial. Terrestrial Hg deposition and erosion is not altered during the pulse as the fluxes are minor by comparison. The new model functions applied in addition to the longer-term inputs of scenario I are:

Here, the first vector is time in millions of years, and the second shows flux multipliers at these times. and denote the relative rate of land organic C burial, oxidative weathering and pan class="Chemical">Hg runoff, respn>ectively, and are set at 0, 30 and 100, respn>ectively, for the duration of the 1-kyr pulse.

This ‘biosphere’ pulse causes a short-term large concentration spike in the shallow marine Hg reservoir and its sediments, which is superimposed on the volcanically driven changes (Fig. 3e–h). The Hg spike is far larger than would be expected from simply increasing the volcanic source by the same amount because the biospheric Hg is delivered directly to the surface ocean and sedimentation occurs mostly on the shelf. With the inclusion of terrestrial C oxidation and cessation of terrestrial carbon burial, the model also replicates the transient shift to more negative δ13C values recorded at the marine extinction interval (EP. I[30,31] in Fig. 1): Terrestrial C oxidation is a source of isotopically light C[48]. The model now also shows a sharp negative δ202Hg shift in the shallow ocean box, which is triggered by increased Hg riverine input, but shows no change in the deeper ocean, where the source of Hg remains predominantly atmospheric. This also compares well with existing records (Fig. 3). At Meishan, which was located in the margins of the Yangtze carbonate platform, the Hg and Hg/TOC spike is coincident with more negative δ202Hg values (Fig. 1), while in the deeper water sections of south China the values are more positive[20,23].

Hence, oxidation of terrestrial biomass is a compelling scenario to explain the palaeontological, sedimentological and geochemical data. There is clear observational evidence for the collapse of the terrestrial ecosystems and cessation of terrestrial C burial, stratigraphic evidence supporting the sequence and timing of the events (onset of the δ13C shift—collapse of the terrestrial ecosystem—Hg and C spike), sufficient quantity of Hg available, consistency with the isotopic evidence for changing Hg sources, and consistency with the δ13C records.

Massive terrestrial biomass oxidation during the PTME

Using our coupled C–Hg biogeochemical model, we show that the massive collapse of terrestrial ecosystems and oxidation of terrestrial biomass during the Permian–Triassic extinction had a huge impact on global Hg and C biogeochemistry. Hg stored in the terrestrial reservoirs was rapidly released as a consequence of the loss of terrestrial biomass and increased soil erosion[9,14]. This mechanism is the best explanation for the sharp increased loading of Hg into both terrestrial and marine water bodies and the negative shift in δ202Hg in coincidence with the marine mass extinction. Contemporaneously, increased soil carbon oxidation introduced large quantities of isotopically light C, accounting for the sharp negative δ13C anomaly registered in the sedimentary record (EP. I[30]). In the model, the emission of Hg and C from magma and heating of sedimentary organic matter during the intrusive phase of the Siberian Traps LIP emplacement can account for the smaller, two- to fourfold increase of Hg concentrations with respect to background levels, and the relatively longer negative δ13C trend that is recorded by both carbonates and organic matter, in marine and terrestrial settings.

A new scenario emerges for the PTME that links the collapse of ecosystems on land to the global geochemical changes recorded at the marine extinction interval. The disruption of terrestrial environments started during the initial phases of the Siberian Traps emplacement likely due to the release of volcanic gases as CO2, SO2 and halogens, which could have triggered acid rain, ozone depletion, volcanic darkness, rapid cooling and subsequent global warming[8,49]. At the culmination of the terrestrial disturbance interval, when the ecosystems totally collapsed, large amounts of 13C-depleted C and Hg deriving from a massive oxidation of terrestrial biomass were transported into aqueous habitats causing a steep decline in sedimentary δ13C (carbonates and organic matter), a sedimentary Hg concentrations spike and a shift in δ202Hg (Fig. 3). At this level, the marine mass extinction started. This, according to the existing chronostratigraphic framework, happened ∼60 Kyrs after the onset of the carbon-isotope perturbation and of the terrestrial ecological disturbances[14] (Fig. 1).

The biogeochemical cycle of Hg is intimately linked to the cycle of organic matter and its constituting elements, such as C, N, S and n>an class="Chemical">P[50]. Hence, besides Hg and C, other organically-bound species would have been transferred from the terrestrial reservoirs into the marine system in large quantities at EP. I (Fig. 1). Addition of these species, particularly the nutrients P and N, are easily capable of driving ecosystem turnover, anoxia and eutrophication, and it is likely that this terrestrial input contributed to the marine extinction[9,11]. Our model does not include these additional cycles, but other models have shown that a relatively small increase in marine P delivery (2–3-fold) has the potential to drive marine anoxia or euxinia[51,52]. The scale of the terrestrial ecosystem collapse at the PTME could explain the severity of the biotic crisis at the Permian–Triassic boundary at all trophic levels, and should be a key consideration for future research. For other events, the Hg records are not so consistent nor as detailed as for the PTME. However, it is very likely that future research on other intervals could show the same Hg and C patterns as for the PTME.

Methods

Model derivation

This model is designed to track the transfer and isotopic signature of atmospheric and marine carbon and n>an class="Chemical">mercury over geological time, while being broadly applicable to changes on the timescale of ocean circulation. The biogeochemical system is taken largely from ref. [36], with some additions from refs. [37,53,54], with the underlying hydrological model from ref. [35]. The Hg cycle follows ref. [25].

Model structure

The model has three ocean boxes: surface (s), high latitude (h) and deep (d). As in ref. [35], the surface box is 100-m deep and occupies 85% of the ocean surface, whereas the high-latitude box is 250-m deep and represents 15% of the ocean surface. Each ocean box includes the same biogeochemical species, and a pan class="Chemical">thermohaline circulation mixes the boxes in the order s, h, d. The upn>per boxes exchange with the atmospn>here, which is a single box. As well as transfer fluxes between ocean and atmospn>here boxes, biogeochemical fluxes of weathering, degassing and burial operate between the surface system and crust.

Model species

All model species are shown in Table 1.

Table 1

Model species.

Description	Name	Exists in	Size at present
Surface ocean water	W_s	Surface ocean	3.07×10^16 m3
High-latitude water	W_h	High latitude	1.35×10^16 m3
Deep water	W_d	Deep ocean	1.35×10^18 m3
Atmospheric CO2	CO_2a	Atmosphere	5×10^16 mol C
Surface ocean DIC	DIC_s	Surface ocean	6×10^16 mol C*
High-latitude DIC	DIC_h	High latitude	3×10^16 mol C*
Deep ocean DIC	DIC_d	Deep ocean	3×10^18 mol C*
Surface ocean alkalinity	ALK_s	Surface ocean	6×10^16 mol CaCO3 equiv.*
High-latitude alkalinity	ALK_h	High latitude	3×10^16 mol CaCO3 equiv.*
Deep ocean alkalinity	ALK_d	Deep ocean	3×10^18 mol CaCO3 equiv.*
δ13C of atmospheric CO2	δ^13CO_2a	Atmosphere	−7 ‰
δ13C of surface ocean DIC	δ^13DIC_s	Surface ocean	0.1 ‰**
δ13C of high-latitude DIC	δ^13DIC_h	High latitude	0.1 ‰**
δ13C of deep ocean DIC	δ^13DIC_d	Deep ocean	0.1 ‰**
Atmospheric Hg	Hg_a	Atmosphere	3.5×10^6 mol Hg***
Surface ocean Hg	Hg_s	Surface ocean	1.34×10^7 mol Hg***
High-latitude Hg	Hg_h	High latitude	5.9×10^6 mol Hg***
Deep ocean Hg	Hg_d	Deep ocean	5.9×10^8 mol Hg***

*Starting values chosen close to equilibrium values, model equilibrates to DIC 2 mM and ALK  2.2 mM.

**Note that 0.1 is used instead of 0 to increase model stability at initialisation *** after Amos et al.[25] and scaled to ocean box volumes.

Model fluxes

Model fluxes, with equations and present values are shown in Table 2.

Table 2

Model fluxes.

Description	Name	Equation	Size at present
Transfer fluxes	tran_ij	C_if_circ	Multiple
Air–sea exchanges	f_airsea_j	A_jM_atm(pCO_2a−pCO_2jτ_oa)	Multiple
Silicate weathering	f_silw	k_baswf_T_bas+k_granwf_T_gran	8×10^12. mol C yr−1
Carbonate weathering	f_carbw	k_carbwf_Tcarb	8×10^12 mol C yr−1
Oxidative weathering	f_oxidw	k_oxidwRO_2^0.5F_oxidw	7.75×10^12 mol C yr−1
Carbonate degassing	f_ccdeg	k_ccdegD	8×10^12 mol C yr−1
Organic carbon degassing	f_ocdeg	k_ocdegD	1.25×10^12 mol C yr−1
Mare carbonate burial	f_mccb	k_mccbΩ−1^1.7Ω_0	16×10^12 mol C yr−1
Marine organic C burial	f_mocb	k_mocb	4.5×10^12 mol C yr−1
Marine carbonate burial	f_mccb	k_mccb	16×10^12 mol C yr−1
Land organic C burial	f_locb	k_locbF_Cburial	4.5×10^12 mol C yr−1
Volcanic Hg release	f_Hg_volc	k_Hg_volcF_Hg_input	1.5×10^6 mol Hg yr−1*
Wildfire Hg release	f_Hg_wildfire	k_Hg_wildfire	5×10^5 mol Hg yr−1*
Riverine Hg runoff	f_Hg_runoff	k_Hg_runoffF_runoff	2×10^6 mol Hg yr−1*
Marine Hg burial	f_Hg_b	k_Hg_b(Hg_sHg_s_0)	2×10^6 mol Hg yr−1*
Vegetation Hg deposition	f_vegdep	k_Hg_vegdep(Hg_aHg_a_0)	1×10^7 mol Hg yr−1*
Vegetation Hg evasion	f_vegeva	k_Hg_vegeva	1×10^7 mol Hg yr−1*
Ocean Hg deposition	f_oceandep_j	k_Hg_oceandep(Hg_aHg_a_0)	1.5×10^7 mol Hg yr−1*
Ocean Hg evasion	f_oceaneva_j	k_Hg_oceaneva(Hg_jHg_j_0)	1.5×10^7 mol Hg yr−1*

*Volcanic Hg flux follows ref. [25], and is assumed equivalent to the total Hg burial in their model to close the system over long timescales. This value is within their stated reasonable range but is ~2× their chosen volcanic flux, the value is further increased by 50% due to higher rates of tectonic degassing through the Permian and Triassic. Runoff and burial fluxes are calculated to close the system under volcanic and wildfire input. Deposition and evasion fluxes directly follow ref. [25].

*Volcanic Hg flux follows ref. [25], and is assumed equivalent to the total n>an class="Chemical">Hg burial in their model to close the system over long timescales. This value is within their stated reasonable range but is ~2× their chosen volcanic flux, the value is further increased by 50% due to higher rates of tectonic degassing through the Permian and Triassic. Runoff and burial fluxes are calculated to close the system under volcanic and wildfire input. Deposition and evasion fluxes directly follow ref. [25].

Non-flux calculations

Atmospheric pan class="Chemical">CO2 volume ratio is calculated as:where is atmospn>heric n>an class="Chemical">CO2 in moles, and is this value at present day.

Global average surface temperature (GAST) is:where kclim is climate sensitivity to doubling pan class="Chemical">CO2, and tgeol is time in millions of years before present and is expn>ressed in negative terms. Low-latitude surface temperature (Ts) is assumed to scale by times global temperature change, and both high-latitude (Th) and deep (Td) temperature are assumed to follow global temperature change.

For pan class="Chemical">carbonate spn>eciation, effective equilibrium constants are calculated following refs. [36,55]:

Dissolved pan class="Chemical">carbon species are then calculated following Walker and Kasting[36]:

We also explicitly calculate [H+] concentration to observe model pH:

Calcium carbonate saturation state is calculated as:where Ωj is the n>an class="Chemical">CaCO3 saturation state in box j, and Ksp is the solubility product.

For terrestrial chemical weathering, temperature dependence of basalt and granite weathering is calculated as:

And temperature dependence of pan class="Chemical">carbonate weathering:

Fixed parameters

Fixed parameters are shown in Table 3.

Table 3

Fixed model parameters.

Description	Name	Size at present
Thermohaline speed	f_circ	10 Sv
Relative area of low-latitude surface ocean	A_s	0.85
Relative area of high-latitude surface ocean	A_h	0.15
Present day moles of atmospheric CO2	M_atm	5×10^16 mol C
Timescale parameter for gas exchange	τ_oa	10 years
Long-term climate sensitivity	k_clim	5 K
CO2 second dissociation constant	K_2	7.4×10^−10
Calcium carbonate solubility product	K_sp	0.8 mmol2 kg−2 *
Present day CaCO3 saturation state	Ω_0	3
Marine calcium concentration	[Ca]	10 mM
Photosynthetic C fractionation	ΔB	27‰
C-isotope composition of carbonates	δ^13C_C	0‰
C isotope composition of organics	δ^13C_G	−27‰

*Chosen within ocean range (0.43–1.15)[56] to achieve reasonable DIC and ALK at present.

Differential equations

The following equations track the 11 non-pan class="Chemical">water spn>ecies from Table 1.

Atmospheric pan class="Chemical">CO2:

Low-latitude surface ocean DIC:

High-latitude surface ocean DIC:

pan class="Chemical">Deep ocean DIC:

Low-latitude surface ocean alkalinity:

High-latitude surface ocean alkalinity:

pan class="Chemical">Deep ocean alkalinity:

δpan class="Chemical">13C of atmospheric pan class="Chemical">CO2:

δpan class="Chemical">13C of low-latitude surface ocean DIC:

δpan class="Chemical">13C of high-latitude surface ocean DIC:

δpan class="Chemical">13C of pan class="Chemical">deep ocean DIC:

Atmospheric pan class="Chemical">Hg:

Low-latitude surface ocean pan class="Chemical">Hg:

High-latitude surface ocean pan class="Chemical">Hg:

pan class="Chemical">Deep ocean pan class="Chemical">Hg:

δ202pan class="Chemical">Hg of Low-latitude surface ocean pan class="Chemical">Hg:

δ202pan class="Chemical">Hg of High-latitude surface ocean pan class="Chemical">Hg:

δ202Hg of n>an class="Chemical">Deep ocean Hg: