Enhanced clay formation key in sustaining the Middle Eocene Climatic Optimum

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 CO2 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 CO2 over this timescale. Here we investigate silicate weathering during this climate warming event by measuring lithium isotope ratios (reported as δ7Li), which are a tracer for silicate weathering processes, from a suite of open-ocean carbonate-rich sediments. We find a positive δ7Li 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.


Site 959
The sediments of Site 959 for the interval explored in this study are principally composed of various forms of porcellanite, with fluctuations in carbonate content as the dominant lithology cycles through the following types: porcellanite foram calcareous chalk, clayey porcellanite, porcellanite micrite chalk, micritic porcellanite, porcellanite with micrite, and porcellanite.For depths associated with the MECO there is also evidence that the sediments are relatively organic-rich, with some dissemination of glauconite and pyrite 1 .The 187 Os/ 188 Osinitial values at Site 959 are marginally more radiogenic than the other two sites and the isotope excursion appears to begin earlier, but these differences are not substantial enough to suggest that Os dynamics operated differently at Site 959 compared to Sites 1263 and U1333.Site 959 also exhibits a positive LIE (Fig. S1), moving from a low of ~8‰ prior to the MECO, to a high of ~15‰ during the MECO.The excursion begins before the start of the MECO as defined by TEX86 data 2 and has a longer duration.However, there has been an assumption that changes to the TEX86 data should be co-eval with changes to the δ 18 O record, but this may not necessarily be true and non-thermal factors may have influenced the record 3,4 .As such Site 959 could have a MECO duration similar to the other two sites, but the age constraints for this site are very uncertain, thus we do not assign specific ages to the data from this site.More importantly, our δ 7 Li data for this site show a much lower starting value, pre-MECO, to the other two sites, with an excursion that is approximately double in size (see compared to Fig. 1).As the sediments from all three sites were deposited in open, deep ocean settings at ~40 Ma, a possible reason for the difference in δ 7 Li data at Site 959 may be that the porcelaneous nature of the sediments has induced a lithium isotope fractionation factor very different from that typical of bulk carbonates, charting at the extreme end of that seen in biogenic aragonites 5 .Hence, we focused the bulk of our main text discussion on Sites 1263 and U1333.

Temperature and CO 2 proxies
For ground-truthing our model outputs, it is crucial to compare them to proxy data.Alongside the δ 7 Licarb and 187 Os/ 188 Osinitial data from this study and van der Ploeg et al. 6 , we also compare our model results against pCO2 and temperature proxy data.Atmospheric CO2 proxy data for the MECO comes almost exclusively from one study 7 , with only one data point from stomatal indices and one data point from alkenones, from just before and just after the MECO, respectively, available 8 .The pCO2 data points (see Fig. S2) are based on δ 11 B from planktic foraminifera from ODP Sites 865, 1260, 1263 and 702.The sites were aligned and ages reported on the 2012 GTS 7,9 but we have updated the ages to the 2020 GTS 10 .The authors converted δ 11 B to a pH and then, using their best estimate of δ 11 Bsw = 38.5 -38.9‰, converted pH to pCO2 values, assuming either a constant ocean saturation state (Ωcalcite) or a constant ocean alkalinity (ALK), which represents, approximately, the two pCO2 endmembers 7 .
Although there are uncertainty ranges for the data (see their Fig.7) 7 , we do not include these in our figures, as our aim is to approximately match the timing and general magnitude inferred from the δ 11 B data.

