Enhanced clay formation key in sustaining the Middle Eocene Climatic Optimum

Abstract

The Middle Eocene Climatic Optimum (around 40 million years ago) was a roughly 400,000-year-long global warming phase associated with an increase in atmospheric CO₂ concentrations and deep-ocean acidification that interrupted the Eocene’s long-term cooling trend. The unusually long duration, compared with early Eocene global warming phases, is puzzling as temperature-dependent silicate weathering should have provided a negative feedback, drawing down CO₂ over this timescale. Here we investigate silicate weathering during this climate warming event by measuring lithium isotope ratios (reported as δ⁷Li), which are a tracer for silicate weathering processes, from a suite of open-ocean carbonate-rich sediments. We find a positive δ⁷Li excursion—the only one identified for a warming event so far —of ~3‰. Box model simulations support this signal to reflect a global shift from congruent weathering, with secondary mineral dissolution, to incongruent weathering, with secondary mineral formation. We surmise that, before the climatic optimum, there was considerable soil shielding of the continents. An increase in continental volcanism initiated the warming event, but it was sustained by an increase in clay formation, which sequestered carbonate-forming cations, short-circuiting the carbonate–silicate cycle. Clay mineral dynamics may play an important role in the carbon cycle for climatic events occurring over intermediate (i.e., 100,000 year) timeframes.

Main

During the course of the Cenozoic Earth’s surface temperatures have exhibited a long-term reduction^1,2. However, amidst the backdrop of this general cooling trend, the Earth experienced a number of warming phases on million year timescales, such as the Early Eocene Climatic Optimum³. Global warming events on timescales of tens of thousands of years (for example, the Palaeocene–Eocene Thermal Maximum or PETM) were accompanied by a massive input of ¹³C-depleted carbon and deep-sea carbonate dissolution due to ocean acidification, leading to a shoaling of the carbonate compensation depth (CCD)^4,5,6. Carbon cycle dynamics across such long- and short-term warming events are relatively well understood^7,8. However, a warming phase on an intermediate timescale (i.e. 100,000’s of years)—termed the Middle Eocene Climatic Optimum (MECO)—around 40 million years ago (Ma) ^9,10,11 has raised interest because it may be inconsistent with the current model of climate–carbon cycle interactions¹². Although the MECO exhibits similar characteristics to the hyperthermal events of the early Eocene, such as increases in atmospheric CO₂ ($p_{{\mathrm{CO}}_{2}}$) and temperature, perturbation of the hydrological cycle¹³ and temporary shoaling of the CCD^10,11,14, its protracted nature (which has been estimated to last between 270 and 500 thousand years (kyr)^9,12,15) appears to suggest that different dynamics were in operation^12,16. Over geologic timescales, the chemical weathering of silicate rocks acts as a negative feedback on rising $p_{{\mathrm{CO}}_{2}}$ and temperature levels, thus regulating the climate⁷. Across the MECO, silicate weathering should have kept pace with increases in $p_{{\mathrm{CO}}_{2}}$ and temperature, leading to an enhanced supply of metal cations and alkalinity to the oceans, with a consequent increase in carbonate burial and a reduction in $p_{{\mathrm{CO}}_{2}}$ levels¹². While other greenhouse warming events were also characterized by transient ocean acidification and shoaling of the CCD (and the lysocline, that is, the depth at which calcite dissolution increases dramatically), the timing from the onset of these events to peak acidification was much more rapid than for the MECO¹⁷. Thus, the MECO currently stands out as a unique event in the Cenozoic^10,12,18.

An increase in CO₂ input to the atmosphere from enhanced volcanic degassing is one hypothesized trigger for the MECO⁹. Conventional carbon cycle theory predicts that, on the timescale of the MECO, injections of volcanic CO₂ should lead to a deepening of the CCD, since the additional CO₂ should be converted to alkalinity via continental weathering. However, the available data suggest that the opposite occurred¹². Osmium (Os) isotope records, combined with carbon cycle box modelling, suggest that a diminished silicate weathering feedback operated during the middle Eocene, enabling the accumulation of volcanic CO₂ in the exogenic carbon pool¹⁶. Indeed, a substantial proportion of the present-day Earth has low weathering rates due to soil shielding or aridity¹⁹. Here we aim to test this hypothesis further by assessing silicate weathering dynamics during the MECO using lithium isotope ratios (δ⁷Li) in marine sediments, followed by model simulations using a newly developed biogeochemical box model.

The dissolution of primary rock shows no lithium isotopic fractionation, whereas the formation of secondary minerals (for example, clays) preferentially take up the light ⁶Li isotope, leaving residual surface waters enriched with the heavier ⁷Li isotope^20,21,22,23. Thus, the δ⁷Li composition of rivers can provide information on the weathering congruency—the ratio of primary rock dissolution versus secondary mineral formation—at that particular time. The importance of clay formation lies not only in its cation-exchange capacity but also in its ability to aid the preservation of buried organic matter, since both organic matter and clays can exhibit a negative surface charge that attracts Mg²⁺ and Ca²⁺, leading to their retention in soils as well as changes in their supply in the dissolved load to the oceans^24,25,26. Here we show that this balance between primary mineral dissolution and secondary clay formation is critical for resolving the MECO conundrum.

