Constraints on Earth System Functioning at the Paleocene‐Eocene Thermal Maximum From the Marine Silicon Cycle

The Paleocene‐Eocene Thermal Maximum (PETM, ca. 56 Ma) is marked by a negative carbon isotope excursion (CIE) and increased global temperatures. The CIE is thought to result from the release of 13C‐depleted carbon, although the source(s) of carbon and triggers for its release, its rate of release, and the mechanisms by which the Earth system recovered are all debated. Many of the proposed mechanisms for the onset and recovery phases of the PETM make testable predictions about the marine silica cycle, making silicon isotope records a promising tool to address open questions about the PETM. We analyzed silicon isotope ratios (δ30Si) in radiolarian tests and sponge spicules from the Western North Atlantic (ODP Site 1051) across the PETM. Radiolarian δ30Si decreases by 0.6‰ from a background of 1‰ coeval with the CIE, while sponge δ30Si remains consistent at 0.2‰. Using a box model to test the Si cycle response to various scenarios, we find the data are best explained by a weak silicate weathering feedback, implying the recovery was mostly driven by nondiatom organic carbon burial, the other major long‐term carbon sink. We find no resolvable evidence for a volcanic trigger for carbon release, or for a change in regional oceanography. Better understanding of radiolarian Si isotope fractionation and more Si isotope records spanning the PETM are needed to confirm the global validity of these conclusions, but they highlight how the coupling between the silica and carbon cycles can be exploited to yield insight into the functioning of the Earth system.


Introduction
The early Paleogene (65 to 50 Ma) is characterized by atmospheric and ocean temperatures up to 15°C warmer than the preindustrial period (e.g., Cramwinckel et al., 2018;Evans et al., 2018), ice-free poles, and atmospheric CO 2 concentrations (pCO 2 ) around 1,000 ppmv (Anagnostou et al., 2016). Hyperthermalstransient periods of extreme warmth occurring on 10 5 year timescales-are superimposed on a long-term warming trend (Lourens et al., 2005;Zachos et al., 2001). The most notable, the Paleocene-Eocene Thermal Maximum (PETM), is particularly interesting due to its potential parallels with modern anthropogenic climate change.
The PETM is characterized by a negative carbon isotope excursion (CIE) of 2-6‰ in both marine and terrestrial archives (McInerney & Wing, 2011). The existence of this globally manifest CIE requires a significant amount-thousands of Gt-of 13 C-depleted carbon to have been added to the exogenic carbon pool. Mass balance dictates that the more 13 C depleted the carbon source, the less carbon is required to produce a given CIE (Dickens, 2001;Pagani et al., 2006). An internally consistent understanding of the PETM requires that the amount and δ 13 C of carbon is compatible with terrestrial and marine observations of the CIE, as well as the magnitudes and rates of warming, changes in ocean pH, and any carbonate compensation depth (CCD) change. Additional unknowns include the rate and duration of carbon addition and the site(s) of its release.
Given the difficulties in constructing precise, millennial-scale age models for early-Cenozoic sediments, and in interpreting proxy records, the source of this carbon has been variously attributed to highly 13 C-depleted, bacterially generated methane (δ 13 C~−60‰) from the thermal dissociation of seafloor gas hydrates (Dickens et al., 1997;Katz et al., 2001), the thermogenesis of methane via sill emplacement in organic carbon rich marine sediments (δ 13 C = −35 to −50‰; Svensen et al., 2004), volcanic degassing of magmatic CO 2 (δ 13 C = −6‰; Gutjahr et al., 2017), or to the decomposition of organic carbon from thawing permafrost (δ 13 C = −27 to −22‰; DeConto et al., 2012;Panchuk et al., 2008), among others. Overall, the source and amount of carbon, the mechanism(s) and/or thresholds that triggered its release, the timescales involved, and the processes involved in the return to preperturbation conditions are still heavily debated almost 30 years after the first identification of the PETM CIE-and is further complicated by the discovery of similar smaller events through the early Eocene (e.g., Lourens et al., 2005).
Conventionally, four phases are recognized. The first consists of pre-CIE conditions, during which δ 13 C values very gradually decrease and reconstructed temperatures increase according to the long-term trend (Zachos et al., 2001). The second-the CIE onset-is taken to be the period between the last sample with pre-CIE values and the peak of the excursion. Estimates of the onset duration include 8-23 kyr from continental sedimentary sections (Aziz et al., 2008;Magioncalda et al., 2004), circa 10 kyr from marine carbonates , and it has recently been modeled as likely less than 5 kyr (Kirtland Turner et al., 2017;Kirtland Turner & Ridgwell, 2016;Zeebe et al., 2016). It is noteworthy that at many locations, temperature appears to start increasing 3-6 kyr before the CIE onset (Frieling et al., 2019;Littler et al., 2014;Sluijs et al., 2007). A CIE "body" constitutes the third phase, with low and stable δ 13 C and temperatures (Bowen, 2013;McInerney & Wing, 2011). The CIE body lasts for about 110 kyr and is followed by the final recovery phase during which δ 13 C and temperatures gradually return to pre-CIE values over~80 kyr (McInerney & Wing, 2011;Murphy et al., 2010).
If the source of the 13 C-depleted carbon was internal to the Earth's surficial carbon cycle (i.e., not predominantly of mantle origin), then some aspect of the Earth system must have crossed a threshold to trigger the release, whether from methane hydrates, sedimentary organic matter, permafrost, or elsewhere. Many triggers have been proposed, including a comet impact (Kent et al., 2003), orbital forcing (Littler et al., 2014;Lourens et al., 2005;Zachos et al., 2010), and volcanism, in particular associated with the formation of the North Atlantic Large Igneous Province (LIP; Dickson et al., 2015;Larsen et al., 2003;Svensen et al., 2004;Wieczorek et al., 2013). The return of climate and by inference the size of the exogenic carbon pool to their pre-excursion states must have been achieved by increases in one or both of the two key sinks in the long-term carbon cycle, namely, (i) the net burial of organic carbon (Bowen & Zachos, 2010) and (ii) the conversion of atmospheric CO 2 to carbonate alkalinity via the chemical weathering of silicate minerals, and subsequent carbonate mineral burial. The relative importance of these two processes is uncertain (Bowen & Zachos, 2010).
The global silica and carbon cycles are intimately coupled across a range of space and timescales in part because they are both involved in the weathering of silicate rocks, both cycle through the biosphere, and both affect the ocean's alkalinity. Thus, many of the processes invoked to explain aspects of the PETM onset, body, or recovery also affect the silica cycle. For example, increased silicate weathering should increase the flux of silica to and burial in the ocean (Penman, 2016), at least over timescales greater than the residence time of dissolved silica in seawater (ca. 10 kyr). Overall, potential changes to the marine silica cycle can be grouped into two categories: a change in the nature and magnitude of dissolved silica (dSi) entering the ocean or a change in some internal aspect of silica cycling within the ocean. Both types of change may be manifest as a change in seawater dSi concentrations and/or 30 Si/ 28 Si ratios (expressed as δ 30 Si in per mille deviation from the standard NBS28). This means that archives of seawater δ 30 Si have previously unexploited potential to provide answers to questions about the relative importance of different processes during the PETM. Such an archive exists in the δ 30 Si of siliceous organisms, particularly diatoms, radiolarians, and siliceous sponges.
It is well established that siliceous microfossils can be used to investigate various aspects of the marine silica cycle (Abelmann et al., 2015;Fontorbe et al., 2016;Fontorbe et al., 2017;Hendry & Brzezinski, 2014;Hendry & Robinson, 2012). Silicon isotope ratios in each group provide information about the silica cycle in the past. Diatom δ 30 Si, for example, is conventionally used as a euphotic zone Si-utilization proxy. Radiolarians live mostly in the upper few 100 m (Botlovskoy 2017; Suzuki & Not, 2015) but do not bloom, so their δ 30 Si is considered more of a passive tracer of water mass δ 30 Si. Sponge spicules have the useful characteristic that their Si isotope fractionation is related to ambient dSi concentrations, making them a proxy for bottom water dSi concentrations (hereafter [Si]) (Hendry et al., 2019).
Here, we present δ 30 Si values from radiolarian tests and sponge spicules spanning the PETM at Blake Nose in the Western Atlantic (ODP Leg 191B Site 1051). Using a simple box model, we investigate the sensitivity of seawater Si isotopes to changes in sources to, and cycling of, Si in the ocean, and thus their ability to support or challenge interpretations based on existing geochemical, paleontological, and sedimentological data.
This period encompassing the PETM in this core has previously published cyclostratigraphic age models based on core-scanning XRF Fe data tuned to orbital cycles (Norris & Röhl, 1999) and 3 He concentrations, assuming a constant extraterrestrial 3 He flux (Farley & Eltgroth, 2003). These models broadly agree and at face value suggest a total PETM duration of 150 ± 20 and 120 kyr, respectively. However, subsequent work suggested Blake Nose may contain a substantial hiatus (of about 55 kyr) around the onset of the CIE (Farley & Eltgroth, 2003) when compared to Site 690 (Weddel Sea, Southern Ocean). This may have been associated with slumping induced by the release of hydrates (Katz et al., 2001). The presence and magnitude of a hiatus in the age model depends largely on which tie points are selected to correlate between Sites 1051 and 690 (Farley & Eltgroth, 2003) and is nonexistent-or considerably shorter-in some cases. Others have shown inconsistencies between the existence of this hiatus and the presence of correlative putative impact ejecta at the onset of the PETM at Site 1051 and others where a sedimentary hiatus is presumably not present (Schaller et al., 2016). We have chosen to align our Si microfossils ages using the extraterrestrial 3 He age model by Farley and Eltgroth (2003), without existence of a major hiatus. Further discussion on the implications of alternative age models for our interpretations can be found in section 4.
Sample preparation followed standard protocols. Briefly, sediments were cleaned with 30% H 2 O 2 and 1 M HCl to remove organic matter and carbonate material (Morley et al., 2004). Biogenic silica (bSi) was separated from detrital material by repeated heavy liquid separation using sodium polytungstate (density 2.1-2.3 g/ml). The light fraction containing bSi was then wet sieved at 53 μm to isolate intact radiolarian tests and sponge spicules that were handpicked under a light microscope. The radiolarians were well preserved and represent a bulk community , whereas only monoaxonic spicules were picked to avoid variability associated with varying morphologies (Cassarino et al., 2018;Hendry et al., 2015;Hendry et al., 2019). Microscopy shows no evidence for silica conversion to, or overgrowth by, clay minerals.
The isolated siliceous microfossils were dissolved in 0.4 M NaOH at 100°C in sealed polypropylene microcentrifuge tubes for a minimum of 3 days following , acidified to pH~2 using Suprapur HNO 3 , and purified chromatographically through a cation exchange resin after Georg et al. (2006).

