Model description

Python code for the model is available open source from the first author’s website. The time evolution of the carbon cycle is described by the following set of equations:

Here C O and C P are the concentrations of carbon (Tmol C kg−1) in the atmosphere–ocean and pore-space, respectively. The carbon concentration in the pore-space is equivalent to the dissolved inorganic carbon (DIC) abundance, C P =DIC P , whereas carbon in the atmosphere–ocean reservoir is equal to marine DIC plus atmospheric carbon, , where s is a scaling factor equal to the ratio of total number of moles in the atmosphere divided by the mass of the ocean, . Similarly, and are the carbonate alkalinities in the atmosphere–ocean and pore-space, respectively (Tmol eq kg−1). The global outgassing flux (Tmol C per year) is specified by F out , whereas the rates of continental silicate weathering and carbonate weathering are F sil and F carb , respectively (Tmol C per year). Seafloor weathering from basalt dissolution (Tmol eq per year) is F diss , and the precipitation flux of carbonates (Tmol C per year) in the ocean and pore-space are given by P ocean and P pore , respectively. The mass of the ocean and the pore space is given by and , respectively20. We assume a range for J of kg per year, which implies the time to circulate one ocean volume through the pore-space is between =20 and =1,000 kyr, consistent with estimates from Johnson and Pruis56 and Caldeira20.

A common simplification in carbon cycle modelling is to neglect carbonate weathering15,22,24. This is justified on the grounds that carbonate weathering does not constitute a net carbon source. On long timescales, the carbon consumption from silicate and seafloor weathering must balance carbon outgassing (plus any imbalance in organic weathering and burial), and this balance determines atmospheric CO 2 and climate. However, the saturation state of the ocean is affected by carbonate weathering, and so we include carbonate weathering to model ocean chemistry. We do not track crustal or mantle reservoirs, and the atmospheric reservoir of CO 2 is set by equilibrium partitioning with the ocean (see below). Because we are not tracking fluxes between the atmosphere and ocean, continental silicate weathering is not a carbon source or sink; the release of cations (alkalinity) from silicate dissolution does not directly add carbon to the combined atmosphere–ocean system. Instead, the cations consume carbon indirectly when they later precipitate as marine carbonates.

There is no surface ocean box in our model. This is because we are only interested in the changes in properties of the bulk ocean, which are determined by varying boundary conditions not by the deep-surface partitioning. In addition, we are focused on timescales >106 years, and so the dynamics of deep-surface ocean mixing are unimportant. If we were to partition the ocean into surface and deep, the main difference would be differential temperature dependence of carbon speciation. We also do not include organic carbon weathering and burial. This is justifiable because the negative feedback from oxidative weathering ensures they are approximately balanced57. Empirically, the carbon isotope record reveals that the organic burial fraction has changed by only ∼10% over the last 100 Myr ago58. Consequently, the change in the organic burial flux is likely modest and so the omission of organic carbon will not affect our conclusions.

The functional forms for the flux terms not already described in the main text are presented below.

Continental silicate weathering

The continental weathering parameterization was described in the main text except for the coefficient α in equation (2). This coefficient is assumed to be 0.2–0.5 (ref. 59). Strictly speaking, soil pCO 2 should replace atmospheric pCO 2 in equation (2) because soil pCO 2 determines soil pH, and therefore silicate dissolution rates. However, if modest changes in maximum biological productivity are allowed, then the range of CO 2 dependencies from the term , with α varying from 0.2 to 0.5, is broadly equivalent to replacing atmospheric pCO 2 with soil pCO 2 (ref. 60; see Supplementary Fig. 6). Consequently, we retain atmospheric pCO 2 and in equation (2). To test the sensitivity of our results to weathering parameterizations, in addition to the pCO 2 power-law dependence in equation (2), we also consider a Michaelis–Menten law in the main text (see results).

The range adopted for relative Cretaceous weatherability, 1+W, is based on literature estimates10,23,43,47,49 of how external factors may have affected weatherability, and extending the range in either direction does not markedly change our results. The weatherability factor, , can also be interpreted as the sensitivity of the weathering response to changes in pCO 2 and temperature. For example, an increase in implies that an increase in surface temperature will result in a larger change in the continental weathering flux, F sil .

Continental carbonate weathering

We assume that carbonate weathering has the same functional form as silicate weathering, except for an additional dimensionless multiplicative factor, ω carb , to allow for the possibility that carbonate weathering is subject to different temperature dependence, CO 2 dependence and weatherability factors:

The carbonate weatherability factor is defined as follows:

We assume a range of values for C WF from −0.9 to 1.5 to allow for large differences between carbonate weathering and silicate weathering. For example, carbonates may have a lower effective activation energy compared to silicates because carbonate weathering is sensitive to runoff, whereas silicate weathering is sensitive to both runoff and a kinetic temperature effect7,25. Carbonate weathering may also have a different response to changes in uplift23, or varying fluxes due to changes in the crustal reservoir of carbonates. It should be noted that changes in carbonate weathering only affect saturation state; the changes in CO 2 , temperature and ocean pH due to carbonate weathering changes are negligible. This is because—as explained above—temperature and CO 2 are set by the balance between outgassing and silicate weathering plus seafloor weathering. Consequently, any conclusions drawn about those variables are unaffected by our formulation for carbonate weathering.

Climate model

To relate ΔT S to pCO 2 , we adopt the following climate model:

Here ΔT 2x is the climate sensitivity parameter, ΔP is a palaeogeography parameter and the second term accounts for solar luminosity changes (Supplementary Note 2). We divide by ln(2) so that ΔT 2x has conventional units of Kelvin warming per CO 2 doubling. Supplementary Fig. 9 compares different climate parameterizations and GCM results from the literature and illustrates why equation (9) is suitable. The parameter ΔP is the secular cooling (in K) since the mid-Cretaceous due to palaeogeography changes. A review of GCM studies concluded that ΔP=0–3.0 K (ref. 52). We assume ΔP=0–5.0 K to be conservative, noting that some models suggest 5 K of warming from an Eocene continental configuration61.

Outgassing

Estimates of Cenozoic and Mesozoic outgassing histories vary substantially. Figure 7 shows a variety of outgassing reconstructions from the literature expressed as crustal production or spreading rates, which are assumed to covary with global outgassing. Generally, these reconstructions suggest that global outgassing at 100 Ma was between 1.5x and 2.5x modern outgassing. This conclusion is based on several independent lines of evidence including reconstructions of plate extent and plate motion, seismic imaging of subducted plates, and reconstructions of seafloor age and depth (see Supplementary Note 7 for a summary of outgassing estimates with references). The outlying reconstruction in Figure 7 (aqua curve) is disputed because it uses a contentious crustal age distribution42.

Figure 7: Global outgassing reconstructions. Reconstructions of crustal production rates or spreading rates, relative to modern. These histories are assumed to reflect global outgassing histories. The studies listed in the legend are given in Supplementary Note 7. The grey region is the range of outgassing histories explored in our model. Full size image

For simplicity, we assume a linear global outgassing history:

Here is modern outgassing (Tmol C per year), t is time (in Myr ago) and V is a dimensionless scaling factor. Using crustal production or spreading rate as a proxy for global outgassing is a simplification because it ignores subaerial metamorphism and hot spot volcanism (for example, ref. 62). Given the uncertainty in these other contributions, we adopt a very broad range of outgassing histories since 100 Myr ago by assuming V=0.2–1.5. Thus, we allow mid-Cretaceous outgassing to range from 20 to 150% greater than modern.

Basalt dissolution and seafloor weathering

The temperature dependence of seafloor weathering uses the following Arrhenius-style expression17:

Here E bas (kJ mol−1) is the effective activation energy of basalt dissolution, R is the universal gas constant and T pore is the pore-space temperature. Coogan and Dosso17 reported an empirically derived activation energy of E bas =92±7 kJ mol−1, whereas experimental studies14,31 of basalt dissolution suggest activation energies between 42 and 109 kJ mol−1. We adopt a range of activation energies from E bas =40 to110 kJ mol−1.

Because Cenozoic and Mesozoic pore-space temperatures are controlled by deep-ocean temperature17, we must determine the link between global mean surface temperatures and deep-ocean temperatures. Figure 8 shows mean global surface temperatures plotted against deep-ocean temperatures using output from fully coupled atmosphere–ocean GCMs and palaeoclimate proxy data (see Supplementary Note 1 for details). The relationship is described by an empirical linear fit:

Here T D (K) is the mean deep-ocean temperature and T S (K) is the mean surface temperature. The best-fit gradient and intercept are a grad =1.02 and b int =−16.7, respectively. However, we assume a broad gradient range a grad =0.8–1.4, whereas the intercept, b int =274.037−a grad × 285 is chosen to ensure consistency with modern conditions. Figure 8 shows the range of possible relationships used in this study.