Middle Eocene lithium isotope ratios

We investigate the response of silicate weathering and clay formation during the MECO by obtaining δ⁷Li data from the same carbonate-rich sediments deposited in pelagic settings that were previously investigated for osmium isotope ratios (¹⁸⁷Os/¹⁸⁸Os)¹⁶, which are together representative of the global ocean. These include the Ocean Drilling Program (ODP) Site 1263 in the South Atlantic (Fig. 1a,b), the Integrated Ocean Drilling Program (IODP) Site U1333 in the equatorial Pacific (Fig. 1a,c) and ODP Site 959 in the equatorial Atlantic (Fig. 1a and Supplementary Fig. 1). We updated existing age models¹⁶ to convert the core depth to age (based on the 2020 Geologic Time Scale (GTS2020)²⁷ and plotted the data for Sites 1263 and U1333 against both depth and age). While Site 959 exhibits the same secular trend for both the lithium and osmium isotopes as Sites 1263 and U1333, indicating that it is recording global changes, its sediment lithology is somewhat anomalous, mainly porcellanite, and the absolute δ⁷Li values are markedly different; thus, we do not include this site in further discussions (see Supplementary Information for more details).

Fig. 1: Lithium and osmium isotope records from (I)ODP Sites 1263 and U1333.

a, Sample site locations and the possible palaeogeography at that time⁶¹. The general trend of the middle Eocene to the middle Oligocene was one of decreasing flooded continental area^42,62,63,64. b,c, Lithium and osmium isotope records for Site 1263 (b) and Site U1333 (c), where the filled blue circles represent the mean of each sample (n = 3) and the error bars represent the ±2σ precision on each δ⁷Li analysis; the ¹⁸⁷Os/¹⁸⁸Os data are from ref. ¹⁶. The grey band represents the estimated duration of the MECO based on oxygen isotope (δ¹⁸O) data, from a number of ODP sites, tied to an age model framework based on ODP Site 702 (ref. ¹⁸). Other data previously used¹⁶ to distinguish the MECO (that is, δ¹⁸O (ref. ¹⁰) and CaCO₃ (ref. ⁶⁵)) are also plotted. All data are plotted against both the adjusted revised metres composite depth (rmcd) and age (GTS2020). Uneven temporal spacing is due to variations in sedimentation rates.

Both Sites 1263 and U1333 exhibit a positive lithium isotope excursion (LIE) of ~3–3.5‰ across the MECO, with Site 1263 increasing from a pre-MECO value of ~22.1‰ to a peak of 25.6‰ during the event, and U1333 moving from ~23.1‰ to 26.1‰. We combined the isotopic data from Sites 1263 and U1333 and then applied a smoothing spline to obtain a general trend for our timeframe of interest for modelling purposes (41.5–39.0 Ma; Fig. 2 and Methods). Based on an average fractionation factor between seawater and carbonate (Δ⁷Li_sw-carb) of 4‰ (2–6‰ full range^28,29,30,31), we surmise that, before the MECO, the δ⁷Li of seawater (δ⁷Li_sw), based on the smoothed average of the two sites, was ~26.1–27.8‰ (full range 24.1–29.8‰), which is in line with previous estimates²⁰, and rose to a peak of ~29.7‰ (27.7–31.7‰) between 40.18–40.10 Ma. Seawater pH declined during the MECO¹⁸, and although pH has been shown to affect δ⁷Li fractionation somewhat³², the change in pH would have resulted in only a minor change in δ⁷Li fractionation, suggesting that an impact on our dataset is unlikely. Furthermore, the absolute effect of pH on δ⁷Li fractionation is still contested³³.

Fig. 2: The general trend of the lithium and osmium isotope records across the MECO.

a,b, The data from Sites 1263 and U1333 have been combined (as with Fig. 1, the filled symbols represent the mean ± 2σ) and smoothing spline fits (mean ± 2σ) have been applied for δ⁷Li_carb (a) and ¹⁸⁷Os/¹⁸⁸Os_initial (b)¹⁶. The grey dotted lines indicate the length of the MECO (~40.425 to ~40.023 Ma) based on the δ¹⁸O data from a number of ODP sites tied to an age model framework based on ODP Site 702¹⁸. The data are plotted against time (GTS2020).

Given the ~1 Myr residence time of lithium in the present ocean²⁰, the shift in δ⁷Li within a few 100 kyrs indicates, geologically speaking, a fairly rapid and strong change in the Earth system conditions. Consistent with the long residence time is the observation that, by 39 Ma, the δ⁷Li of carbonates (δ⁷Li_carb) is still higher than before the MECO. The beginning of the positive LIE appears to occur at approximately the same time as the negative osmium isotope excursion (OsIE), with both isotope systems exhibiting changes at ~40.5 Ma just before the start of the MECO, while the seawater ¹⁸⁷Os/¹⁸⁸Os_initial (reflecting the osmium isotope ratio of marine sediments at the time of their deposition) recovers to pre-MECO values at, or just after, 40 Ma. The faster recovery time of the OsIE may be explained by the shorter residence time of osmium in the ocean (5–54 kyr)³⁴. Furthermore, the opposite signs of the LIE and OsIE are contrary to previous work, where both isotope systems change in tandem; whereas a positive LIE has hitherto been associated with cooling events, warming events have produced negative LIEs^29,35,36. A negative OsIE could be explained by an increased flux of high-temperature hydrothermal fluids, which would supply unradiogenic osmium to the ocean^16,34. However, this would also deliver isotopically light δ⁷Li (refs. ^20,37), causing a negative LIE, in contrast to our results from all three (I)ODP sites.