Silicon Isotope Analysis
The silicon isotopic composition of the purified Si solutions was determined during two analytical sessions. During the first session, silicon isotopic compositions were analyzed using a Neptune MC-ICP-MS (Thermo Scientific) at the Pole Spectrometrie Ocean (Ifremer, Brest). Samples and standard solutions were diluted to 2 ppm Si with 1% HNO 3 and introduced to the Neptune through an APEX-HF (ESI) desolvating nebulizer.
These samples were doped with a matching concentration of Mg for mass-bias stabilization (Oelze et al., 2016) and to correct the isotopic ratios for internal mass bias (Cardinal et al., 2003). During the second session, we used a NuPlasma II MC-ICP-MS (Nu Instruments Ltd) at the Vegacenter in the Natural History Museum, Stockholm. Similarly, the Si solutions were diluted with 1% HNO 3 to 1-2.5 ppm and introduced in wet plasma mode through a glass nebulizer. Secondary standards (Diatomite and BHVO-2) on both instruments yield results in good agreement with accepted values (Reynolds et al., 2007) with long-term reproducibility <0.15 ‰ (2σ). A subset (n = 5) of full procedural replicates were analyzed during both sessions with a difference between each analysis lower than the long-term reproducibility. In all figures depicting δ 30 Si values, we use the long-term reproducibility as our uncertainty estimate, unless the internal measurement error of individual samples were higher.
Silicon isotope ratios are expressed relative to the standard NBS28 (NIST RM 8546) as δ 30 Si, corresponding to an average of several series of bracketed measurements (i.e., blocks of standard-sample-standard): where R sam and R std are the x Si/ 28 Si ratio, with x = 29 or 30, of the sample and the standard, respectively (after Mg correction, where appropriate). All values fall on the mass-dependent kinetic fractionation line δ 30 Si = 1.93 × δ 29 Si ( Figure S1 in the supporting information), demonstrating the successful removal of polyatomic interferences during measurement.
Between 6 and 15 sponge spicules from 37 samples were handpicked, embedded in epoxy and polished before being analyzed for δ 30 Si with a CAMECA-IMS 1280 secondary ion mass spectrometer (SIMS) at the Nordsim facility in the Natural History Museum, Stockholm. A critically focused primary Cs + ion beam with a current of 2.5 nA was used to sputter the secondary Si − ions with a raster size of 10 μm. A low-energy electron gun was used to minimize charge buildup. The Si − ions were introduced into the multicollector mass spectrometer, where 28 Si, 29 Si, and 30 Si were monitored on Faraday cups L1, H1, and H2, respectively. Samples were bracketed with analysis of grains of the quartz standard NBS28 under the assumption that matrix effects during material sputtering are negligible (Marin-Carbonne et al., 2012). Data reduction followed standard protocols detailed in Kleine et al. (2018).