Figure S2 | Atmospheric CO2 concentrations during the MECO. The concentrations are from
Henehan et al. 7 , where they calculated CO2 by using their best estimate of δ 11 Bsw obtained from foraminifera and assumed either a constant calcite saturation (black dots) or a constant alkalinity (blue squares).The grey line is the result of using MATLAB's curve fitting app to fit a smoothing spline to the CO2 estimations from constant calcite saturation and is plotted here to show the general trend of the changes to CO2.
For the temperature proxy data, we use the δ 18 O data from benthic foraminifera (ODP Sites 865 and 1260) and from bulk sediments (ODP Sites 1263 and 702), which Henehan et al. 7 collated and tied to the 2012 GTS 9 , and as with the δ 11 B data, we have updated to the 2020 GTS 10 (see Fig. S3).Hansen et al. 11 derived a set of equations that can be used to convert δ 18 O data to temperature, with different equations to be employed when the Earth was ice-free (and therefore δ 18 O values are <1.75‰) and when the Earth had ice sheets and thus changes in δ 18 O reflect changes to the ice mass as well as temperature (δ 18 O data is >1.75‰).As all the δ 18 O data here 7 is <1.75‰, we convert δ 18 O to surface temperatures following Hansen et al. 11 using a surface temperature of 22.07℃ at 40.884 Ma: (1) We then used MATLAB's curve fitting app to fit a smoothing spline to our data (see Fig. 3b).Although Henehan et al. 7 generated δ 18 O data from their planktic foraminifera samples, we do not use these in this study because even in the modern day it is difficult to disentangle seasonal changes to shell δ 18 O from depth migration effects 12 , thus making it more difficult to convert δ 18 O to surface temperatures.

Figure S3 | Global deep ocean δ 18 O data across the MECO and the resultant estimated surface
temperature.In a) δ 18 O data is from benthic foraminifera at ODP Sites 865 and 1260, and from bulk sediments at ODP Sites 1263 and 702, with ages tied to the 2020 Geologic Time Scale 7,10 .See text for more details.In b) δ 18 O data is converted to a global surface temperature using the equations from Hansen et al. 11 , and a smoothing spline fit to the data (green line) is added.The 'start', 'peak' and 'post' MECO times are based on the δ 18 O data from Site 702 7 .

Schematic
The schematic for the CARLIOS model can be seen in Fig. S4.Much of the information about the model can be found in the Methods section of the main text, however we also discuss some of the key parameters below.

Long-term tectonic degassing
The long-term tectonic degassing rate (fG) is a significant variable in long-term box models [13][14][15] , controlling the return rate of carbon (and sulfur etc.) species from the crust to the ocean/atmosphere via both degassing and metamorphism.This parameter, based on seafloor spreading rates, has also been used in simple box modelling to control the rate of hydrothermal input/output of other elements, such as lithium or strontium 14,16 .In terms of CO2 input to the exogenic system (i.e.oceans and atmosphere), basing fG on seafloor spreading rates likely does not fully encapsulate the effects of continental rifting and continental arc volcanism, both of which can be significant sources of CO2 to the atmosphere, especially if lithospheric or crustal carbonates are remobilised as a result [17][18][19] .Indeed, incorporation of continental rift lengths into fG has been shown to somewhat close the gap between pCO2 estimates from geochemical proxies and those generated by models such as GEOCARBSULF 20,21 for the Cenozoic.
Such detailed deep carbon cycling work is beyond the scope of this study.Instead, fG prior to the MECO is estimated to be 1.25 -1.43x the present-day level based on seafloor spreading rates 13,22 but tectonic degassing may have been up to 3x the present-day level when rift lengths are taken into consideration 23 .
The use of the upper value of the range (when taking into account the additional atmospheric CO2 (see below)) allows our model to produce a reasonable fit to pre-MECO CO2 and temperature estimates, and the input of elements into the ocean sourced from seafloor spreading also results in a general match to our isotopic proxy data.However, because of the uncertainty range on seafloor spreading rates, and tectonic degassing in general, we incorporate fG into our Monte Carlo simulations with bounds of 1.25 -1.65x the present-day level.We then model different timings and amounts of additional CO2 input directly to the atmosphere.