Modelling Earth system changes during the MECO

To assess the background Earth system state before the MECO and the changes that occurred during the event, we explored a number of scenarios using our newly developed CARLIOS biogeochemical box model (see Methods and Supplementary Information for further details about the model). For all scenarios we attempted to approximate the $p_{{\mathrm{CO}}_{2}}$ and temperature proxy data (see Supplementary Information for details on the proxies) as well as our δ⁷Li and ¹⁸⁷Os/¹⁸⁸Os data, and to reproduce some general characteristics of the MECO, such as a decrease in the surface ocean pH and a shoaling of the lysocline (Table 1).

Table 1 CARLIOS model simulations assessed against key MECO metrics

All scenarios (1–8) in Table 1 can capture one or more of the main MECO characteristics (Fig. 3 and Extended Data Figs. 1–7); however, although Scenarios 1, 3 and 4, for example, produce an increase in the $p_{{\mathrm{CO}}_{2}}$ (and temperature) and a decrease in the ocean pH, none of these scenarios shoal the lysocline, and Scenarios 1 and 4 do not produce a good fit to the δ⁷Li and ¹⁸⁷Os/¹⁸⁸Os data. An increase in normalized uplift only (Scenario 2) does raise the lysocline depth and matches the δ⁷Li data, but reduces $p_{{\mathrm{CO}}_{2}}$ while the surface pH increases, contrary to the proxy data. Shifting more of the carbonate burial from deep sediments to the continental shelves (Scenario 5), to reflect platform flooding due to thermal expansion of seawater, barely makes a difference to the background (pre- and post-MECO) results of the model, thus in previous modelling¹² this parameter change may have made little difference to the results. Combining the parameters of Scenarios 1 and 2, we model an injection of CO₂ directly to the atmosphere that occurs concurrently with a big increase in uplift (Scenario 6). Although this raises both $p_{{\mathrm{CO}}_{2}}$ and the lysocline, and approximately matches our δ⁷Li data, the pH also increases. None of the changes in Scenarios 1–6 produce a match to the ¹⁸⁷Os/¹⁸⁸Os data. For Scenario 7, we again added CO₂ to the atmosphere, albeit with a different functional form to Scenarios 1 and 6 (compare Extended Data Figs. 6a and 7a), modelled a relatively modest increase in uplift, allowed the global erosion-to-silicate-weathering ratio to vary (this ratio is used to help convert normalized uplift to erosion in tonnes per year), and increased the flux of unradiogenic osmium from continental eruptions (see Supplementary Information for more details about the osmium cycle). Scenario 7 produced results which were a good fit to almost all key aspects of the MECO; however, the lysocline deepened. Thus, for Scenario 8 we used all of the parameters from Scenario 7 but hypothesized that there may have been two additional changes during the MECO: a decline in the concentration of magnesium in seawater ([Mg]_sw), which affects the Mg/Ca ratio of seawater and the calcite solubility product (see below); and a reduction in the total areal exposure of carbonate rocks to weathering¹². This scenario (Fig. 3) successfully captured the key aspects of the MECO as listed in Table 1. Although the calcite saturation state (Ω_calcite) in this scenario does increase slightly, it remains at <1, meaning that carbonate dissolution will still outweigh sedimentation, and hence we use Scenario 8 as our chosen scenario for further discussion.

Fig. 3: Key parameter changes to the CARLIOS model and its results for Scenario 8.

a–g, Parameter changes for: an injection of CO₂ to the atmosphere (a), a change in the normalized uplift (b), a change to the global erosion-to-silicate-weathering ratio (E:W) (c), a pulsed continental eruption of unradiogenic osmium which is then subaerially weathered (d), scalings for the burial of marine organic carbon (F_bg-sea, blue), reverse weathering (F_rw, dashed gold) and the carbonate land area (f_L, red) (e), a change to the Mg/Ca ratio in the oceans ([Mg]_sw/[Ca]_sw) (f) and the factor that apportions carbonate burial (F_bc) to shelf environments, where a lower number indicates more burial on shelves (g). h–p, Key model results for atmospheric $p_{{\mathrm{CO}}_{2}}$ (black) plotted against proxy estimations from boron isotope (δ¹¹B) data¹⁸ (see Supplementary Information) (h), global average surface temperature (GAST) plotted against proxy estimations from δ¹⁸O data (see Supplementary Information) (i), ocean pH (for a low-latitude surface (gold), a high-latitude surface (blue) and a deep ocean (green), and where the red dots and line denote individual pH estimations from δ¹¹B data¹⁸ and a smoothed fit, respectively) (j), modelled lysocline depth (k), deep-ocean Ω_calcite saturation (l), the modelled mean δ⁷Li_sw offset by a fractionation factor (3–5‰, dark grey band; 2–6‰, light grey band) (m), the modelled weathering intensity through time (n), the contribution of rivers versus the combined sinks to the δ⁷Li_sw signature (o) and the modelled ¹⁸⁷Os/¹⁸⁸Os_initial of seawater (p). In all panels (apart from e, j and m) the black line is the model average and the grey band denotes ±1σ. In m and p, the blue line and band are from Fig. 2. ALK, ocean alkalinity.