Ocean Si Cycle Box Model
A two-box model, developed by De La Rocha and Bickle (2005) and modified by Frings et al. (2016), is used to track ocean [Si] and δ 30 Si as a function of time ( Figure 2). We use this model to test a series of scenarios, outlined in Table 1, related to aspects of the Earth system that are known or suggested to have altered before or during the PETM. We categorize these changes into four nonexclusive groups: (i) modification of the global ocean overturning patterns, (ii) modification of the delivery of Si to the ocean (both in magnitude and δ 30 Si), (iii) modification of ocean biogeochemistry, and (iv) volcanism. Given uncertainties in the ODP 1051 age model, the occurrence of transient states in the ocean Si cycle, and likely nonuniqueness of any solution, we chose to use a forward rather than inverse model approach. The aim is to create a framework that shows which processes Si isotope records are sensitive to, rather than to precisely reconstruct PETM silicon cycling.
The model simulates independent mass balances for the isotopes 28 Si and 30 Si in the ocean using an explicit, forward-in-time approach. Fluxes of Si are defined for total Si and then partitioned between the two isotopes according to the 30 Si/ 28 Si ratio in the source pool and any associated isotope fractionation. The two boxes schematically represent the euphotic zone and the deep ocean, respectively. Exchange of Si between the surface and bottom layers is achieved via a parameterization of physical upwelling and downwelling based on a prescribed oceanic mixing time (ca 1 kyr; De La Rocha & Bickle, 2005), and by the sinking and dissolution of bSi produced in the surface layer. The production of bSi in the surface layer follows Monod kinetics. This is functionally equivalent to a Michaelis-Menten formulation for enzyme kinetics. Production is thus a function of dSi concentration in the surface box: where Prod is the annual production of bSi in mol/y, V max is a maximum rate of bSi production in mol/y, K M is the half-saturation constant (i.e., the dSi concentration in mol/l at which the production of bSi occurs at half the maximum rate), and [dSi] surface is the dSi concentration in the surface waters in mol/l. This simple formulation, rooted in empirical observations and enzyme kinetics, is sufficient to capture plankton growth rates at a variety of space and timescales (e.g., Martin-Jezequel et al., 2000) and has been successfully used at global scales (Yool & Tyrrell, 2003. The total dissolution of bSi (Diss surface or Diss deep , in mol/y) in each box is a function of the residence time of bSi in the given box, the dissolution rate R diss , and bSi production Prod: where D is the depth of the box in m and V is the sinking velocity of particles in m/year. The dissolution rate of bSi (R diss ) is scaled linearly as a function of the degree of dSi undersaturation (Loucaides et al., 2012) and is expressed as ( 3) where k is a fitted constant that bundles together reactive surface area and opal intrinsic reactivity.
[dSi] eq is the apparent solubility of bSi in seawater in mol/l. The k value is tuned to fit observations that today about 3% of bSi production is ultimately buried in marine sediments and that about 50% of bSi dissolution taking place in the euphotic layer (De La Rocha & Bickle, 2005;Nelson et al., 1995). This formulation thus considers dissolution during sinking in the deep ocean and during early diagenesis together. Based on the above, a mass balance for 28 Si and 30 Si is calculated in each box at each time step (1/16 years). All material required to use the box model (MATLAB R2018b) is available from the corresponding author upon request. Table S1 summarizes the formulations, constants, and parameters used in the box model to test the scenarios outlined in Table 1.

Results
Using the age model from Farley and Eltgroth (2003) based on the assumption of a constant flux of extraterrestrial 3 He, the radiolarian test samples span from circa −70 to +346 kyr relative to the onset of the CIE. The sponge spicule samples range from −70 to +351 kyr relative to the onset of the CIE. The silicon isotopic composition of radiolarian samples ranges from +0.16‰ to +1.41‰ (Figures 3 and S2). Over the entire period, radiolarian δ 30 Si (δ 30 Si rad ) decreases slowly with time (δ 30 Si rad = 0.8635-0.0018 × time; r 2 = 0.48; n = 64). Focusing on the period around the CIE (−50 to +100 kyr relative to the onset of the CIE), there is an abrupt change in δ 30 Si rad coeval with the CIE (Figure 3 and S2). Pre-CIE, δ 30 Si rad averages +1.06‰ (1 sd = 0.15‰; n = 15). At the PETM onset, δ 30 Si rad decreases from about +1.2‰ to +0.5‰, after which they are constant with time (+0.65 ± 0.17‰, n = 24).
Sponge δ 30 Si (δ 30 Si sponge ) values range from −0.18‰ to +0.62‰. No significant trend at the α = 0.05 level could be identified over the entire period (r 2 = 0.21; n = 39) nor across the PETM specifically (r 2 = 0.08; n = 23). Insufficient sponge spicules and radiolarian tests between +18 and +36 kyr relative to the onset of the CIE precluded analysis in this period: the possibility thus exists that we may not have recovered the full magnitude of the δ 30 Si excursion associated with the CIE. Together with a potential sedimentary hiatus at Site 1051B (Röhl et al., 2000), the magnitude and duration of the shift in radiolarian δ 30 Si must be considered as minimums.
The 375 sponge spicule samples analyzed for δ 30 Si by SIMS produced values between −2.10‰ and + 2.62‰. Mean values ( Figure S3) for individual sample depths range between −0.62‰ and +1.69‰ but do not vary

Discussion
The competing explanations for the various enigmas of the PETM, including the carbon source, trigger, and removal mechanisms, can be recast as a series of hypotheses about the marine Si cycle-and thus predictions about the δ 30 Si of biogenic silica in the sedimentary record. A comparison of predictions from a model with measured δ 30 Si values is thus a direct test of the hypotheses and provides unique insights into the functioning of the Earth System during the PETM.
Several aspects have to be considered when investigating changes to the marine Si isotope budget, including the origin and magnitude of the change and its timing. At Blake Nose, the latter is difficult to address. Several constraints on the age model at Site 1051 exist, including cyclostratigraphy (Norris & Röhl, 1999), extraterrestrial 3 He fluxes (Farley & Eltgroth, 2003), and the presence of correlative impact ejecta (Schaller Paleoclimate context is given with carbon (black line) and oxygen (green line) stable isotope ratios from Site 1051 (Bains et al., 2000). Barite concentrations (yellow line) thought to indicate changes in paleoproductivity (Ma et al., 2014). Time is given relative to the onset of the PETM. et al., 2016). By comparison with Site 690 (Weddell Sea, Southern Ocean), it is plausible that part of the record at Site 1051 is missing at the onset of the PETM (Farley & Eltgroth, 2003;Röhl et al., 2000) and would result in a hiatus of up to 55 kyr. This is potentially due to a slumping caused by methane release in the vicinity of Site 1051 (Katz et al., 2001). Interestingly, the presence of spherules attributed to impact ejecta found at Site 1051 and two other locations north of Blake Nose challenge a sedimentary hiatus (Schaller et al., 2016). Schaller et al. (2016) argue that these spherules originate from an air-fall event and are contemporary. At all three locations, the peak abundance of spherules is coeval with the CIE such that if a hiatus were substantial, spherules at Site 1051 would be up to 55 kyr younger than elsewhere. The sedimentation history at Site 1051 is complex, and no consensus has yet emerged. Therefore, we focus our discussion explanations to the Si cycle without a large focus on replicating the timing. Interpretations based on the timing and transient effects at Site 1051 must be taken with caution and may evolve as more data from elsewhere become available. In the following, we summarize how the changes to the Earth system at the PETM may have affected the Si cycle.

Ocean Overturning
Changes to the overturning circulation have been invoked as a mechanism for the warming of deep ocean during the PETM (Bice & Marotzke, 2002;Lunt et al., 2010). The rationale is that a weakening of deepwater convection in the Southern Ocean displaced the loci of deepwater formation toward northern high latitudes. This had the potential effect of increasing the temperature of the deep ocean and destabilizing methane hydrates, releasing 13 C-depleted carbon into the atmosphere. In this way, ocean circulation could have acted as a positive feedback on an initial small warming (Bice & Marotzke, 2002;Nunes & Norris, 2006) or simply been a consequence of climate warming.
Empirical reconstructions are inconsistent, and in any case most focus on the Pacific, Southern, and South Atlantic oceans. Model studies generally suggest a weakening of deepwater formation in the Southern Ocean (Bice & Marotzke, 2002) that is supported by the carbonate ion interbasinal gradient (Zeebe & Zachos, 2007). However, there is little evidence for reorganization of the global ocean circulation over the PETM from Nd isotopes (Thomas et al., 2003). Nunes and Norris (2006) advance the possibility of multiple sources of deepwaters during the PETM.
If there were a reorganization of the ocean circulation, as suggested by the reversed carbonate ion gradient (Zeebe et al., 2009), the dSi concentration gradient should also be reversed relative to today (Penman, 2016). This would result in dSi-depleted North Pacific deepwater and dSi-enriched North Atlantic deep/intermediate water. Sponge spicule Si isotope ratios are well suited to substantiate this, given the relationship between their magnitude of Si isotope fractionation and the dSi concentration of the water they grow in (Hendry et al., 2019;Hendry & Robinson, 2012;Wille et al., 2010). This takes the form of an inverse relationship that asymptotes at a maximum fractionation of approximately −5‰ at high [Si] (Hendry et al., 2019;Hendry & Robinson, 2012). The inverse nature of the relationship also imparts greater sensitivity at lower dSi concentrations.
Previous analyses of sponge spicule δ 30 Si, also from ODP Site 1051, indicate that dSi concentrations in this region were low (<20 μmol/L) through the Cenozoic , a result robust to assumptions about seawater δ 30 Si. Today, there is a strong relationship between deepwater ventilation age and dSi concentration as diatom and radiolarian opal progressively dissolves (Sarmiento et al., 2007 increase, making the predicted decline even larger. Second, the required change of −1‰ to 2‰ is unrealistic: The modern Atlantic-Pacific deepwater gradient is about 0.3‰ (de Souza et al., 2014;Holzer & Brzezinski, 2015), and there is no known mechanism to drive deepwater to higher silicon isotope ratios.
We cannot conclusively rule out the possibility that the sponge spicules we analyzed have been laterally transported from shallower waters, which would equate to lower dSi concentrations and higher δ 30 Si values. If the fraction of transported spicules increased during the CIE, this could explain the stasis in the sponge spicule record. The unlikely coincidence required for this scenario to work notwithstanding, as discussed in Fontorbe et al. (2016), the sedimentology and foraminifera assemblages do not support substantial transport of biosiliceous material from upslope of Site 1051. The in situ SIMS δ 30 Si also argues against the occurrence of downslope transport. If the spicules do indeed derive from a range of depths, we would expect the distribution of silicon isotope ratios to reflect this. Instead, the range of values is no greater than that predicted by a typical SIMS δ 30 Si analytical reproducibility of approximately ±0.6-0.7‰ (2 σ) (Basile-Doelsch et al., 2005;Stefurak et al., 2015). Additionally, local reworking (cf. Katz et al., 2001) would smooth but not obliterate a signal.
The sponge data overall suggest that waters bathing the Western North Atlantic during the PETM remained dSi depleted, and so the source of bottom water likely remained nearby. Tripati and Elderfield (2005) argue that the emergence of a land bridge in the North Atlantic during the latest part of the Paleocene would have restricted water flow from the North Sea to the North Atlantic. It is therefore likely that the waters bathing Site 1051 originated from the North Atlantic from depleted surface waters.

Global Climate Change and Silicate Weathering
During the recovery from the PETM, δ 13 C and temperatures return to pre-PETM values. Even with constant input and output fluxes, δ 13 C would return to pre-event baseline values as the system is flushed (Dickens et al., 1997), but the return of climate (as proxied by δ 18 O) to pre-event conditions requires a net reduction in radiative forcing from atmospheric carbon to pre-event levels. The mechanisms proposed to remove carbon can be broadly classified into (i) increased carbon storage in the terrestrial biosphere (e.g., Beerling, 2000;Bowen & Zachos, 2010); (ii) increased marine organic carbon burial (Bains et al., 2000); and (iii) increased silicate weathering and subsequent carbonate-mineral burial (Ravizza et al., 2001). An increase in the size of the terrestrial biosphere is relatively neutral with regards to global silica cycling and is not discussed further. Increased marine productivity is addressed in section 4.3 below. Here, we focus on the potential for continental silicate weathering to explain our observations.

Increased Silicate Weathering: Expectations and Observations
Silicate weathering removes carbon from the atmosphere via a well-understood mechanism. CO 2 dissolves in rainwater to produce carbonic acid, a source of protons for silicate decomposition, epitomized by the hydrolysis of wollastonite: where the bicarbonate and dissolved cations subsequently recombine in the ocean to produce carbonate minerals, often biologically mediated. Thus, for every mole of silicate-hosted cation that is weathered, one mole of carbon is removed from the Earth's surface carbon cycle. Because the rate of these reactions is largely temperature, rainfall, and vegetation dependent, it functions as a negative feedback on atmospheric CO 2 (Maher & Chamberlain, 2014;Walker et al., 1981). dSi production by weathering should thus be climate dependent.
Weathering-derived dSi is ultimately transported to the ocean by rivers and groundwater, representing about 75% of the total inputs in the ocean Si budget Treguer et al., 1995;Tréguer & De La Rocha, 2013). Because the ocean Si budget can be conceptualized as a balance between higher δ 30 Si river water, and lower (near Bulk Silicate Earth (BSE) value) fluids resulting from the hydrothermal alteration of oceanic crust, mean ocean δ 30 Si reflects the balance between the two, modified by the fractionations associated with Si outputs (predominantly biogenic Si burial). A shift in the ratios of the two inputs would therefore be reflected in the Si isotope ratio of seawater-and of silica precipitated from it. A complicating factor is that the river endmember has variable δ 30 Si because silicate mineral weathering reactions are rarely as simple as the wollastonite example above (equation 5). They typically involve only the partial decomposition of the primary mineral to produce a secondary phase, as in the weathering of K-feldspar to kaolinite: where the proton is derived from carbonic acid or some other source. Experimental, theoretical, and empirical evidence all suggest that the partitioning of Si between secondary solid (in this case, kaolinite) and solute (i.e., dSi) is associated with Si isotope fractionation. For this reason, river water δ 30 Si is variable but has almost exclusively higher δ 30 Si than silicate rocks. The corresponding lower δ 30 Si counterpart required by mass balance is found in the secondary phases. River water δ 30 Si therefore varies with weathering congruency-the extent to which the conversion of primary minerals to solutes is complete. Changes in ocean δ 30 Si thus convolve both the magnitude and isotope ratio of the river dSi flux.
It is interesting to note that the direction of PETM radiolarian δ 30 Si change is opposite to that observed for biogenic silica records spanning the last deglacial period Sutton et al., 2018). The two periods are comparable in magnitude and rate, i.e. several degrees of global warming over several millennia. Yet a compilation of deglacial data from diatoms, radiolarians, and sponges suggests an average rate of change of +0.04‰ ka −1 (Sutton et al., 2018). Given the difficulties involved in defining precise age models for older sediments, a rate of change at the PETM is harder to calculate, but the direction is opposite in sign. For the deglacial, the two simplest explanations are an increase in river dSi flux at constant δ 30 Si, or an increase in river δ 30 Si at constant flux: most likely is some combination of the two. Both require an increase in silicate weathering rates as the planet warmed: the first to sustain an increased flux and the second because higher river δ 30 Si requires relatively more incorporation of weathered Si into secondary phases (i.e. clay minerals) and therefore greater initial solubilization of Si from parent rocks to maintain a constant flux. The face value interpretation of the PETM δ 30 Si is therefore a decrease in silicate weathering rates, which is improbable, since increases in the three key factors that govern silicate weathering rates-temperature, hydrology, and rock supply-are well documented at the PETM. For an activation energy of 60 kJ mol −1 , typical for silicate mineral dissolution reactions (Brantley, 2003;Oelkers et al., 2018;West et al., 2005), a 5°temperature increase should result in a reaction rate constant increase of about 50%, assuming an Arrhenius-style temperature dependency (i.e., k ∝ e − Ea =RT ð Þ , where k is the reaction rate constant, Ea the activation energy, R the universal gas constant, and T temperature in Kelvin). Recent reactive-transport models and empirical studies have emphasized the crucial role of water supply in weathering rates, particularly for so-called chemostatic systems where water-mineral interaction times are sufficiently long that the attainment of saturation can limit any further weathering (Maher & Chamberlain, 2014;Torres et al., 2015;Winnick & Maher, 2018). Finally, the most consistent observation from modern weathering studies is that weathering rates tend to increase with erosion rates and thus reactive mineral supply (Gaillardet et al., 1999;Riebe et al., 2001;West et al., 2005). There is good sedimentological evidence for increased erosion at the PETM (e.g., Pujalte et al., 2015;Schmitz et al., 2001). Overall then, there is good reason to expect a PETM increase in silicate weathering rates.
Yet strong evidence for an increase in silicate weathering rates over the PETM is surprisingly limited. Previous work with osmium isotopes from offshore of Svalbard (Wieczorek et al., 2013), the Indian and North Atlantic Oceans (ODP Sites 213 and 549; Ravizza et al., 2001), and the Arctic and peri-Tethys (Dickson et al., 2015) has produced inconsistent results. The ocean Os isotope budget is a balance between continental (radiogenic, high 187 Os/ 188 Os, and mantle (unradiogenic, low 187 Os/ 188 Os) end members. The prediction for the PETM would thus be for an increase in marine 187 Os/ 188 Os. While this is largely borne out by the data, converting 187 Os/ 188 Os ratios to CO 2 consumption fluxes is not easy. Interpretation of Os isotope ratios is complicated by the fact that osmium is sensitive to both flux and provenance changes: Organic-rich shales are rich in Re and can produce highly radiogenic dissolved Os. Recycling of Os in reducing marine sediments further complicates their interpretation as does the potential for an increased mantle contribution from North Atlantic volcanism (see section 4.4 below).
A second line of evidence put forward for enhanced silicate weathering during the PETM is the observation that kaolinite-a clay mineral devoid of cations and indicative of intense weathering-is found in greater abundance in PETM aged deposits than pre-or post-event. Localities with increased kaolinite include the north east American margin (Gibson et al., 2000;John et al., 2012), the North Sea (Kemp et al., 2016), the Basque Basin (Schmitz et al., 2001), the Tethys (Bolle & Adatte, 2001), and Antarctica (Robert & Kennett, 1994). However, there is a strong argument that clay mineralogy is insensitive to millennial-scale climate change (Frings, 2019;Thiry & Dupuis, 2000) because kaolinite-rich soil formation rates are slow. Sedimentological and palaeohydrological evidence for more intense rainfall events and physical erosion suggests the kaolinite peak could reflect transport of preexisting deposits (Schmitz & Pujalte, 2007). John et al. (2012) demonstrate with a mixing model that the clay δ 18 O from the New Jersey margin is compatible with this interpretation; that is, the kaolinite end member δ 18 O did not change. The emerging consensus is that the presence of kaolinite in marine sediments is more an indicator of reworking of previously formed laterites than contemporary clay formation, making the clay mineralogical evidence for increased silicate weathering rates during the PETM weak.
One important constraint on the magnitude of carbon released during the PETM is the degree to which the CaCO 3 preservation depth (i.e., the CCD) shallowed as the ocean carbonate saturation decreased. Enhanced supply of alkalinity and DIC from silicate weathering during the recovery should act to increase saturation and thus bury more carbonate. This may transiently deepen the CCD but could also be manifest as greater burial at shallow depths (cf. Greene et al., 2019). The presence of an CCD "overshoot" during the recovery period should be clear evidence for greater silicate weathering alkalinity supply and is a prediction of carbon cycle models (e.g., Zeebe et al., 2009). Penman et al. (2016) showed that the overshoot is present in IODP Site U1403 from the north Atlantic at a palaeo-depth of~4,400 m. Here, sediment devoid of CaCO 3 transitions to sediment with~30 wt% CaCO 3 , seemingly implicating silicate weathering as an enhanced alkalinity source. However, the lack of globally distributed sites with a CCD overshoot limits its use in quantifying any increase in silicate weathering rates. The CCD deepening at U1403 also appears to be permanent, questioning the extent to which it reflects a transient weathering increase versus a longer-term redistribution in the carbon cycle. Finally, Luo et al. (2016) showed how CaCO 3 overshooting can alternatively be explained by feedbacks internal to the ocean biogeochemical system. Enhanced carbonate burial rates above the CCD may also influence interpretations. Penman et al. (2016) showed that several sites globally exhibit increased CaCO 3 accumulation but extrapolating to the global ocean is difficult. A quantification of enhanced weathering from CaCO 3 burial rates awaits further work.
A final line of evidence that has been used to suggest increased silicate weathering rates over the PETM concerns the Si cycle itself. Penman et al. (2016) identified chert (microcrystalline quartz) layers on the Southeast Newfoundland Ridge (IODP Expedition 342) at the Paleocene-Eocene boundary that were interpreted to reflect the additional influx of 4.2-8.3 × 10 17 mol dSi. This mass was calculated to balance a total carbon input of 5,000-10,000 Gt (section 1), assuming a constant cation/silicon weathering ratio. The number of sites with chert layers in the Atlantic was subsequently expanded (Penman et al., 2019). However, their box-model calculations that incorporate a reversal in ocean overturning circulation-and thus a focusing of dSi in the North Atlantic (rather than the Pacific, as today; see above)-are largely sufficient to explain these observations without invoking an extra input of Si to the ocean. It is also conceivable that chert formation is promoted by the removal of any carbonate dilution effects, and the lack of pH buffering at higher values typical of carbonate systems (which would promote silica dissolution). While the presence of the chert layers is suggestive, it is also not conclusive evidence of enhanced continental silicate weathering. Penman et al. (2016) derived an increased river dSi flux of 2.8-5.6 × 10 12 mol Si year −1 , by distributing the equivalent carbon consumption flux over a 150 kyr period. The modern river dSi flux is around 6 × 10 12 mol year −1 (Beusen et al., 2009;Dürr et al., 2011), so this implies an increase in silicate weathering rates by 50-100% over the PETM, assuming a constant cation/silicon weathering ratio. The key finding of Penman et al. (2016) is the demonstration that the Si cycle is substantially impacted by changes to the carbon cycling across the PETM. Nevertheless, can the "back-of-the-envelope" estimates be refined? Caves et al. (2016) present a framework to predict the increase in weathering rates for a given pCO 2 (and thus climate) perturbation by defining the strength of the weathering feedback as a constant of proportionality k relating the weathering flux to climate:

How Much of a Weathering Increase Is Expected?
where a change in k reflects a change in the "weatherability" of the Earth surface, that is, the efficiency with which the silicate weathering flux, F silw , is produced at R CO2 , a pCO 2 normalized to a modern value, taken here to be 300 ppmv. By combining a reconstruction of silicate weathering fluxes with atmospheric pCO 2 estimates, Caves et al. (2016) derive a late Paleocene k of approximately 50% the modern value (cf. Froelich & Misra, 2014), explaining the early Cenozoic high pCO 2 . Taking a modern F silw of 10 × 10 12 mol C year −1 (which Caves et al., 2016, show cannot have varied by more than a few percent, assuming constant degassing), and a late-Paleocene pCO 2 of 800 ppmv (i.e., R CO2 ≈ 2.7) (Anagnostou et al., 2016), then late Paleocene k ≈ 3.3 × 10 12 mol CO 2 year −1 . Ignoring any short-term change in global weatherability associated with altered rainfall distributions and lithological effects (discussed in Ibarra et al., 2016), we can predict the expected change in F silw over the PETM. For peak PETM atmospheric pCO 2 of~2,000 ppmv (i.e., R CO2 ≈ 6.7) (Gutjahr et al., 2017;Schubert & Jahren, 2013), then peak F silw would equal 15 × 10 12 mol C year −1 , an~50% increase. This is less than the doubling that Penman et al. (2016) invokes, and a time integrated average for the duration of the CIE would be lower still, since proxy-driven model reconstructions of atmospheric pCO 2 over the PETM predict CO 2 does not remain elevated for the entire duration.
Assuming that the other inputs to the ocean Si budget (dust and sediment dissolution, hydrothermal fluid recirculation, and submarine groundwater discharge) remain constant, Figure 4 shows the sensitivity of steady-state ocean δ 30 Si to changes in river dSi flux and δ 30 Si. The contour depicting the required~0.6‰ drop is highlighted in red. The takeaway message is that it is hard to reproduce the observations just by tweaking the river parameters. An almost complete cessation of river dSi delivery would explain the observations but is clearly impossible. A simple weathering rate increase, at constant congruency, would act in the wrong direction and increase ocean δ 30 Si. These conclusions are corroborated by our ocean Si box model (Figure 5a; Scenarios A and B in Table 1), which allows nonsteady-state responses to be recognized. In line with previous observations, increased river dSi fluxes result in 5-10% and 3-4% increases in dSi concentrations in the surface and deep ocean, respectively Penman, 2016). These scenarios also predict an increase in δ 30 Si of dSi in the surface waters by about 0.10-0.17‰. This response is incompatible with the radiolarian δ 30 Si measurements, meaning a simple increase in the continental weathering rate cannot explain the observations. Apart from changes in delivery fluxes, another aspect of altered continental chemical weathering considers a shift in weathering congruency, that is, the completeness of chemical weathering. Figure 4 suggests the most effective pathway would be to maintain a constant CO 2 consumption flux but with more congruent reaction stoichiometries (Line 1 in Figure 4b, where the slope is related to the mean fractionation associated with clay formation). Any increase in silicate mineral weathering rates (and therefore initial Si mobilization) would act to lessen the gradient of Line 1 and reduce the efficiency by which increased weathering congruency would lower the weighted mean δ 30 Si of all ocean Si inputs.
In the modern world, there is a latitudinal and altitudinal zonation of clay mineralogies, with the most Si-deplete clays-i.e., those formed under the most congruent weathering conditions-occurring in warm and wet environments (Chamley, 1989;Fagel, 2007). The overall expectation of PETM weathering would thus be a shift toward more congruent weathering, producing more Si-deplete secondary phases. Assuming a constant fractionation factor associated with clay formation, mass balance dictates that the dissolved pool must become closer to the parent material. This is observed today in the lowland, "blackwater" tributaries of the Congo and the Amazon Hughes et al., 2013), where the net weathering reactions are extremely congruent with respect to silicon, and so river dSi has δ 30 Si values close to 0‰ (with bedrock typically lying around −0.2‰). The important implication here is that river dSi fluxes can increase without necessitating an increase in CO 2 removal by silicate weathering.
A further nuance comes from considering the state of the pre-PETM landscape. It is generally considered that large parts of the Paleocene world were blanketed in thick, lateritic deposits (Bolle & Adatte, 2001;John et al., 2012). These contain large pools of kaolinite, thought to be the source of the kaolinite peak in 10.1029/2020PA003873