Additional atmospheric CO2 injection estimates
Prior to and during the MECO, volcanic CO2 degassing from the India -Asia collision (in Tibet -the Himalayas), may have been over an order of magnitude higher than is presently released in this area 18 .
It has also been hypothesised that a flare up of continental arc volcanism in the Neotethys suture zone (in what is present day Iran and Azerbaijan) could have been a trigger for the MECO 19 .Elsewhere along the Neotethys suture zone, the volcanic arc extends westwards (into what is present day Armenia, Georgia and Turkey) 24 , and there is evidence for regional uplift and a switch from submarine to subaerial volcanism (which is important for increasing CO2 degassing to the atmosphere without supplying isotopically light δ 7 Li to the ocean) from the Lutetian to Bartonian Ages 25 .There is also a minor increase in the areal extent of continental large igneous province (LIP) eruptions from 41 to 40 Ma 26 , which may be linked to some initial eruption of flood basalts in Southern Ethiopia/Northern Kenya prior to the main eruptive phase of the Afar-Arabian LIP 27,28 .While some of the CO2 degassed from these locations may have been accounted for in fG, we assume that overall fG is an underestimate of the total amount of CO2 delivered to the atmosphere from tectonic activity and thus model an additional input of CO2.For each scenario (ScX, where X indicates the scenario number) we assume that there is a baseline of additional CO2 which is steadily degassing from Tibet 18 , and we incorporate the uncertainty regarding the magnitude of this flux into our Monte Carlo model simulations.Then for Sc1 and Sc6-8 we model a further injection of CO2 into the atmosphere during the MECO sourced, primarily, from Iran 19 , but also possibly from other volcanic arcs and continental LIPs as described above.Unlike for Iran, however, we currently have no estimations as to how much CO2 may have degassed from other volcanic arcs during this time.

Long-term uplift
The long-term uplift of sediments and their subsequent erosion (fR) is another key variable in biogeochemical box models.Based on sediment deposition reconstructions, with a smoothing of 100 Myrs, during the mid-Eocene the long-term uplift value is estimated to be approximately 0.5x the present-day 29,30 .Increasing uplift is a common feature in inter-eruptive or pre-eruptive stages for volcanoes 31 , while eruptions themselves increase erosion, either via lavas or lahars 32 .Additionally, due to increasing atmospheric CO2 levels (and thus temperature), rainfall patterns, magnitudes and intensity can all change to produce large variations in erosion on geologically short timescales that may not be evidenced in the long-term record 33,34 .Thus, we might expect an increase in erosion due to an increase in volcanic activity, alongside complex topographical changes in Tibet producing some uplift 35 and for example, evidence for intermittent but substantial uplift in the Isle of Wight 36 .Also, at 40 Ma some of the remaining LIP materials were in areas of higher altitude (Fig. S5), such as northern South America, the Alps, China and Mongolia, and parts of Tibet, meaning that some volume of basaltic rock would have been uplifted to zones more susceptible to high rates of erosion and/or weathering.As such, we model increases in fR in Sc2 and Sc6-8.

Figure S5 | Vestigial large igneous provinces during the MECO. The paleogeographic distribution of
LIP-associated sediments (in yellow) remaining from previous emplacements, at ~40 Ma.This map was made in GPlates with the base paleogeography from Scotese 37 and the LIP distribution from ref. 26 and refs.therein.
In Sc8 (our preferred scenario, see main text) and Sc7 (see below) we set the minimum value for the pre-MECO to 0.5x and the maximum value during the MECO to 1x present day uplift.These values are chosen as the sediment reconstruction and smoothing of the data 29,30 , as well as a separate polynomial fit to terrigenous sediments 38 indicate that the long-term fR has not exceeded ~1.1x the present-day rate throughout the Phanerozoic, albeit rapid erosion events are likely missed out in the data.Although it may appear as if we are utilising a small increase in erosion, if the pre-anthropogenic erosion flux is ~2x10 10 tonnes yr -1 39-41 , then we are modelling an increase of 1x10 10 tonnes yr -1 at the peak of the MECO.We do, however, model two scenarios (Sc2 and Sc6) whereby fR reaches a peak of 3x present-day to test what a large erosion event might do to the Earth system.