A global shift towards enhanced clay formation

Our model scenarios indicate that, as with previous work¹², no single mechanism in isolation can reproduce the carbon cycle characteristics of the MECO and match the δ⁷Li and ¹⁸⁷Os/¹⁸⁸Os data. However, individual scenarios are useful for understanding the various impacts on the Earth system. For example, Scenario 2 indicates that the global weathering intensity (W/D)—the ratio of chemical weathering W to total denudation D (which is W plus erosion E)—decreased across the MECO, implying that there was a shift in the Earth’s weathering regime³⁰ from clay dissolution to clay formation, on a global average.

Calcium (Ca) and magnesium cations are attracted to the negatively charged surface layers of some clays and are structural components of others (for example, smectites). With increased terrestrial clay formation as suggested by the δ⁷Li data, it follows that more calcium and magnesium will be retained in soils, leaving a deficit of oceanic calcium and magnesium to form carbonates and close the carbonate–silicate cycle^24,26. Variable calcium and magnesium concentrations are observed in modern rivers (for example, ref. ²²), and more work needs to be conducted on the relationship between riverine δ⁷Li and cation concentrations. Many marine clay-formation reactions (reverse weathering) also take up calcium and magnesium, as well as bicarbonate ions, releasing CO₂ to the oceans³⁸. It is currently unknown by how much terrestrial calcium and magnesium retention scales with clay formation, and thus the subsequent impact on atmospheric CO₂ is equally unknown. It has been shown that magnesium is less mobile than calcium³⁹ and is therefore trapped by clays more efficiently; magnesium is also taken up in higher molar ratios during reverse weathering reactions³⁸. Thus, clay formation, either on land or in the oceans, will probably impact upon the magnesium cycle and its availability in the oceans for carbonate-forming reactions more than calcium. Although magnesium currently has a long residence time (~13 Myr), [Mg]_sw at ~40 Ma was much lower than today (~38.4 mmol compared with ~53 mmol)⁴⁰. The inclusion of this reduction in the Mg/Ca ratio of seawater in our model (Scenario 8) helps to shoal the lysocline (compared with Scenario 7) because it affects the calcite solubility product, which then affects the calcite saturation with subsequent knock-on effects for the carbonate system in the oceans.

We plotted the weathering intensity (W/D) and δ⁷Li of rivers (δ⁷Li_riv) results for Scenario 8 against the ‘Dellinger boomerang’²¹ (Fig. 4a) to determine more precisely the changes in the global W/D regime. Before the onset of the MECO, owing to a relatively low background erosion rate⁴¹, a steadily declining sea level and areal extent of flooded continent⁴², our modelling suggests that the Earth predominately had a mixture of mature floodplains in the high latitudes (Fig. 4a,b, position 4) with areas of tropical weathering (Fig. 4a,b, position 5), although the Earth would still have had mountainous areas with rapid erosion (Fig. 4b, position 1). This may have resulted globally in a net amount of secondary mineral dissolution³⁰. Data from the continental interior of North America⁴³ and from South America⁴⁴ suggest that warm and wet conditions prevailed at the beginning and end of the MECO. Fully coupled climate model simulations for 38 Ma, validated by field reconstructions, suggest near-surface air temperatures of >20 °C from the 30° to 50° latitudes and >30 °C from the equator to the ±30° latitudes⁴⁵. This would indicate that the background state before (and after) the MECO provided ideal conditions for extensive tropical weathering, although, just as in the present day, there would have been some areas that experienced aridity (for example, refs. ^19,46). Hence, we postulate there was, globally, more intense weathering, although more climatic studies of the MECO are required to validate this hypothesis.

Fig. 4: Possible changes to the Earth system as inferred by δ⁷Li data and model results.

a, The modern-day relationship between W/D and δ⁷Li_riv (the boomerang shape), with which our modelled average W/D and δ⁷Li_riv results from Scenario 8 are compared, as well as those inferred from the modelling of other climatic events (where the W/D and δ⁷Li_riv data decline for Oceanic Anoxic Event 2 (OAE2) and the PETM (red arrow) but increase for the Hirnantian glaciation (blue arrow)); the present-day global average is also shown^{21,29,30,35,36}. b, The inset shows the pre- and post-MECO conditions (grey and blue triangles near position 4), the peak MECO conditions (yellow cross near position 3) and the present-day conditions (yellow star) and how they correspond to the different surface environments (main image) that affect δ⁷Li_riv values³⁰.

