Constraining climate sensitivity and continental versus seafloor weathering using an inverse geological carbon cycle model

The relative influences of tectonics, continental weathering and seafloor weathering in controlling the geological carbon cycle are unknown. Here we develop a new carbon cycle model that explicitly captures the kinetics of seafloor weathering to investigate carbon fluxes and the evolution of atmospheric CO2 and ocean pH since 100 Myr ago. We compare model outputs to proxy data, and rigorously constrain model parameters using Bayesian inverse methods. Assuming our forward model is an accurate representation of the carbon cycle, to fit proxies the temperature dependence of continental weathering must be weaker than commonly assumed. We find that 15–31 °C (1σ) surface warming is required to double the continental weathering flux, versus 3–10 °C in previous work. In addition, continental weatherability has increased 1.7–3.3 times since 100 Myr ago, demanding explanation by uplift and sea-level changes. The average Earth system climate sensitivity is  K (1σ) per CO2 doubling, which is notably higher than fast-feedback estimates. These conclusions are robust to assumptions about outgassing, modern fluxes and seafloor weathering kinetics.

We see that for the range of temperature variations we are considering, any weathering function with a runoff term can be fitted with a single exponential, and hence it is appropriate to model the overall temperature dependence of continental weathering using an exponential. Supplementary

4). Grey and red
shaded regions represent the model output 90% confidence obtained from 10,000 forward model runs using the parameter ranges described in Table 1. The grey and red solid lines are the median model outputs. Black and red dots represent binned geochemical proxy data, and error bars denote the range of binned proxy estimates (see main text for references and explanation). Here, the model envelopes marginally encompass the proxy data. The upper end of the temperature and seafloor envelopes fit proxies, pCO 2 is an excellent fit, and the saturation state and pH proxies are slightly outside the envelope.
Supplementary show how each pair of variables co-varies. The joint distributions show degeneracies that could be resolved with better data. For example, there is a positive correlation between Cretaceous weatherability and Cretaceous outgassing; if outgassing at 100 Ma was high, then the weatherability change was necessarily modest, and vice versa. This degeneracy highlights uncertainty surrounding Cretaceous climate: whether high pCO 2 levels were caused by enhanced outgassing or reduced CO 2 sinks. If either variable could be constrained by data, then the retrieval would be more informative.
Supplementary soil pCO 2 is used instead, adopting the model of Volk 1 to link soil and atmospheric pCO 2 , and assuming a maximum biosphere productivity of 4 times the modern productivity. We see that for our purposes it is valid to use atmospheric CO 2 in equation (2) Here, t is in years since the present, and [Ca 2+ ] is in mol kg -1 . This expression was used for the calcium abundance in both the ocean and the pore-space.
Supplementary Figure 9: A selection of GCM outputs (colored symbols) and simple climate models (lines) from the literature. The grey shaded region is the range of climate parameterizations considered in our model, equation (9), where the climate sensitivity, . The red circles are the 10 Ma binned data used in this study. The uncertainty envelope for each binned data point is taken to be the range of proxy values within each 10 Ma bin.
Supplementary Figure 11: Black squares show proxy estimates of atmospheric pCO 2 from the literature [10][11][12][13][14][15][16] . The red circles are the 10 Ma binned data used in this study. The uncertainty envelope for each binned data point is taken to be the range of proxy values within each 10 Ma bin.
Supplementary Figure 12: Black and red circles show binned mean temperature values for the deep ocean and the surface, respectively. The proxy data used to construct these ranges is described in Supplementary Methods. Black crosses show Cenozoic deep ocean temperatures from Hansen et al. 17 .
Supplementary Figure 13: Colored lines show estimates of Calcite Compensation Depth (CCD) for various ocean basins from the literature 3,[18][19][20][21][22] . The red circles are the 10 Ma binned data used in this study. The uncertainty envelope for each binned data point is taken to be the range of reconstructions within that 10 Ma bin.
Supplementary Figure 14: Comparison of dynamical model (solid lines) and steady state calculations (circles). Supplementary Note 3 describes steady state calculations in full detail. There is excellent agreement between the two methods.
Supplementary Figure 15: An additional alkalinity source has been introduced in the pore space to account for K-feldspar formation (see Supplementary Note 4). Grey and red shaded regions represent the model output 90% confidence obtained from 10,000 forward model runs using the parameter ranges described in Table 1. The grey and red solid lines are the median model outputs. Black and red dots represent binned geochemical proxy data, and error bars denote the range of binned proxy estimates (see main text for references and explanation). The fit with pore-space precipitation is improved at the expense of a temperature and pH mismatch.
Supplementary Figure 16: Results of Bayesian inversion with limited data set. Each variable is fitted to a single mid Cretaceous data point (Cenozoic data point used for pH since no Cretaceous data exist). 1σ uncertainties in each data point are plotted, and shaded regions represent 95% credible intervals.
Supplementary Figure 17: Selected posterior probability distributions for the Bayesian inversion in Supplementary Fig. 16. Dotted lines represent the median value with 1σ error bars. The offdiagonal elements are joint probability distributions, which show how each pair of variables covaries. Even when fitting a single mid Cretaceous data point for each variableand thereby limiting the potential biases from underfitting and linearizing parameters -the general conclusions reported in the main text remain unchanged. Specifically, e T is likely to be large, relative Cretaceous weatherability, 1+W, is approximately half modern weatherability,