The erosion to silicate weathering ratio, δ 7 Liriv and FLiriv
In their modelling, Caves-Rugenstein et al. 42 calculate the lithium isotopic signature of rivers (δ 7 Liriv) and the riverine flux of Li (FLiriv) dynamically.They estimate a weathering intensity (W/D, where W is weathering and D is denudation, which is the total of erosion (E) + W) by first converting the model silicate weathering flux from moles yr -1 into tonnes yr -1 .Using: this flux; the normalised uplift (fR); and a global erosion to weathering ratio (E:W), they produce an erosion flux in tonnes yr -1 .They then use their W/D values to generate a value for δ 7 Liriv.The δ 7 Liriv can then be used to derive a value for the fraction of lithium partitioned from the bedrock into the dissolved load via a Rayleigh distillation function, which is then used to calculate FLiriv.An uncertainty in this multistep approach is the present day global silicate weathering flux 43 which is used to convert between Mol yr -1 to tonnes yr -1 .However, the magnitude of silicate weathering is typically outweighed by the erosion flux when calculating denudation rates, such that at the present day a 28% reduction in the silicate weathering rate only decreases total denudation 41,43,44 from 2.25x10 10 tonnes yr -1 to 2.04x10 10 tonnes yr -1 .A more significant unknown is the global erosion to weathering ratio (E:W), which is only used in the Li and not the C cycle calculations.It is assumed that E:W at 16 Ma ranges between 3-40 42,45 , but this may have been different at ~40 Ma, as fR is lower during the Eocene than the Neogene 29,30 .
As our data show a rise in the δ 7 Lisw it is likely that there was some increase in the riverine δ 7 Li (some of the seawater excursion may have also been sourced from changes to the Li sinks -e.g.see Extended Data Fig. 3), but without modelling we do not know whether the Earth moved from a kinetically-limited, or a supply-limited version of congruent weathering, to incongruent weathering [46][47][48] .We used the parameterisation of Caves-Rugenstein et al. 42 to calculate the δ 7 Liriv and FLiriv dynamically, but for Sc1-6 we kept E:W fixed throughout the model runs at a value of 5 which we found generated a pre-MECO δ 7 Liriv signature that leads to the model δ 7 Licarb record being consistent with our sample data (see Extended Data Figs. 1 -7m).When combined with large increases in fR (such as in Sc2 and Sc6), this produces a reasonable fit to the δ 7 Licarb data.For Sc7 and Sc8 we included the E:W as a variable in the Monte Carlo simulation, setting the boundaries to incorporate the full Neogene range (3-40) 42,45 and used a more modest increase in fR which produces a very good fit to the δ 7 Licarb data.

Additional information on the failed scenarios -Extended data figures 1 -7
In Extended Data Figs. 1 -7 we plot the key parameter changes plus model results for Sc1-7: the failed scenarios.Some of the scenario hypotheses (e.g.Sc1 and Sc4) were based upon previous work 49 but parameter values for, as an example, the input of additional CO2 to the atmosphere (Sc1) or a near-total shutdown in marine organic carbon (Corg) burial were chosen to try and match our pCO2 model outputs to the proxy data.Table 1 in the main text illustrates why these model runs are considered failures (e.g. the lysocline does not shoal, or the model does not reproduce the δ 7 Licarb data etc.).In Sc3 we apply a scaling to the reverse weathering flux to artificially force an increase in its magnitude, while holding fR and the additional CO2 input at their pre-MECO values.This scenario was conducted to test the effect of a pulse of reverse weathering (and thus a large increase of marine authigenic clay (MAAC) formation) on the modelled estimate for δ 7 Lisw (and thus δ 7 Licarb).However, we do not fix the riverine flux or its δ 7 Li signature in this scenario, instead it varies in response to the CO2 released to the atmosphere from reverse weathering reactions.Likewise, in scenarios such as Sc1 or Sc2 where we change the additional CO2 input or the fR value, which affects terrestrial weathering and erosion, we do not hold the lithium sink fluxes (MAAC formation and the alteration of oceanic crust) at a constant size.Rather, these fluxes vary in response to the amount of silicon in the ocean.Thus, in Fig. 3o in the main text, and Extended Data Figs. 1 -7o, the rivers versus sinks panel reflects changes to both the riverine input and the combined outputs (both flux size and isotopic value) in each modelled scenario, but showcases whether the rivers or the combined sinks hold the greater control over the δ 7 Lisw record in our model (<1 for sinks, >1 for rivers).