An increase in CO₂ input to the ocean-atmosphere from volcanic CO₂ degassing, as well as from organic carbon oxidation and sulfide oxidation (leading to sulfuric acid weathering of carbonates) due to surface exposure after erosion⁴⁷, would have increased $p_{{\mathrm{CO}}_{2}}$ levels and consequently global surface temperatures. This in turn would probably have resulted in changes to the hydrological cycle and thermal expansion of surficial bodies of water. Several papers have suggested a sea-level rise across the MECO (for example, ref. ⁴²) and several recorded transgressive sediment sequences seemingly correspond to the MECO, but the age uncertainties associated with all of these data are greater than the duration of the MECO itself⁴⁸. Nevertheless, an enhanced hydrological cycle¹³ and thermal expansion would have increased both the areal extent of silicate weathering and the water–rock interaction time via the formation of new floodplains⁴⁹. At the same time, thermal expansion may have flooded shelf environments, minimizing the land area available for carbonate weathering but increasing the areal extent of shallow seas from which carbonates could precipitate (for example, ref. ¹² and references therein), which would have reduced the alkalinity supply and thus carbonate accumulation rates in the deep ocean.

An increase in volcanic activity globally during the MECO would have increased the global erosion rate due not only to a warming-driven increase in rainfall frequency and intensity but also from incision of the underlying bedrock by lava flows, and a combination of the two (increased precipitation and bedrock incision) can manifest itself in lahar events^50,51,52. Such an increase in erosion rates would have provided fresh primary minerals to the newly formed floodplains, which in turn provided an ideal environment for weathering and the subsequent formation of secondary minerals (Fig. 4b, position 3), leaving residual riverine waters enriched in the ⁷Li isotope^21,30.

An increase in secondary mineral formation would have led to more entrainment of Ca²⁺ and Mg²⁺ within these locales, thus reducing the flux of carbonate-forming cations to the oceans. This likely led to a disruption in the carbonate–silicate cycle between the dissolution of terrestrial silicate rocks and the formation of marine carbonates (that is, silicate weathering became less efficient⁵³), resulting in a positive feedback loop whereby the retention of Ca²⁺ and Mg²⁺ in the terrestrial realm led to an accumulation of CO₂ in the atmosphere. Consequently, rates of silicate weathering increased, liberating more Ca²⁺ and Mg²⁺ to be taken up by secondary minerals. The positive feedback loop may only have abated once the supply of primary minerals had declined sufficiently such that clay dissolution once more outweighed clay formation, leading to a net release of Ca²⁺ and Mg²⁺ back to the oceans and potentially to an overshoot of the CCD, although the currently available data are not of sufficient resolution to be able to evaluate this proposition. Alternatively, a short-lived pulse of CO₂ centred around the MECO ‘peak’, as suggested by the pH decline¹⁸, may have induced rapid hydrological changes such that the global weathering regime became congruent and the cation flux to the oceans increased considerably. Our model suggests that the global average weathering regime on Earth during the MECO was not dissimilar to that at present (Fig. 4). Assuming that the amount of magnesium and calcium in soils before the MECO is the same as at present (6.10 × 10¹⁴ t), then by the end of the MECO, assuming a duration of 400 kyr, we calculate that the amount of magnesium and calcium retained in soils would be 6.227 × 10¹⁴ t, an increase of ~2%, which we believe is a feasible amount of cations stored (see Supplementary Information for derivation). Reverse weathering may also have had an important part to play in terms of reducing the concentration of magnesium and calcium in the oceans while also contributing more CO₂ to the ocean-atmosphere system^38,54. Our model suggests that, pre- and post-MECO, the influence of sinks outweighed rivers (rivers/sinks < 1; Fig. 3o), but during the MECO the rivers had a greater control on the δ⁷Li signature of the oceans (rivers/sinks > 1). However, this is not to say that the sinks of lithium declined; indeed, there is evidence for an increase in abyssal smectite precipitation during the MECO^46,55, but proportionally these sinks were dwarfed by riverine changes.

All previously analysed global warming events, both long-term (OAE2)²⁹ and transient (PETM)³⁶, exhibit a negative LIE, which is interpreted as a relatively greater increase in erosion rates over weathering rates (causing a decrease in W/D). This is attributed to an accelerated hydrological cycle, consistent with modern river observations (Fig. 4a)⁵⁶. Hence, the MECO is the only warming event identified so far with a positive LIE. However, the mechanisms underlying changes to the weathering regime remain the same across all warming events analysed so far: the erosion rate increases relative to the weathering rate. The direction of the LIE then depends on the pre-event weathering regime (Fig. 4), which for the MECO is at a vastly higher W/D than the other hyperthermal events analysed so far.