Supplementary Note 1: Climate and deep-ocean temperature parameterizations
The Global Circulation Model (GCM) outputs used in Fig. 8 were taken from Li et al. 23 , Stouffer and Manabe 24 , and Danabasoglu and Gent 25 . These studies were chosen because they used fully coupled atmosphere-ocean GCMs with complete ocean circulation. Each GCM was run for thousands of years, which is sufficient time for the deep ocean to reach equilibrium.
The peak-PETM mean surface temperature, 27.3°C, was obtained by averaging all the reported temperatures in Table 1 of Jones et al. 26 , and weighting each measurement by the cosine of its latitude. The pre-PETM mean surface temperature, 23.03°C, was obtained by subtracting the best-fit warming (4.3°C) from the peak-PETM temperature, although directly averaging pre-PETM measurements produced a similar result. Deep ocean pre-PETM and peak-PETM temperatures were estimated from Table 2 of Jones et al. 26 . Observed deep ocean temperatures were weighted by the areas of their respective ocean basins, and the pre-PETM and peak-PETM averages were found to be 12.1°C and 16.5°C, respectively. Error bars for PETM proxy temperatures were estimated as follows: McInerney and Wing 27 and Higgins and Schrag 28 estimate 5-8°C surface warming during the PETM, and so we adopted an error of 2°C in both the peak-PETM and pre-PETM surface temperatures.
Last Glacial Maximum (LGM) temperatures were estimated as follows. Clark et al. 29

Supplementary Note 2: Climate model
The solar luminosity term in our climate model (equation (9)) was derived as follows. Doubling atmospheric CO 2 is equivalent to a 3.7 W m -2 radiative forcing 39 . Taking into account the geometry of insolation and assuming an albedo of 0.3, the solar luminosity decrease required to offset a 3.7 W m -2 radiative forcing is: Following Catling and Kasting 40 the evolution of relative solar luminosity, L , is approximated by: We seek the time, However, mean surfaces are slightly more sensitive to pCO 2 changes than to luminosity changes because the CO 2 forcing is more effective at high latitudes 41 . Hansen et al. 41 reported CO 2 forcings to be 1.26 times more effective than luminosity forcings and so we adopt this correction in our climate equation (neglecting paleogeography): The linear paleogeography term is described in the main text.