Calculating changes to Mg and Ca retention in soils
Our modelling, based on the W/D results (Fig. 4), suggests the Earth transitioned (pre-MECO) from a global regime with a mix of secondary mineral formation and dissolution, with dissolution slightly outweighing formation, to one more dominated by secondary mineral formation (even if only slightly) during the MECO.It is possible that the early Cenozoic warmth 50 coupled with the eruption of large igneous provinces (e.g. the North Atlantic Igneous Province) 51 , which generated lots of highly weatherable basalt, resulted in high rates of weathering and erosion but eventually in widespread secondary mineral (e.g.clay) formation well before the MECO.Over time, as clays became thicker and much of the basalt from the LIPs was eroded, the Earth slowly transitioned, globally, to a more supply limited regime, where clays started to dissolve.As noted in the main text, there will have been spatial heterogeneity, with some areas experiencing high erosion and low chemical weathering (and vice versa).This ties in with the δ 7 Li record of seawater in the Paleogene 52 where δ 7 Lisw rises by a few ‰ after the Early Eocene Climatic Optimum (~50 Ma) to ~44 Ma whereafter the limited available data suggests δ 7 Lisw starts to decline and is backed up by our low δ 7 Li pre-MECO data.
The present day global δ 7 Liriv average is 23‰ 48 and the W/D, based on estimates of total rock weathering 44 of 2.13x10 9 t yr -1 and total erosion 41 of 2x10 10 t yr -1 , is 0.096 (see main text Fig. 4a), although the weathering estimate does not include contributions from oceanic island basalts.It has been proposed that the silicate weathering flux may have been overestimated 43 , however, at most, this reduces W/D to 0.090.Meanwhile, a modelled estimate for W/D for the pre-anthropogenic present is ~0.063, but the δ 7 Liriv is higher at ~28‰ 42 .Our modelling suggests that during the MECO, the global δ 7 Liriv average was ~23-24‰ with a W/D of ~0.065-0.085,which is very similar to the present day.We attempt to work out what the change in weathering regime may have meant for the amount of Ca and Mg cations retained on land, and whether such an increase is feasible.
The general trend of soil formation rates appear to follow changes in climate over time, with higher rates occurring during periods of higher temperatures 53 .Assuming that surface temperatures have declined by 6℃ over the last 50 Myrs [54][55][56] , while the soil formation rate 53 has declined from 2.2×10 10 to 1.3×10 10 tonnes yr −1 , we calculate an increase in soil formation at the MECO peak to 2.280x10 10 tonnes yr -1 compared to the pre-MECO background formation rate (~1.830×10 10 tonnes yr −1 ), assuming a global average surface temperature increase of 3℃ 57 .
If the Earth's land surface is approximately 148,940,000 km 2 and assuming that approximately 94% of this area could host soil formation 58 , with a soil depth of 2 m, this gives an approximate value of 280,007,200 km 3 or 2.800x10 17 m 3 soil volume, split as 4.200x10 16 m 3 topsoil and 2.380x10 17 m 3 subsoil.A tonne of a standard, moist soil, equals approximately 0.670 m 3 thus at the present day, there is 6.269x10 16 tonnes of topsoil and 3.552x10 17 tonnes of subsoil (4.179x10 17 tonnes in total).
To work out how much of the soil is composed of clay-size minerals we use the clay fraction (% of the soil that is composed of clays) per unit area data (2° grid resolution) for the present day Earth 59 .Any values that were less than 0 (indicative of grid cells with no data, such as the ocean) were converted to NaNs in MATLAB.Out of 16,200 cells for the topsoil, 5,117 of these were NaNs (the same was true for the subsoil dataset).The clay fraction data for the topsoil for all grid cells was summed together (see Table S7) and then divided by 11,083 (the number of cells with data, i.e. non-NaNs) to give a global average clay fraction for the topsoil.The same was done with the subsoil dataset.As the average clay fraction of topsoil and subsoil is 9.42% and 10.86% respectively, using our estimates of soil tonnage, we can thus surmise that clay-size minerals make up 5.905x10 15 tonnes and 3.858x10 16 tonnes respectively (4.448x10 16 tonnes in total).
However, not all clay-size minerals contain Ca and/or Mg, so we work out the fraction of total clay that is kaolinite, smectite, vermiculite or illite, which are the most abundant clay-size minerals in soils and are Ca or Mg containing 59 .Again, any areas of the Earth with no data are converted to NaNs and omitted.For each mineral type the topsoil and subsoil are summed together, to give a global value (see Table S8).The kaolinite, smectite, vermiculite and illite data are also added together.us, using our tonnage of clays in soils (4.448x10 16 tonnes) and the percentage of Ca and Mg in clays, we can calculate the amount of these cations in total in soils at Present: 6.103x10 14 tonnes.
If there was an increase from 0.5x to 1x present day erosion across the MECO (corresponding to 1x10 10 tonnes yr -1 to 2x10 10 tonnes yr -1 ), as per our favoured Sc8 (see main text), then using our soil formation rates across the MECO (an increase from ~1.830×10 10 tonnes yr −1 to 2.280x10 10 tonnes yr -1 ), this would indicate a change of net soil formation from 0.17x10 10 tonnes yr -1 to 0.28x10 10 tonnes yr -1 .Assuming that the percentages of the various clay-size minerals, and thus the amount of Ca and Mg, in soils has not changed considerably over time (and we acknowledge that this is an area of considerable uncertainty), we calculate the net addition of Ca and Mg to soils across the MECO would increase from 2.332x10 7 tonnes yr −1 to 3.841x10 7 tonnes yr −1 .
If the MECO lasts 400 kyrs, as per our modelling, using linear interpolation between the two net addition of Ca and Mg values, and assuming they are added via an equally distributed 'witch's hat' shape, the total addition of carbonate forming cations to the soils is 1.235x10 13 tonnes.If the pre-MECO soil Ca and Mg reservoir was similar to today, then this addition would make the Ca and Mg reservoir equal to 6.227x10 14 tonnes at the end of the MECO, which represents a ~2% increase in the amount of Ca and Mg retained in soils.Of course, the pre-MECO Ca and Mg soil reservoir may have been much smaller than at Present, which means a greater percentage increase.Alternatively, because the W/D calculated by the model at the MECO peak is approximately equal to that at the present-day (see Fig. 4a), we could deduce that prior to the MECO the terrestrial reservoir of Ca and Mg was equal to 5.98x10 14 tonnes (6.103x10 14 -1.235x10 14 tonnes).For comparison, the amount of Mg and Ca in seawater is approximately as follows: Mg: 1.29x10 21 l * 1.272 g/l = 1.64x10 21g, which is 1.64x10 15 t Ca: 1.29x10 21 l * 0.411 g/l = 5.3x10 20 g, which is 5.3x10 14 t Summed together, the total amount is 2.17x10 15 t which suggests that our estimates for Ca and Mg retained in soils is feasible in the context of the present-day Earth system and also highlights the potential power that clay retention of cations may have on ocean chemistry and climate.