Our results appear to indicate that silicate weathering continued at pace but was dwarfed by the increase in erosion. Global net secondary mineral formation increased, trapping carbonate-forming cations on the continents, effectively making the link between silicate weathering and $p_{{\mathrm{CO}}_{2}}$ drawdown less efficient. In addition, warmer deep-ocean waters possibly increased the reverse weathering rate, trapping more metal cations and exchanging bicarbonate ions for CO₂. Consequently, only a modest injection of CO₂ would have been required to kickstart the MECO, which is in line with other work¹⁸, and would help explain why there is no global negative δ¹³C excursion associated with the onset of the MECO^12,57. A spike in CO₂ degassing during the MECO, as modelled in Scenarios 1, 6, 7 and 8, may explain the rise in $p_{{\mathrm{CO}}_{2}}$ suggested by δ¹¹B (ref. ¹⁸). Once the fresh supply of primary minerals started to decline, the warm climate^14,45 likely facilitated a return to net clay dissolution, providing the oceans with some Ca²⁺ and Mg²⁺, leading to an increase in carbonate burial¹⁰. Importantly, a flux of partially dissolved clays to the oceans may have provided microsites for enhanced marine organic carbon burial, whilst also liberating clay-bound PO₄³⁻ that would have stimulated primary productivity, resulting in an overall decline in $p_{{\mathrm{CO}}_{2}}$ (ref. ⁵⁸). In addition, it has been shown that bacterial metabolic rates are sensitive to temperature changes, with a greater efficiency for remineralizing organic matter at higher temperatures^59,60. This may have created a positive feedback, which accelerated both MECO warming and post-MECO cooling.

For CO₂ drawdown to be determined correctly, changes in the silicate dissolution rate must be tempered with changes to the weathering regime. It is the latter that not only changes due to rapid alteration of the climate but is also highly dependent on the pre-perturbation palaeogeography, climate and vegetation regimes. While this is true for all climatic events, it appears that on intermediate timescales, such as those of the MECO, the dynamics of clay minerals may play a particularly important role in the global carbon cycle, effectively disrupting the negative feedback system that regulates $p_{{\mathrm{CO}}_{2}}$ and temperature. Hence, we recommend further investigations of the impact of clay dynamics on intermediate-timescale events in a model with more complex terrestrial and marine carbonate and clay chemistries.

Methods

Samples

Carbonate samples were analysed as detailed in previous work^31,36. Briefly, the samples were leached in 1 M sodium acetate buffered to pH 5 with acetic acid at room temperature. A split (~20%) of the subsequent leachate was analysed via quadrupole inductively coupled plasma mass spectrometry, using matrix-matched calibration standards and the JLs-1 reference material as an accuracy standard, to determine the element/Ca ratios. This is in part to determine ratios such as Al/Ca and Rb/Ca, to ascertain that negligible silicate material has been leached. Cutoffs for the silicate contribution on Al/Ca ratios are suggested as <0.8 mmol mol⁻¹ (ref. ²⁹) and on Al/(Ca + Mg) as <0.45 mmol mol⁻¹ (ref. ⁶⁶). Meanwhile, a cutoff of around <30 μmol mol⁻¹ for Rb/Ca has been suggested⁶⁷. All of our samples fit within these parameters (see Supplementary Tables 2 and 3). The remaining 80% of the sample was purified through a two-step cation-exchange method, using 0.2 M HCl as the eluent⁶⁸, and analysed relative to the lithium isotopic reference material IRMM-016 using a Nu Plasma 3 multi collector ICP mass spectrometer (Nu Instruments) at the LoGIC laboratories of University College London. All samples were renormalized to the LSVEC reference material, and the analytical uncertainty was propagated to account for this⁶⁹. Using an Aridus II desolvator, a signal intensity of ~130 pA (~13 V) on ⁷Li for a solution of 5 ng ml⁻¹ was achieved (using 10¹¹ Ω resistors and an uptake rate of ~100 μl min⁻¹, and between 1 and 5 ng lithium was analysed), which is much greater than the background and total procedural blank of 0.01 pA (<1 mV)³⁶. The results of several different rock standards analysed using this method have been reported previously⁶⁹, and seawater gives δ⁷Li = 31.18 ± 0.38‰ (2σ; n = 29). Of particular relevance for low-concentration carbonate samples is that aliquots of seawater (~3 µl) that were purified and then analysed at concentrations of 0.5 ng ml⁻¹ (n = 6) yield a long-term analytical uncertainty of ±0.4‰ (2σ)³⁶.

In addition to assessing whether or not our samples had been affected by silicate leaching, we used elemental ratios to evaluate if the samples from Sites 1263 and U1333 had been subject to diagenetic alteration or other forms of contamination (see Supplementary Tables 2 and 3). Both sites show little correlation between the elemental ratios (Mn/Ca, Rb/Ca, Al/Ca) and δ⁷Li, indicating that there is no effect from manganese oxyhydroxides and, as mentioned above, little to no contamination from the leaching of silicates. Although there is a short-lived peak in Mn/Ca during the MECO at Site U1333, there is no clear change in ratios during the event at Site 1263, and no identifiable trends in Rb/Ca across the MECO at either site. Likewise, there is little correlation between Li/Ca and δ⁷Li, and no clear changes in Li/Ca across the MECO. The Li/Ca values are lower than values from the PETM³⁶, but are very similar to those measured in the very pure marine carbonates at OAE1a and OAE2^29,70. Comparing our trace element analyses to recent work on assessing diagenetic affects on δ⁷Li in carbonates⁶⁶ indicates insignificant diagenetic alteration of our samples from Sites 1263 and U1333, and, combined with the similar δ⁷Li signals from globally disparate cores, this strongly suggests that these marine carbonates are recording original seawater δ⁷Li values, with a carbonate fractionation factor.

Data treatment and availability