Our parameterization of deep-ocean temperature improves upon Brady and Gíslason14 because ours is based on an ensemble of GCM results and globally averaged palaeoclimate proxies rather than a single climate model. In addition, Brady and Gíslason14 overestimated the dependence of deep-ocean temperature on surface climate because their parameterization is based on a single near-equatorial latitude of 6.7°. The relationship between globally averaged abyssal temperatures and surface climate is more muted than the relationship with equatorial abyssal temperatures. This can be seen in our Fig. 8 and in Fig. 12 of Manabe and Bryan63, despite the model of the latter being the basis for the Brady and Gíslason14 parameterization.

To relate the pore-space temperature to the deep-ocean temperature, we adopt empirical results17. Oxygen isotopes indicate that for both the Cenozoic and Mesozoic, the mean pore-space temperature of seafloor carbonate precipitation is consistently ∼9 K warmer than the minimum temperature of seafloor carbonate precipitation (deep-ocean temperatures). Consequently, we assume T pore =T D +9. This modification has a very minor effect on the model output because it is largely the change in temperature, and not its absolute value, that controls variations in the seafloor weathering flux.

Quantifying the pH dependence of seafloor weathering is more challenging because most experiments either fail to separate the pH and direct CO 2 effect14, focus on individual minerals rather than whole-rock dissolution rates64,65, or do not explore the full pH range relevant to seafloor weathering66.

Gudbrandsson et al.67 measure whole-rock crystalline basalt dissolution rates for 2<pH<11. They find that Ca release is a U-shape function of pH (their Fig. 8), where the minimum of the ‘U’ at 25 °C is somewhere between pH=7 and pH=9, depending on the assumptions made about the reactive surface area (experimental results are sparse and so are fitted with an analytic model of dissolution). Given the uncertainty in this dissolution curve, it is difficult to predict the sign of the dissolution change for a modest change in ocean pH. For example, a decline in pH from 8.2 to 7.4—which is approximately the change from the modern ocean to mid-Cretaceous—predicts a 20% decrease in Ca release according to one Gudbrandsson et al.67, their Fig. 8f fit, and a 7% increase in Ca release according to the alternative Gudbrandsson et al.67, their Fig. 8e fit. Either way, the change in dissolution is minor, and so a possible first order approximation is to assume basalt dissolution on the seafloor is independent of pH for the Mesozoic and Cenozoic.

However, the Gudbrandsson et al.67 experiments may not accurately capture the pH dependence of seafloor weathering because they were not done in seawater and did not include carbon chemistry. Some experimental studies65,68 show that olivine dissolution is CO 2 -dependent at high pH values, ostensibly because abundant carbonate ions protect Si–O bonds, thereby decreasing dissolution with increasing DIC. In contrast, Golubev et al.64 studied the effect of CO 2 on dissolution rates for a range of pH values and found that forsterite, diopside and hornblende do not have CO 2 -dependent dissolution rates at any pH. Unfortunately, however, Golubev et al.64 did not consider plagioclase, so their results cannot easily be extrapolated to the basaltic seafloor. Wolff-Boenisch et al.66 showed that the CO 2 dependence of crystalline basalt is independent of CO 2 at low pH levels, but did not repeat the experiment at high pH values.

We allow for the possibility of pH dependence by setting the rate of dissolution proportional to , where γ varies from 0 (no pH dependence) to 0.5 (strong pH dependence dominated by basaltic glass dissolution):

Here k diss is a proportionality constant chosen to match the modern flux, is the hydrogen ion molality in the pore-space and is the modern molality. Because dissolution is dependent on crustal production and crustal production is proportional to global outgassing, we assume dissolution is dependent on outgassing with some unknown power–law relationship, defined by β=0–1.

Better knowledge of the pH dependence, temperature dependence and CO 2 dependence of basalt dissolution would improve our constraints on the carbon cycle. Specifically, whole-rock dissolution experiments performed at high pH that separate the effects of pH and DIC would allow for more precise parameterizations of seafloor weathering.

Precipitation fluxes

The precipitation flux of marine carbonates, P ocean , is the sum of the fluxes of shelf carbonates, P shelf , and pelagic carbonates, P pelagic :

Following Ridgwell69, the shelf precipitation flux is given by:

Here Ω O is the saturation state of the ocean (defined below by equation (20)), A shelf is the area of continental shelf available for carbonate precipitation, with denoting the modern shelf area, and k shelf is a proportionality constant. Shelf area is approximated by a polynomial fit (Supplementary Fig. 7) to reconstructed tropical shelf area from Walker et al.70. The exponent n defines the proportionality between the saturation state of the ocean and the precipitation flux. This is typically taken to be 1.7 based on the latitudinal dependence of carbonate accumulation and saturation state71. Rather than consider calcite and aragonite precipitation separately, we instead allow n to vary widely from 1.0 to 2.5.