Other possible MECO triggers
Although in this study we have invoked an increase in tectonic activity and associated volcanism to explain the MECO, there are other potential triggers for the MECO that have been suggested.One of these is the influence of orbital changes, as it has been found that the MECO occurs coincidentally with the minima of ~400 kyr and 2.4 Myr eccentricity cycles 61 and the timing of these correlate to variations in the δ 13 C and δ 18 O records 62 .It is uncertain as to whether eccentricity changes alone could have triggered the MECO (perhaps via changing ocean and hydrological circulation patterns, affecting weathering and primary productivity rates), or whether they perhaps enhanced the effects of some other trigger, such as enhanced volcanism.
Another possibility is the impact of bolide(s) upon the Earth's surface.As an extraterrestrial body, a bolide would likely contain unradiogenic Os 63 and depending on the lithology of the impact location, could expose silicate minerals to the surface and/or vaporise carbonate rocks.In Eurasia, several bolides potentially of MECO age, have been documented 64 , with the Beenchime-Salaaty impact crater exposing Cambrian age carbonate rocks, the Logancha uncovering Triassic basalts, and the Sunak exposing primarily Devonian rhyolites.However, the age constraints on these impact events remains poor.
Additionally, bolides have low δ 7 Li values 65 which could potentially lower δ 7 Lisw dependent on whether the impact site was sufficiently far enough away from the coastline to allow for clay formation and isotopic fractionation.Thus, in light of all available data we favour changes in tectonic activity to explain the MECO.Conversion from Sv to m 3 yr -1 for thermohaline circulation thconv 3.15x10 13 Dal Corso et al. 66  This study (all sources combined minus reverse weathering) *N.B.Osmium fluxes are taken to be the average of the range from Lu et al. 63 .