Paleoceanography and Paleoclimatology
FONTORBE ET AL. marine sediments described above. Today, typical δ 30 Si values for kaolin-group minerals are around −2‰ Ziegler, Chadwick, Brzezinski, & Kelly, 2005;Ziegler, Chadwick, White, & Brzezinski, 2005). If some fractions of these low-δ 30 Si secondary solid phases were to dissolve, it would drive the riverine Si isotopic composition toward lower values. Dissolution of bSi produced by terrestrial vegetation (phytoliths) and stored in soils may also be a source of low-δ 30 Si silicon. Importantly, because there are no cations to be lost during the decomposition of clay minerals or phytoliths, there is no associated CO 2 consumption, and so the silicon and carbon cycles become partly decoupled.
Overall, both greater weathering congruency and/or dissolution of preexisting clay minerals or soil bSi are plausible mechanisms for lowering the δ 30 Si of river dSi. It is difficult to specify likely values for this lowering, so we prescribe river dSi δ 30 Si decreases of 0.5‰ and 1.0‰ in our model for illustrative purposes. The δ 30 Si of ocean surface water decreases with lower river δ 30 Si (Figure 5b): A drop of 0.5‰ produces about a third of the total magnitude of the radiolarian δ 30 Si excursion, while the drop of 1.0‰ reproduces the full event. Note that we do not attach great significance to replicating the timing of our records, given the uncertainties in the age model and the potential hiatus at Site 1051 (Katz et al., 2001). Si. White star represents the pre-PETM steady state (control values in Table 1). Contours are calculated using a modern-day ocean budget as outlined in Frings et al. (2016). The approximate magnitude of ocean δ 30 Si decline at the PETM as recorded in Blake Nose radiolarians (0.6‰) is highlighted in red. The implication is that to decrease by ocean δ 30 Si by the required amount river δ 30 Si is a more effective parameter than river Si flux. Shaded panels below the plot indicate the approximate weathering rate increases suggested by the Caves et al. (2016) framework and by Penman (2016). (b) Conceptual expectations of vectors to be followed given scenarios described in the main text. Line 1 describes a change in weathering congruency at a fixed solubilization flux, while Line 2 describes a change in solubilization (weathering) flux at fixed congruency.
Overall, if we wish to interpret our data purely in terms of continental weathering, the simplest interpretation requires that the increase in CO 2 consumption by silicate weathering in response to the PETM was small and that the major changes were more to the style of weathering. This is in line with (i) the prediction of a weak early Cenozoic thermostat , (ii) model considerations that suggest the Paleocene silicate-weathering feedback was too weak to account for the rapidity of recovery from the CIE (Bowen & Zachos, 2010;Gutjahr et al., 2017), and (iii) recent evidence that suggests feedbacks in the organic carbon cycle are strongly implicated in the recovery (Bridgestock et al., 2019;Frieling et al., 2019). The greater the CO 2 consumption rate increase by silicate weathering, the more problematic the interpretation of the radiolarian data becomes. Thus, we favor a scenario where river Si fluxes increase in magnitude and decrease in δ 30 Si due to a combination of increased weathering congruency and leaching of Si from preexisting laterite deposits (Figures 5a and 5b). Note that this does not negate the existence of the silicate-weathering thermostat. Rather, we suggest future model-data inversions could use silicon isotope records as a constraint on the efficiency of the thermostat.