The pelagic carbonate flux depends on the calcite compensation depth (CCD), Z CCD (km), which can be calculated using the following equation72:

Pelagic carbonate deposition is proportional to the fractional area above the CCD, f(Z CCD ), which can be approximated by an exponential fit to hypsometric data73:

Equations (14, 15, 16, 17) can then be combined to give the total ocean precipitation flux:

The proportionality constants k shelf and k pelagic are chosen to reproduce the modern partitioning between shelf and pelagic carbonates (see below). Here the precise functional form of equation (18) only matters for determining ocean saturation state; carbon fluxes, pCO 2 and temperatures are unaffected.

The pore-space carbonate precipitation flux is analogous to shelf precipitation except that there is no area dependence:

Here Ω P is the saturation state of the pore-space and k pore is a proportionality constant chosen to reproduce the modern flux. The exponent n is the same as for shelf precipitation. Repeating the inverse analysis allowing different exponents for shelf and seafloor carbonate precipitation does not change results substantially (not shown).

Finally, the saturation state of the ocean and the pore-space are defined as follows:

Here K sp =K sp (T) is the temperature-dependent solubility product from Pilson74. Supplementary Methods explain how the solubility product is calculated.

Ocean chemistry

Alkalinity and DIC have the following standard definitions in our model, where ALK is often referred to as ‘carbonate alkalinity’ in the literature:

Given carbon and alkalinity in the atmosphere–ocean (C O , ALK O ) or the pore-space (C P , ALK P ), we can calculate ocean chemistry using the following set of equations74:

Here H CO2 is the Henry’s law constant for CO 2 , [CO 2 aq] is the sum of the concentrations of free CO 2 and H 2 CO 3 , and and are the first and second apparent dissociation constants of carbonic acid, respectively. Temperature-dependent expressions for these constants can be found in Supplementary Methods. The set of equations described above must be solved separately for the ocean and the pore-space by substituting the generic carbon concentration and alkalinity (C, ALK) for (C O , ALK O ) and (C P , ALK P ), respectively. The scaling factor, s, is defined above with respect to equation (6). Equation (25) is derived by combining equations (21), (23) and (24) (Supplementary Methods). This quadratic can be solved to find [H+]. Once this is known, then equations (22), (23) and (24) define the remaining carbon chemistry variables.

Rather than attempt to model the complexities of calcium and magnesium cycling in our model, we impose observed changes in [Ca2+] abundances from seawater inclusions75. In our model, changes in alkalinity are driven purely by weathering and carbonate precipitation. Thus, by imposing [Ca2+] variations, we are effectively assuming that observed [Ca2+] changes are offset by changes in other cations and anions, such that they have no direct effect on alkalinity, for example, magnesium exchange with the seafloor. We fit a third-order polynomial to the [Ca2+] reconstruction in Tyrrell and Zeebe75 to achieve the fast computation times necessary for Bayesian inversion (Supplementary Fig. 8).

Initial conditions and numerical solution

Table 2 shows all the initial values assumed in our model or ranges for variables that are uncertain. All other initial values are fully determined by the variables in this table.

The system of differential equations describing the carbon cycle (equation (6)) was solved in Python using the ordinary differential equation integrator in the SciPy module. Model outputs were compared with equivalent steady-state calculations and were always in agreement to within a few per cent or better (Supplementary Note 3). This validates the numerical integration and implies that the time-dependent model is always in quasi steady state. In addition, the integrated flux imbalance over 100 Myr ago equals the change in the carbon reservoirs to within ∼2% or better in every case, confirming that mass is being conserved in our model.

Proxies

The geochemical proxies plotted in Figs 2, 3, 4, 5, and Supplementary Figs 2, 15 and 16 are described in Supplementary Methods. For each variable, we searched the literature for the broadest possible range of proxy estimates. Proxies were typically binned into 10 Myr ago intervals, and for each interval the best estimate was taken to be the midpoint of the full range of proxy estimates, while the 1σ uncertainty in the best estimate spanned the full range (Supplementary Figs 10–13). This conservative approach helps ensure that our conclusions are robust to uncertainties in different proxy methods.

Data availability

The binned proxy data used as inputs for this analysis along with the carbon cycle model code are available on the website of the first author. www.krisstott.com