Supplementary Note 3: Validation of dynamical model with steady-state calculations
To validate our model we compared outputs to equivalent steady state calculations. To make the steady state calculations analytically tractable, it is necessary to simplify the dynamical equations somewhat. Many different cases were validated, and here we present one illustrative example. In this example the only sources/sinks are carbon outgassing, pore-space precipitation, and pore space basalt dissolution, i.e. no carbonate or silicate weathering. We also assume 0  s . Equation (6) can therefore be simplified to: The model was forced by increasing outgassing from modern to 4.5x modern levels over several billion years. Rather than specify calcium abundances, in this example we let 2+ Ca   evolve dynamically with alkalinity: Additionally, we ignore the pH dependence of basalt dissolution and assume dissolution is purely a function of temperature: . Selected model outputs from this dynamical calculation are plotted in Supplementary Fig. 14.
The same problem was also solved by computing successive steady states. Starting from the initial (modern) steady state for the ocean and pore space, outgassing was incrementally increased up to 4.5x modern levels, and the new steady state was found at each outgassing level. This calculation is described in full below. In the equations that follow, P-subscripts denote pore-space variables, and O-subscripts denote ocean variables.  (9) and (12).
Similarly, since pore P is known, P  can be calculated from equation (19). In other words the steady state conditions allow us to determine 2 pCO and P  for any given level of outgassing, and from these two variables the complete carbon chemistry of the ocean and pore space can be solved simultaneously. By substituting equations (23) and (24) into equations (21) we obtain: Next, the first expression in equation (S7) can be rearranged and substituted into the third expression as follows: Denote the solution to this quadratic   (S11), and the remaining carbon chemistry is trivially solved using equations (21) to (26) in the main text. The results from this steady state calculation are plotted in Supplementary Fig. 14 and compared to equivalent dynamical model outputs. They are in agreement, which implies that the numerical integration is working correctly, and that the carbon cycle is in quasi-steady state.

Supplementary Note 4: Modified models -K-feldspar uptake
A caveat on our results is that the assumed functional relationship for seafloor weathering (equation (13) For typical changes in crustal production, mod out out 1.5 FF  , the alkalinity source due to Kfeldspar formation will be ~2 Tmol eq yr -1 at 100 Ma and 0 Tmol eq yr -1 in the modern porespace.
Including this alkalinity source does not affect our conclusions. For the base case (no weatherability change, W=0, and conventional temperature dependence of silicate weathering, 5 15 K e T  ) the seafloor precipitation flux at 100 Ma is a better fit with data, however the fit with temperature and CO 2 is worsened. The combination of a strong silicate weathering feedback and a large seafloor sink due to K-feldspar formation draws down too much pCO 2 and makes the Cretaceous climate unreasonably cold ( Supplementary Fig. 15). If either the temperature sensitivity of silicate weathering is weakened or a large weatherability increase is imposed, then the fit with temperature and CO 2 is marginal. The best fit is achieved by assuming both a weak climate sensitivity ( ) and a large weatherability change (W=0.4-0.6), not shown. In short, a sizeable change in weatherability and/or a low temperature dependence of silicate weathering is still required to fit proxy data.

Supplementary Note 5: Modified models -Fitting the mid Cretaceous mean state
One potential criticism of our methodology is that by imposing linear trends on some variables in our model we do not fully capture shorter timescale fluctuations. By fitting binned time series data in this way we are arguably underfitting and potentially underestimating the true uncertainty in unknown variables. To test this we repeated the MCMC inversion but simplified the data to include only one mid Cretaceous data point for each variable. By only fitting the mid Cretaceous endpoint the misfit from imposing linear trends is minimized. In this calculation we are no longer fitting a time series but rather fitting mean Cretaceous conditions with broad uncertainties.
The model outputs for this calculation are shown in supplementary figures 16 and 17. Although the uncertainties in fitted parameters are somewhat larger than the nominal case, key conclusions are unchanged: a large weatherability increase is probable, the temperature sensitivity of continental weathering is weak, and Earth system climate sensitivity is likely to be higher than fast-feedback estimates. A high Earth system climate sensitivity is supported by numerous paleoclimate studies [44][45][46][47] .

Supplementary Note 6: Modified models -Different functional forms for continental weathering
The precise functional form for continental weathering is uncertain, and so we repeated our inverse analysis using five different functional forms for continental silicate and carbonate weathering to see how the choice of function affected the probability distributions for key parameters). This expression represents the effect of CO 2 fertilization on vascular land plants. The unknown exponent, , represents the efficiency of CO 2 fertilization. We vary across the full range from 0 to 1 in our inverse analysis such that it encompasses the endmember case where there is no direct pCO 2 dependence. For most of this range, the Michaelis-Menton law provides a weaker pCO 2 dependence than the original power law. For simplicity we adopt the same functional form for carbonate weathering as for silicate weathering, but recall the carbonate weathering function contains an additional free weatherability parameter to account for any differences with silicate weathering. Selected case 1 results are presented in the main text in Table 1 and Fig. 6.
Case 2: No direct pCO 2 dependence and linear runoff dependence. The direct pCO 2 dependence of silicate weathering adopted in the main text is uncertain and potentially overstated. Here, we omit any pCO 2 dependence and instead use the functional forms for silicate weathering adopted in field studies (e.g. 50,51 ) with an Arrhenius temperature dependence and linear runoff dependence. We relate runoff to temperature using an expression from the COPSE model 52 Here ranges from 0.2 to 0.5 as in the main text, and RUN ranges from 0.025 to 0.045 as above.
Case 5: No direct pCO 2 dependence and linear runoff dependence with unknown exponent. Rather than adopt the runoff expression adopted in the COPSE 52 model, we instead assume the runoff exponent is unknown and allow it to vary freely from 0 to 1. This allows for a wide range of runoff sensitivities to climate: Here =0-1 is the runoff exponent and RUN=0.025-0.045 is the runoff sensitivity. Finally, a smaller increase in weatherability since 100 Ma is more probable in cases 1, 2, and 5 than in the original results. This can also be understood in terms of a weaker or absent direct pCO 2 dependence: without a pCO 2 dependence decreasing continental weathering since 100 Ma, the counterbalancing increase in weatherability does not need to be as large in order to fit proxies.
In summary, the choice of functional form for continental weathering does not change any of our qualitative conclusions, namely climate sensitivity is large, the temperature dependence of continental weathering is weak, and there has been a sizeable increase in weatherability since 100 Ma due to factors other than climate. However, our quantitative estimates of key variablesparticularly the weatherability changecould be refined by a better mechanistic understanding of continental weathering on the global scale.