Ocean Biogeochemistry
Sequestration of atmospheric carbon during the recovery phase of the PETM may have been achieved-at least partially-via the increased productivity of marine photosynthetic organisms (Bains et al., 2000). and 10-fold, 5.6 Tmol/year, and −1‰ (purple). Reference to the given scenario (Table 1) is given as a letter in the legend.

Paleoceanography and Paleoclimatology
However, different groups of organisms exhibit differing responses to the PETM. Most notable is the extinction of many benthic foraminifera, with estimates of diversity decreases as large as 50%. Probable culprits include decreasing oceanic pH and oxygen saturation (Dickson et al., 2012), increasing deep ocean temperature and modifications to the food supply (Alegret et al., 2009;Alegret et al., 2010;Aref & Youssef, 2004;Thomas, 1989). Conversely, there is no planktic foraminifera extinction event, though poleward migration and rapid evolutionary diversification are known (Arenillas et al., 1999;Kelly et al., 1996;Pardo et al., 1999). Overall, increased productivity has been argued for on the basis of an increase in nondetrital Ba and barium isotope ratios, showing an increase in calcifying nannoplankton production (Bridgestock et al., 2019) or increased export efficiency in general (Ma et al., 2014), albeit with inconsistent timings of production increase. This is corroborated by Mo isotopes, which suggest a slight decrease in ocean oxygenation over the PETM as the ocean warmed and more organic carbon was oxidized (Dickson et al., 2012).
Although data are limited, there is enough to resolve regional discrepancies in the responses of siliceous organisms to the PETM. No radiolarian mass extinction is observed during the PETM, but sites appear to differ in species turnover. For example, the Mead Stream section (New Zealand) exhibits marked radiolarian faunal changes over the PETM with migration of warm-water species toward higher latitudes (Hollis, 2007) with assemblages indicating the expansion of "weakly eutrophic" waters into the Southern Pacific Ocean. In a low-resolution Cuban section, preboundary radiolarian taxa dominate until the middle of foraminiferal biozone P6b, well into the Ypresian (Sanfilippo & Hull, 1999). At Blake Nose, Sanfilippo and Blome (2001) did not report any major change in the composition of the radiolarian fauna over the PETM apart from a transient decrease in total abundance and preservation. In general, the migration and assemblage changes better known from other groups also affect radiolarian taxa to various degrees. Is it plausible that the radiolarian δ 30 Si excursion reflects ecological change? Simplifying, radiolarian δ 30 Si signatures are a function of (i) ambient water δ 30 Si and (ii) the fractionation associated with radiolarian test production. This means that any biological change (whether related or not to radiolarian biology) that impacts the silicon isotopic composition of surface waters has the potential to be recorded in radiolarian tests.
Diatoms were likely key players in the ocean Si cycle in the Early Cenozoic Renaudie et al., 2018), so we can use insights from the modern ocean to guide our interpretations. Today, sea surface δ 30 Si is primarily driven by diatom Si uptake and precipitation of their frustules. This is associated with a fractionation, canonically placed at around −1‰ to −2‰ (De La Rocha et al., 1997;Sutton et al., 2013). Increasing bio-utilization of dSi thus drives the residual dSi toward higher δ 30 Si values. The general understanding is that as heterotrophs, radiolarian bSi production is only a small fraction of total siliceous production, making them a passive recorder of ambient water δ 30 Si. Within this framework, a large and rapid decrease in diatom productivity could shift surface water δ 30 Si closer to the weighted mean of the inputs and serve to explain the measured radiolarian δ 30 Si excursion. A modern weighted average of the Si inputs to the ocean is~0.75‰ , while average ocean surface water δ 30 Si is~2.0‰-roughly one diatom fractionation offset. As mentioned above, an increased loading of dSi from the continents during the recovery of the PETM would enhance diatom bSi burial. Our box model includes formulations of the feedbacks between bSi production and dissolution, and dSi concentration. In order to recreate a surface water δ 30 Si drop of~0.7‰, the model requires a collapse in gross siliceous production.
The problem with this hypothesis is that there is no known extinction event for diatoms or radiolarians over the PETM, though intensified species turnover is suggested for sites on the epicontinental Russian Sea (Speijer et al., 2012). If anything, diatom production likely contributes to a general increase in organic carbon burial as required by carbon cycle models (e.g., Gutjahr et al., 2017) and perhaps evidenced by biogenic barium records (e.g., Ma et al., 2014), though Frieling et al. (2019) emphasize this could also be driven by increased Ba inputs to the ocean from methane hydrate dissociation. By increasing the relative utilization of dSi, this would act to increase surface water δ 30 Si (see, e.g., Fripiat et al., 2011, for a modern example) and thus increase radiolarian δ 30 Si, all else being equal. Overall, we find that ecological scenarios producing a decrease in surface water δ 30 Si of a similar magnitude to the radiolarian record very unlikely. The most likely is in the opposite direction: Enhanced diatom production would increase surface water δ 30 Si.
An alternative explanation would be that some facet of radiolarian Si isotope fractionation changed over the PETM. Only a handful of studies have tackled Si isotopic fractionation during precipitation of radiolarian tests, and little is known about the mechanisms controlling this fractionation (Abelmann et al., 2015;. Abelmann et al. (2015) showed that two different radiolarian size fractions (125-250 μm and >250 μm) in the Southern Ocean had an offset in their isotopic composition of~1‰. The cause of this offset has yet to be resolved: Potential explanations include a systematic difference in the depth habitat of the different size fractions, or to different fractionations associated with different radiolarian species. A compilation of plankton tow and sediment trap data suggests that most living radiolarians are found today in the upper 400 m, with both species richness and cell densities dropping sharply below this (Boltovskoy, 2017). This suggests there is little scope for a large change in depth habitat, and thus ambient δ 30 Si being recorded. Additionally, Sanfilippo and Blome (2001) do not report a major shift in radiolarian assemblages during the PETM at Blake Nose (IODP Hole 1051A). The PETM interval lies within the radiolarian zone Bekoma bildartensis RP7, which is not marked by a major reorganization of the composition of the radiolarian fauna (i.e., abundant species remain abundant during the majority of the zone, and rare species remain rare).
Other possibilities that would require detailed paleoecological work to begin to address include the potential for (i) a change in seasonality of radiolarian growth, (ii) a change in radiolarian depth habitat, or (iii) a change in silicification and thus fractionation, perhaps associated with micronutrient availability. While we conclusively cannot rule out the influence of ecological change on the radiolarian δ 30 Si record, the absence of evidence for a major reorganization of the radiolarian community at our study site argues against an ecological cause for the rapid decrease in radiolarian δ 30 Si. Future work could assess what drives interspecies and intraspecies Si isotope fractionation variability as a first step toward assessing how much radiolarian paleoecology may be influencing our observations.

Solid Earth Interactions
Increased volcanic activity has been invoked to explain various aspects of the PETM. These include (i) volcanic/magmatic CO 2 as a cause for pre-CIE warming and thus an indirect trigger for the PETM (e.g., Frieling et al., 2019), (ii) volcanic/magmatic activity as a direct trigger for 13 C-depleted carbon release (e.g., Svensen et al., 2004), and (iii) volcanic CO 2 as the dominant cause of the CIE (e.g., Gutjahr et al., 2017). The specific implications these suggestions have for the carbon and silicon cycles depend on whether subaerial or submarine volcanism is invoked.

Pre-CIE Warming
There is good evidence for warming of about 2°C that precedes the CIE by millennia (e.g., Frieling et al., 2019;Sluijs et al., 2007;Thomas et al., 2002). This evidence takes the form of sea surface temperature (SST; reconstructed from δ 18 O or the TEX 86 paleothermometer) increases that stratigraphically precede the CIE. This was interpreted as evidence for the crossing of orbitally driven insolation thresholds that led to a positive carbon cycle feedback (e.g., DeConto et al., 2012). Recent orbital calculations provide inconsistent support for the onset of the PETM occurring in phase with the dominant 405 ka eccentricity cycle Zeebe & Lourens, 2019), and a high-resolution, orbitally tuned benthic foraminifera δ 13 C record seems to argue against forcing in the 405 ka band (Littler et al., 2014;), but thẽ 5 ka lead of temperature (i.e., δ 18 O) of carbon cycling (δ 13 C) is typical of the late Paleocene. Regardless, the apparent global nature of the temperature-CIE lead/lag relationship remains strongly suggestive of a CO 2 forcing (Frieling et al., 2019). If this reflects input to the Earth surface, rather than repartitioning among the ocean-atmosphere carbon pools, the consistency in the pre-PETM carbon isotope record requires that the carbon source is isotopically similar to the mean Earth surface carbon pool (approximately −2‰; Dickens, 2001); attention has therefore focused on volcanically derived mantle carbon, with a δ 13 C of approximately −5‰. In particular, the North Atlantic Igneous Province (NAIP), associated with the mid-Paleocene to early-Eocene opening of the Greenland-Norwegian Sea, has been implicated.
Subaerial volcanism (Wieczorek et al., 2013) would be relatively neutral with regards to ocean Si isotope budgets, minor effects from ash weathering notwithstanding. Submarine volcanism including magmatic intrusion into sediments induces large thermal gradients that can cause extensive hydrothermal activity . Svensen et al. (2004) identified thousands of hydrothermal vent complexes associated with the Vøring and Møre basins offshore Norway. These complexes developed from contact aureoles around sill intrusions by explosive release of fluids and sediments within decades of sill emplacement , and their temperature gradients cause seawater to circulate through them. The interaction of high-temperature fluids with basalts or sediments generates high [Si], low δ 30 Si fluids (De La Rocha et al., 2000;Geilert et al., 2015), so the stasis in North Atlantic radiolarian δ 30 Si pre-CIE places limits on the allowable increase in hydrothermally sourced Si. Taking ±0.2‰ in δ 30 Si as a detection limit, our box model constrains any pre-PETM hydrothermal Si flux increase to below a factor of 3 ( Figure 5, Scenario K). This factor still holds for a restricted Atlantic basin, if the ratio of continental to hydrothermal Si is the same as the global value. Converting this to a CO 2 degassing rate is not straightforward, but taking modern values for the input fluxes from mid-ocean ridges in the long-term carbon cycle (3 × 10 12 mol C year −1 ; Dasgupta & Hirschmann, 2010), then a threefold increase is approximately 0.1 Pg C year −1 , much less than PETM onset rates (Gutjahr et al., 2017), human fossil fuel burning (Le Quéré et al., 2018), or LIP emplacement (Schaller et al., 2011) This argues against mantle-sourced carbon degassing through seafloor vent systems as the cause of pre-PETM warming. Reconciling the observed pre-PETM warming with the lack of a CIE remains a challenge. 4.4.2. Volcanically Triggered Carbon Release Progressively more precise age determinations of individual sills or ash layers demonstrate that pulses of volcanic activity in the North Atlantic region were coeval, within error, of the PETM (e.g., Storey et al., 2007;Svensen et al., 2004;Svensen et al., 2010). One vent complex has been dated to within the body of the PETM (Frieling et al. 2016). Most work has focused on the generation of greenhouse gases in contact aureole as sills were emplaced into organic rich sediments (e.g., Svensen et al., 2004), leading to organic carbon remineralization and/or thermogenic methanogenesis (Storey et al., 2007), causing the CIE and most of the PETM warming. Jones et al. (2019) recently showed that this mechanism can plausibly produce carbon at a sufficient rate. Converting this to a dSi flux would require knowledge of the compositional change of intruded sediment, but our box model ( Figure 5, Scenario O) suggests a~100 times increase in hydrothermally sourced Si is required to explain the radiolarian δ 30 Si excursion. This factor would increase if late Paleocene river δ 30 Si were lower than today, as might be expected. Jones et al. (2019) calculate peak excess magmatic carbon release rates from the NAIP of~2 × 10 12 mol C year −1 . This is of the same order of magnitude as modern magmatic MOR degassing estimates (Dasgupta & Hirschmann, 2010; see above). We thus suggest that a >100 times increase in hydrothermally sourced Si is hard to support, under the assumption that Si/C release ratio of intruding sills and MORs are not orders of magnitude different. The radiolarian δ 30 Si excursion is therefore compatible with volcanically triggered CO 2 of CH 4 release, but this volcanic activity itself is unlikely to explain the excursion. 4.4.3. A Dominantly Volcanic Carbon Source for the CIE Gutjahr et al. (2017) suggest that volcanic carbon from the NAIP was the dominant driver of the CIE. By inverting the cGENIE model with a boron isotope pH record, they derived peak C emission rates of 0.58 Pg C year −1 with a relatively high δ 13 C (approximately −11‰). This is interpreted as reflecting 75-90% C from mantle sources, with more depleted Earth surface reservoirs constituting the remainder. This is a factor of 16 times greater than modern MOR CO 2 flux. Increasing the hydrothermal flux (i.e., if the C emission were all submarine) by this factor in our box model cannot reproduce the observed radiolarian δ 30 Si excursion (Figure 5d, Model Scenarios L and M). As with a volcanic trigger, a dominant volcanic source of carbon is compatible with the silicon isotope record but likely does not explain it.