We updated the depth-age models for Sites 1263 and U1333 to the GTS2020²⁷; this was also done for Site 702 from ref. ¹⁸, such that the proxy data (see Supplementary Information) we plot our model results against are also updated. The tiepoints used for this can be found in Supplementary Table 1. The isotopic data for all three sites can be found in Supplementary Tables 2–4, which can be found in the accompanying Excel file and is also available for download from the Figshare data repository at https://doi.org/10.5522/04/23749197. Supplementary Tables 2 and 3 also include trace element data for Sites 1263 and U1333. The isotopic data for Sites 1263 and U1333 are plotted in Figs. 1 and 2, and the Site 959 data are plotted in Supplementary Fig. 1.

To make it easier to compare the general trend of the isotopic data with our modelling results we combined the data points from Sites 1263 and U1333, and using MATLAB’s curve fitting app we applied a smoothing spline to the average data points and also to the fully propagated uncertainties of ±2σ to create a window of possible isotopic values for our timeframe of interest (41.5–39.0 Ma). For both the lithium and osmium isotope data, the piecewise polynomial is computed from p, which was set to 0.999, and x was normalized by the mean of 40.17 with a σ-value of 0.5538. The individual isotopic data for Sites 1263 and U1333 and the results from the smoothing treatment applied to the combined sites (as seen in Fig. 2) can also be found in a bundle with the model code (see the Model section below).

Model

The biogeochemical box model we use in this study (CARLIOS) is an adaptation of the CARMER model⁷¹, which was used to investigate carbon and mercury cycling during the Permo–Triassic mass extinction event, but has also been used to investigate the influence of calcium (via evaporites) on climate in the Cenozoic⁷²; thus, no processes specific to the Permo–Triassic mass extinction are included. Here, we strip out the mercury cycle and incorporate the lithium and osmium cycles (Supplementary Fig. 4). The operation of the lithium and osmium cycles is not based on the mercury cycle and, in this model, neither of these elemental cycles have a direct influence on the carbon cycle. Another introduction to this model is the inclusion of a simple silicon cycle. This enables us to investigate the influence of reverse weathering (that is, the formation of marine authigenic clays) on the carbon and lithium cycles.

The model has a crust and three ocean boxes (consisting of a low-latitude surface ocean (s), a high-latitude surface ocean (h) and a deep ocean (d), to enable the thermohaline mixing of oceanic components (for example, dissolved inorganic carbon)), as well as an atmosphere box. The carbon, silicon and osmium cycles also have atmosphere boxes, which enable air–sea gas exchange for CO₂ and the deposition of dust (for example, silicon-bearing aeolian, and osmium-bearing aeolian and cosmic) to the two surface ocean boxes. For simplicity, the schematic for the osmium cycle combines the dust fluxes and hydrothermal fluxes (Supplementary Fig. 4c), but these are treated separately within the model itself.

An important distinction from the CARMER model is that we do not assume that carbonates precipitate (that is, the burial flux of carbonates) entirely out of low-latitude shallow waters. Instead, since pelagic calcifiers had evolved by this point, much of the carbonate burial takes place in the deep ocean with some proportion of total carbonate burial taking place in warm shelf environments, following ref. ⁷³. In addition, we do not assume that the calcium carbonate solubility product has the same value for the different ocean boxes and that it was the same during the MECO as it is during the present day. Instead this is calculated throughout, following the equation utilized in ref. ⁷⁴ and using the possible oceanic concentration of magnesium at ~40 Ma (refs. ^40,75). For the lithium cycle the following are adopted: the parameterization of ref. ⁷⁶ is followed to derive the flux size and the δ⁷Li signature of rivers at time t; the present-day flux sizes and isotope ratios are used from refs. ^37,67; and the split of the sinks between marine authigenic clays and the alteration of oceanic crust is used from ref. ²⁰. In all scenarios we include in our modelling a temperature-dependent effect on the lithium isotope fractionation factor during uptake into secondary minerals (both marine and terrestrial), which, as with other warming events, will decrease the δ⁷Li of rivers (δ⁷Li_riv)^36,77. To compare our model results with our δ⁷Li_carb data, we applied a Δ⁷Li_sw-carb of 3–5‰ (the dark grey band in Fig. 3m) and 2–6‰ (the light grey band in Fig. 3m) to the average δ⁷Li_sw produced by the model, to convert it to a δ⁷Li_carb record^28,29,30,31.

We run the model forwards in time and start the model runs at 55 Ma. We run the model from 55 Ma so that the reservoirs can reach a steady state (in terms of the magnitude and isotopic values) before introducing perturbations to the system at the various time points indicated in the different scenarios detailed in the main text and the Supplementary Information. We highlight the fact that this is not a steady-state model (that is, we are not setting the derivatives in Supplementary Table 15 to zero). CARLIOS is run 5,000 times as a Monte Carlo simulation, with the model choosing values within defined ranges for certain parameters (see Supplementary Table 13). Table 1 lists the parameter changes for each scenario.

CARLIOS was constructed in MATLAB using a desktop PC with an eight-core Intel Core i7-9700 CPU @ 3.00 GHz. The model uses MATLAB’s parallel computing toolbox to utilize all eight cores to decrease the run time (typically <2 h for 5,000 model runs).