Supplementary Note 7: Silicate weathering and outgassing
From equation (2) Table 1 were used to compute the doubling temperatures reported in the abstract.
The outgassing reconstructions in Fig. 7 are sourced from Van Der Meer et al. 54 (blue), Vérard et al. 55 (red), Cogné and Humler 56 (aqua), Seton et al. 57 (purple), Hansen and Wallmann 58 (green), and Berner 48 (GEOCARB, yellow). These reconstructions were chosen because they rely on several independent lines of evidence: reconstructions of plate extent and plate motion based on field geology, paleogeography and paleomagentic data 55 , seismic imaging of subducted plates to infer variations in subducted plate length over time 54 , and reconstructions of seafloor age and depth 56,57 . As noted in the main text, the outlying Cogné and Humler 56 reconstruction is disputed 57,59 .

Ocean chemistry
The saturation state for calcite is calculated using the following expression from Pilson 60 Here, T is the temperature of the ocean or pore-space. A constant salinity of 35 parts per thousand is assumed in this expression. Note that we have not included a correction for changing Mg 2+ abundances (e.g. ref 3 ). Hain et al. 61 showed that including an Mg 2+ correction factor introduced significant error to sp K . This is because the empirical correction factor is normally offset by changes in * 2 K (defined in equation (S23)). The most accurate approach, short of implementing the Pitzer equations, is to use the standard, fixed thermodynamic constants 61 . We omit the pressure dependence of the solubility product because pressure is not varying through time (multiplication by a constant would not change model outputs since fluxes are scaled to fit modern fluxes).
For simplicity and computational efficiency we used a constant T=291.15 K in equation (S23), which is the modern surface ocean temperature (surface oceans are warmer than the continents plus oceans average of 285 K). We conducted test model runs with large changes in T and found that it doesn't change observable parameters (pH, CO 2 , carbon fluxes, saturation state) appreciably. Using temperature-dependent constants results in only moderate changes in alkalinity and carbonate speciation.
Derivation of equation (25) From the definition of carbon abundance (Tmol C kg -1 ) in the atmosphere-ocean box we have: This quadratic can be solved to find the + H   molality in either the ocean or the pore space.