Conclusions
There are many facets of the PETM that are generally accepted: There was some pre-event warming that likely pushed some aspect of the Earth system over a threshold and triggered the release of a large amount of 13 C depleted carbon. This carbon led to global warming before it was eventually removed by negative feedbacks within the carbon cycle. The PETM is an interesting analog for anthropogenic climate change that has generated many hypotheses to explain various aspects. Here, we exploited the observation that many of the proposed hypotheses make testable predictions about the marine silicon cycle. Our data from Blake Nose in the western North Atlantic suggest that it is unlikely the pre-event warming resulted from enhanced submarine degassing of magmatic CO 2 (otherwise, enhanced hydrothermal fluid recirculation would lower ocean δ 30 Si). However, it is compatible with a magmatic trigger for, or source of, CO 2 for the main event. The data do not support a change in regional oceanography. Finally, it suggests that the strength of the weathering feedback was weak (otherwise, we would expect the ocean Si isotope budget to shift to a more continentally derived state) and that organic carbon burial had an important role in the recovery, though this organic carbon was likely not produced by silicifying organisms. The challenge of the PETM and the other Cenozoic hyperthermals is to find an internally consistent explanation that can reconcile the apparently discrepant ecological, mineralogical, and geochemical evidence. Given this, and the age-model issues at ODP Site 1051, a globally expanded data set of marine silicon isotope ratios over the PETM or other Cenozoic hyperthermals would help address many of the outstanding issues.