Further information on the successful model Scenario 8

We include here some additional rationale behind the decision to model a decrease in [Mg]_sw and an increase in unradiogenic osmium sourced from continental eruptions for Scenario 8. There is evidence of increasing glauconite (as seen at Site 959), smectite and palygorskite clays during the approximate MECO timeframe^46,55,78, and other research suggests that, generally, the magnitude of reverse weathering (a notable sink of magnesium) increases as temperatures increase^38,79,80. This may be a reason why [Mg]_sw has increased by ~15 mmol over the past ~40 Myr, as the deep ocean cooled and the reverse weathering uptake of magnesium reduced, although a global reorganisation of the silicon cycle may have also played a part^55,81. Combined, this information implies that the formation of marine clays during the MECO increased, taking up more magnesium and reducing its concentration in seawater. The LOWESS (locally weighted scatterplot smoothing) bootstrapping of data does hint at a decline in the Mg/Ca ratio during the MECO, as data from a ridge-flank CaCO₃ vein trends towards a low ratio, with coral data showing a large post-MECO increase⁸². The Mg/Ca data from the study used for our CO₂ proxy data (see Supplementary Information for CO₂ data)¹⁸ indicate minor decreases in Mg/Ca across the ‘warming phase’ of the MECO, with an increase at the MECO peak, although the absolute values of Mg/Ca vary substantially across the four ODP sites sampled. It may be that clay formation (marine and terrestrial) altered the Mg/Ca ratio during the warming phase (resulting in the steady accumulation of CO₂ in the atmosphere and resultant warming), but a pulse of additional CO₂ input to the atmosphere during the peak MECO led to a large amount of congruent weathering, releasing more magnesium to the oceans. Eventually this would have aided in carbonate formation, but on a shorter timescale resulted in the sharp decline in pH¹⁸. Post-MECO, as the deep ocean began to cool, the reverse weathering flux would have started to decline, leading to less CO₂ formation.

Although silicon is also trapped by clays, thermal expansion of the oceans due to the warming of the climate may have increased sand-grain dissolution⁸³. Therefore, although clays would have retained some silicon liberated during weathering, this may have been balanced, at least for a while, by enhanced quartz dissolution. It has been shown that aragonite dissolution helps to preserve calcite deposited below the CCD as it is more soluble, yet Mg-calcites are even more soluble⁸⁴. Thus, if there is a reduction in magnesium in the oceans, there is likely to be a reduction in magnesium calcites precipitating and then being transported to the seafloor. This may then reduce the protection from dissolution given to aragonite—and subsequently calcite—below the CCD, and hence the decline in the CaCO₃ content of marine sedimentary rocks deposited during the MECO¹². This may have continued until the surface layers were covered in clays, reducing the CaCO₃ exposure to acidity⁷⁴. These changes to ocean chemistry will undoubtedly have had effects on the pH and nutrient availability, and thus may have had some influence on ocean biota, but not large enough to cause extinction crises. Indeed, although there is evidence of change, there is no record of large-scale turnover, and benthic foraminifera records suggest that environmental conditions after the MECO returned to their pre-MECO state^{11,85,86,87,88}.

For the osmium cycle, simple mass-balance calculations indicate that the ¹⁸⁷Os/¹⁸⁸Os of rivers must have been considerably less radiogenic than they are today (at 1.2–1.5)³⁴ to achieve the ¹⁸⁷Os/¹⁸⁸Os_initial signature of the pre-MECO seawater¹⁶. Remnants of previous large igneous provinces (LIPs) were in or transited through the humid zone before the MECO⁸⁹ (Supplementary Fig. 5) and the ¹⁸⁷Os/¹⁸⁸Os_initial of seawater post-emplacement of these LIPs shows distinctly unradiogenic values (for example, ref. ⁹⁰). It is possible that osmium was leached from rocks during weathering of the LIPs and was then subsequently trapped by clay-associated organic matter. This osmium was probably then released during clay dissolution in the humid zone, and thus we surmise that the weathering of these LIPs supplied unradiogenic osmium (compared with the present day), even when accounting for increasing radiogenicity due to the time elapsed since emplacement. This is supported by the relatively stable ⁸⁷Sr/⁸⁶Sr values before the MECO as well⁹¹. For the MECO itself, Scenarios 1–6 suggest that enhanced weathering, either due to greater $p_{{\mathrm{CO}}_{2}}$ levels or erosion, is not sufficient to replicate the excursion seen in the data, and thus in Scenarios 7 and 8 we invoked a moderate increase (~75% of the present-day combined hydrothermal fluxes) in the delivery of osmium to the ocean with a mantle-like ¹⁸⁷Os/¹⁸⁸Os, as data from various LIPs suggest that eruptive material was very unradiogenic⁹⁰. However, the trend towards a more unradiogenic seawater signal could instead have been the result of a diminished silicate weathering flux¹⁶ or a mixture of the two processes. As we do not prescribe any short-lived changes within the MECO timeframe, our model misses the small mid-MECO increase in ¹⁸⁷Os/¹⁸⁸Os_initial, which may have occurred due to greater weathering of more radiogenic materials or a temporary decrease in the amount of mafic-associated volcanism.