Proxy records
This section describes the geochemical proxy data plotted in Fig. 2

Ocean pH
Ocean pH estimates were taken from Anagnostou et al. 4 , Tripati et al. 5 , Foster et al. 6 and Bartoli et al. 7 . Data from Edgar et al. 8 and Pearson et al. 9 recalculated using the methodology of Anagnostou et al. 4 were also used. Supplementary Fig. 10 shows all the data considered in this study and the 10 Ma binned data that were used in the forward modeling and inverse analyses described in the main text.
Ocean pH is typically inferred from boron isotope data. Boron speciation is pH dependent, and so ocean pH can be estimated from boron isotopes in carbonates, assuming the boron isotopic composition of seawater and the fractionation factor are known 62 . Anagnostou et al. 4 rigorously constrained Cenozoic ocean pH using a broad range of estimates for boron isotope vital effects and seawater isotopic composition. The resulting data form the basis of the pH dataset in this study. However, the uncertainty in these ocean pH estimates may still be underestimated if the computed values do not reflect a globally integrated signal. In the modern ocean, surface pH varies by at least 0.2 log units 63 .
Atmospheric CO 2 Supplementary Fig. 11 shows the CO 2 proxy data considered in this analysis and the 10 Ma binned data that were used in the forward modeling and inverse analyses described in the main text. Cenozoic pCO 2 proxies were taken from Beerling and Royer 10 . This exhaustive compilation includes stomatal, paleosol, phytoplankton, liverwort, and boron-based estimates and represents a consensus reconstruction.
Cretaceous pCO 2 proxies were taken from Hong and Lee 11 (pedogenic carbonates), Franks et al. 12 (stomatal model), Fletcher et al. 13 (fossil bryophytes), Retallack 14 (stomatal index), Quan et al. 15 (stomatal index) and Barclay et al. 16 (stomatal index). Traditionally, isotopic methods such as those using pedogenic carbonates yield higher pCO 2 estimates than those from stomatal methods 64 . This discrepancy was partially resolved by Breecker et al. 65 who revaluated soil CO 2 concentrations downward yielding lower atmospheric pCO 2 estimates from pedogenic carbonates. However, Cretaceous pCO 2 estimates based on the Breecker et al. 65 methodology 11 remain significantly higher than recent contemporaneous stomatal estimates 12 . These recent stomatal estimates are based on a mechanistic model of leaf exchange rather than an empirical fit to stomatal indices and are arguably more accurate at high pCO 2 12 . However, the sensitivity of this model-based approach to parameter assumptions was highlighted by McElwain et al. 66 , who concluded that Franks et al. 12 could have underestimated atmospheric pCO 2 . In short, the discrepancy between stomatal and paleosol pCO 2 estimates remains unresolved. Consequently, we include data from both methods in Supplementary Fig.  11 and use the full range of proxy estimates to define the uncertainty in each 10 Ma interval.