Full model equations
Convert reservoirs at present day/model initiation from moles to mol/m 3 Convert reservoirs (at time = t) from moles to mol/m 3 ‫ܿ݊ܿܺ‬ = ܺ ‫ݎ݁ݐܽݓ‬ Thermohaline circulation (m 3 yr -1 ) Low latitude surface to high latitude surface ocean DIC transfer Deep ocean temperature ‫݉݁ݐ‬ ௗ = ݉ܽ‫ݔ‬ሺ275.5+ ‫ܶܵܣܩ‪ሺ‬‬ − ‫݉݁ݐ‬ ሻ, 271ሻ Change in GAST from present day ‫݉݁ݐܽ‬ = ‫ܶܵܣܩ‬ − ‫݉݁ݐ‬ Change in low latitude surface ocean temperature from present Change in deep ocean temperature from present ‫݉݁ݐ݀‬ = ‫݉݁ݐ‬ ௗ − 275.5 Temperature dependency for basalt weathering Temperature dependency for granite weathering Temperature dependency for carbonate weathering Temperature dependency for seafloor weathering CO2 dependency for all weathering Combined CO2 and temperature dependency for basalt weathering

Figure S1 |
Figure S1 | Isotopic data for ODP Site 959.δ 7 Licarb values are in blue (this study) with the filled circles representing the mean of each sample (n = 3) and the error bars represent the ± 2SD precision, 187 Os/ 188 Osinitial values are in gold 6 .

Figure S4 |
Figure S4 | Schematic of the CARLIOS model.Panel a) The carbon cycle, b) the lithium cycle, c) the osmium cycle.The ocean is split into three boxes: a low latitude surface ocean (s), a high latitude surface ocean (h) and a deep ocean box (d), and a thermohaline circulation (fcirc) mixes componentsbetween the boxes.We assume lithium has no atmospheric component in this model.In a) an air-sea gas exchange mixes carbon between the atmosphere and surface ocean boxes, degassing supplies inorganic (fmc) and organic carbon (fmg) to the atmosphere, as does the oxidation of organic carbon (fwg), while the weathering of silicates (fws) and carbonates (fwc) transfers carbon from the atmosphere and land to the ocean.The burial of carbonates on the shelf (fbc_shelf) and in the deep ocean (fbc_deep), seafloor weathering (fsfw) and terrestrial and marine organic carbon (fbg_land, fbg_sea) returns carbon to the sediments.In b) lithium is supplied to the oceans by rivers (friv) and hydrothermal fluids (fhyd) and is returned to the sediments via marine authigenic clay formation (fmaac) and alteration of oceanic basalt (faoc).In c) the sources of Os include the dust flux (fdust), which includes both cosmic and aeolian dust, the riverine flux (friv) which includes that from rivers and groundwaters, the hydrothermal flux (fhyd) including both low-temperature and high-temperature hydrothermal activity, and the terrestrial weathering of juvenile crust (fman).The sinks of Os are combined (fsinks).In d) silicate weathering (fws), aeolian dust (fdust), seafloor weathering (fsfw) and hydrothermal fluids (fhyd) deliver silicon to the ocean and the two sinks are via reverse weathering (frw) and SiO2 formation (amorphous and opal, fopal).

Table S8 | Fraction of Ca and or Mg in the major clay-size minerals.
* As montmorillonite, ^From Tardy et al.

Table S9 | Model constants. Parameter Algebraic Representation Value Source
63B.187Os/ 188 Os of aeolian dust is the average of the range in Lu et al.63. *