Temperature
Paleocene-Eocene Thermal Maximum (PETM) temperatures are described in Supplementary Note 1. Cretaceous temperatures were estimated as follows. Deep ocean temperatures were taken from Fig. 3 in Huber et al. 67 . The four timespans reported in Huber et al. 67 17 . We omit plotting Cenozoic surface temperatures for lack of a consensus time-series for globally averaged temperatures. Supplementary Fig. 12 summarizes all the temperature data used in this study, and shows the binned ranges used in the inverse analysis and forward modeling.
Saturation state Supplementary Fig. 13 shows the variety of CCD reconstructions considered in this study 3,[18][19][20][21][22] . The ranges of CCD reconstructions in each 10 Ma interval were used as the consensus CCD in this study. The CCD was converted to a saturation state using equation 4 in Jansen et al. 72 73 convert this to a carbon flux by assuming the modern crustal formation rate is 3 km 2 yr -1 , and then multiplying the resultant flux by a factor of 1.5 to correct for the lower crust contribution. Conveniently, the net effect of this conversion is equivalent to multiplying by (10 12 /1.5)×1.5, and so wt% CO 2 can be easily converted to Tmol C yr -1 by multiplication by 10 12 . Thus the initial value for the Cenozoic porespace precipitation flux is 0.45 Tmol C yr -1 ( Table 2). In contrast, Mills et al. 74 assumed a modern seafloor sink of 1.72 Tmol C yr -1 in their carbon cycle model. However, this flux was obtained from indiscriminately averaging the carbonate content in both Cenozoic and Mesozoic oceanic crust. Overestimating the modern seafloor weathering flux will inflate the importance of the seafloor weathering sink at earlier times in Earth's history.
The upper 300 m of Cretaceous-aged drill cores have CO 2 contents of 2.4±0.7, 1.9±0.3, which gives a best-estimate of 2.35±0.75 wt% using the range method adopted for the other proxies. The cretaceous seafloor precipitation flux will therefore be (2.35±0.75)×(1+V) β Tmol C yr -1 , where we have multiplied by the Cretaceous crustal production rate relative to modern. The error bars for the Cretaceous pore-space flux in Fig. 2, 3, 4, 5, and supplementary figures 2, 15, and 16 were obtained by sampling our assumed range for V and calculating the inferred range of fluxes. For the Bayesian analysis, the observed precipitation sink was calculated for each forward model call by multiplying by (1+V) β .

Isotopic constraints
Isotopic records provide potentially powerful constraints on carbon cycle processes e.g. ref 75 . However, we have deliberately avoided using Sr, Os, or Li isotope records to constrain our model due to uncertainties in their respective interpretations. For example, the upward trend in

86
Sr Sr over the Cenozoic has multiple interpretations (see main text), and the Li isotope record cannot be straightforwardly related to continental weathering fluxes [76][77][78] .

Bayesian inversion
The log-likelihood function used in the Bayesian inversion is defined as follows 79 :  are the observed (proxy) value, model value, and uncertainty in the observed value of the i -th data point of the k -th variable, respectively. The k summation is over the six variables for which we have proxy records: mean surface temperature, mean deep ocean temperature, atmospheric pCO 2 , ocean saturation state, ocean pH, and pore-space precipitation flux. The i summation is over the binned data plotted in the time series figures in the main text (Fig. 2, 3, 4, 5, and supplementary figures 2, 15, and 16), and described above.
The Bayesian analysis was implemented using the 'emcee' package in Python 80 and posterior distribution figures were created using the 'corner' module in Python. The emcee package implements an affine-invariant MCMC ensemble sampler. We used 1000 walkers and 10,000 model stepsthat is a total of 10 million forward model callsto build posterior distributions for our parameters. After accounting for autocorrelation in the walkers, the effective sample size was ~100,000. The initial walker positions were randomized, and a 1000 step burn-in was discarded. Approximately three quarters of walkers have crossed the median parameter value by 1000 steps, indicating that this burn-in is adequate. The 95% credible intervals plotted in Fig.  5 are the 2.5 th -97.5 th percentile range in the distribution of walkers.
A small number of studies have applied Bayesian methods to carbon cycle models. For example, Royer et al. 81 and Park and Royer 45 adopted GEOCARBSULF 82 as a forward model and pCO 2 proxies to constrain Earth system climate sensitivity. Because these studies used pCO 2 proxies for the entire Phanerozoic they could estimate climate sensitivity separately for glacial and non-glacial periods. However, their retrieval only used pCO 2 to constrain model parameters and did not make use of temperature, saturation state, pH, and seafloor carbonate proxies. Additionally, the retrieval in Park and Royer 45 is dependent on the detailed parameterizations within GEOCARBSULF. Although some parameters such as silicate weathering activation energy and biological modifiers on weatherability were allowed to vary in Park and Royer 45 , assumptions about the carbon and strontium isotope records and continental land area through time are 'hard-wired' into GEOCARBSULF. In contrast, we have deliberately designed our forward model to be as general as possible: isotope records were not used to force our model, and outgassing, carbonate and silicate weatherability, and modern-day fluxes are all free parameters. Consequently, the conclusions from our Bayesian analysis are more robust to model assumptions.