Calcium silicate

The geological carbon cycle model builds on that described in Krissansen-Totton and Catling (36 (link)). Here, we summarize its key features, and additional details are provided in the SI Appendix. The Python source code is available on GitHub at github.com/joshuakt/early-earth-carbon-cycle.
We model the time evolution of the carbon cycle using two separate boxes representing the atmosphere–ocean system and the pore space in the seafloor (Fig. 1 and SI Appendix A). We track carbon and carbonate alkalinity fluxes into and between these boxes, and assume that the bulk ocean is in equilibrium with the atmosphere.
Many of the parameters in our model are uncertain, and so we adopt a range of values (SI Appendix, Table S1) based on spread in the literature rather than point estimates. Each parameter range was sampled uniformly, and the forward model was run 10,000 times to build distributions for model outputs such as pCO₂, pH, and temperature. Model outputs are compared with proxy data for pCO₂, temperature, and carbonate precipitation (SI Appendix D).
Continental silicate weathering is described by the following function: $F_{\text{sil}}=f_{bio}{f}_{\mathrm{land}}{F}_{\text{sil}}^{\mathrm{mod}}{\left(\frac{{\text{pCO}}_2}{{\text{pCO}}_2^{\mathrm{mod}}}\right)}^{\alpha }\exp \left(\mathrm{\Delta }{T}_{\text{S}}/T_e\right)$
Here, $f_{\text{bio}}$ is the biological enhancement of weathering (see below), $f_{\text{land}}$ is the continental land fraction relative to modern, $F_{\text{sil}}^{\mathrm{mod}}$ is the modern continental silicate weathering flux (Tmol y⁻¹), $\mathrm{\Delta }{T}_{\text{S}}=T_{\text{S}}-T_{\text{S}}^{\mathrm{mod}}$ is the difference in global mean surface temperature, $T_{\text{S}}$ , relative to preindustrial modern, $T_{\text{S}}^{\mathrm{mod}}$ . The exponent $\alpha$ is an empirical constant that determines the dependence of weathering on the partial pressure of carbon dioxide relative to modern, ${\text{pCO}}_2/{\text{pCO}}_2^{\mathrm{mod}}$ . An e-folding temperature, $T_e$ , defines the temperature dependence of weathering. A similar expression for carbonate weathering is described in SI Appendix A.
The land fraction, $f_{\text{land}}$ , and biological modifier, $f_{bio}$ , account for the growth of continents and the biological enhancement of continental weathering, respectively. We adopt a broad range of continental growth curves that encompasses literature estimates (Fig. 2A and SI Appendix A). For our nominal model, we assume Archean land fraction was anywhere between 10% and 75% of modern land fraction (Fig. 2A), but we also consider a no-land Archean endmember (Fig. 2B).
To account for the possible biological enhancement of weathering in the Phanerozoic due to vascular land plants, lichens, bryophytes, and ectomycorrhizal fungi, we adopt a broad range of histories for the biological enhancement of weathering, $f_{bio}$ (Fig. 2C). The lower end of this range is consistent with estimates of biotic enhancement of weathering from the literature (37 –39 (no links found)).
The dissolution of basalt in the seafloor is dependent on the spreading rate, pore-space pH, and pore-space temperature (SI Appendix A). This formulation is based on the validated parameterization in ref. 36 (link). Pore-space temperatures are a function of climate and geothermal heat flow. Empirical data and fully coupled global climate models reveal a linear relationship between deep ocean temperature and surface climate (36 (link)). Equations relating pore-space temperature, deep ocean temperature, and sediment thickness are provided in SI Appendix A.
Carbon leaves the atmosphere–ocean system through carbonate precipitation in the ocean and pore space of the oceanic crust. At each time step, the carbon abundances and alkalinities are used to calculate the carbon speciation, atmospheric pCO₂, and saturation state assuming chemical equilibrium. Saturation states are then used to calculate carbonate precipitation fluxes (SI Appendix A). We allow calcium (Ca) abundance to evolve with alkalinity, effectively assuming no processes are affecting Ca abundances other than carbonate and silicate weathering, seafloor dissolution, and carbonate precipitation. The consequences of this simplification are explored in the sensitivity analysis in SI Appendix C. We do not track organic carbon burial because organic burial only constitutes 10–30% of total carbon burial for the vast majority of Earth history (40 ), and so the inorganic carbon cycle is the primary control.
The treatment of tectonic and interior processes is important for specifying outgassing and subduction flux histories. We avoid tracking crustal and mantle reservoirs because explicitly parameterizing how outgassing fluxes relate to crustal production and reservoirs assumes modern-style plate tectonics has operated throughout Earth history (e.g., ref. 12 ) and might not be valid. Evidence exists for Archean subduction in eclogitic diamonds (41 (link)) and sulfur mass-independent fractionation in ocean island basalts ostensibly derived from recycled Archean crust (42 (link)). However, other tectonic modes have been proposed for the early Earth such as heat-pipe volcanism (43 (link)), delamination and shallow convection (44 (link)), or a stagnant lid regime (45 ).
Our generalized parameterizations for heat flow, spreading rates, and outgassing histories are described in SI Appendix A. Fig. 2D shows our assumed range of internal heat flow histories compared with estimates from the literature. Spreading rate is connected to crustal production via a power law, which spans endmember cases (SI Appendix A). These parameterizations provide an extremely broad range of heat flow, outgassing, and crustal production histories, and do not assume a fixed coupling between these variables.
We used a 1D radiative convective model (46 ) to create a grid of mean surface temperatures as a function of solar luminosity and pCO₂. The grid of temperature outputs was fitted with a 2D polynomial (SI Appendix E). We initially neglect other greenhouse gases besides CO₂ and H₂O, albedo changes, and assumed a constant total pressure over Earth history. However, later we consider these influences, such as including methane (CH₄) in the Precambrian. The evolution of solar luminosity is conventionally parameterized (47 ).
Our model has been demonstrated for the last 100 Ma against abundant proxy data (36 (link)) and it can broadly reproduce Sleep and Zahnle (12 ) if we replace our kinetic formulation of seafloor weathering with their simpler CO₂-dependent expression (SI Appendix B). Agreement with ref. 12 confirms that the omission of crustal and mantle reservoirs does not affect our conclusions.

Sediment profiling was performed using commercial microsensors operated with a motorized micromanipulator (Unisense A.S., Aarhus, Denmark). Sediment was brought level to the core liner surface. Cores were placed in an aquarium containing water collected from the study site and held at in situ temperature and constantly bubbled with air. During microsensor profiling, an airstream was additionally provided at the water surface upstream of the sediment cores to ensure constant water flow over the sediment, estimated at 2–4 cm s⁻¹ at 0.5 cm above the sediment surface. Oxygen and H₂S microelectrodes had a tip diameter of 50 μm, and pH microsensors had a tip diameter of 100 or 200 μm. Vertical oxygen microprofiles were recorded at 100 μm steps, beginning at least 1 mm above the sediment–water interface. A two-point O₂ calibration was made in air-saturated seawater and the anoxic zone of the sediment (100% and 0%, respectively). Vertical H₂S and pH microprofiles were made concurrently, recorded at 200 μm steps in the oxic zone, and 200 or 400 μm steps below. For H₂S, a five-point standard curve was made using Na₂S standards. ΣH₂S was calculated from H₂S based on pH measured at the same depth using the R package AquaEnv (Hofmann et al., 2010 ). For pH, a three-point calibration using NBS standard buffers was made, followed by a salt correction with TRIS buffer (Dickson et al., 2007 ). pH measurements were performed with a Ag/AgCl reference electrode and values are reported on the total scale. Porosity was determined from water content and dry density measurements upon drying to constant weight at 60 °C.
Reactive transport modelling was used to calculate reaction rates from the microsensor data (Table 1). The diffusive oxygen uptake (DOU) of the sediment was calculated from the oxygen depth profile near the sediment–water interface using Fick's first law as below:

where ϕ is the measured porosity, the term (1−2ln ϕ) is a correction for sediment tortuosity (Boudreau, 1996 ), and D_i is the diffusion coefficient for oxygen calculated at salinity (S) and temperature (T) calculated using the R package marelac (Soetaert et al., 2012 ). The maximum concentration gradient, ∂C_i/∂x, was calculated as the linear slope immediately below the sediment–water interface.
The cathodic proton consumption (CPC) is the proton consumption corresponding to the cathodic half-reaction, O₂+4e⁻+4H⁺→2H₂O. The CPC was calculated as the sum of the upward and downward alkalinity fluxes immediately above and below the pH maximum in the oxic zone. Alkalinity fluxes were calculated as the sum of the fluxes of the major individual porewater compounds that contribute to the alkalinity:

The speciation of individual ions (carbonate, borate, sulphide and water equilibria) were calculated from H₂S and pH microsensor data using the R package AquaEnv (Hoffman et al., 2010 ). For this calculation, we assumed a constant depth profile for dissolved inorganic carbon and borate. The contribution of other macronutrients (ammonium, silicate) to alkalinity fluxes are likely small and were not considered in this study.
The concentration gradient ∂C_i/∂x for each individual species was estimated by a noise robust numerical differentiation using the Savitsky–Golay algorithm, optimized with a N=21 point stencil (Savitsky and Golay, 1964 ). The current density was calculated from the CPC using a conversion factor of one electron per proton (1A=1.036 × 10⁻⁵ mol e⁻ s⁻¹). The percentage, p, of oxygen consumed electrochemically was calculated as P=25 (CPC/DOU), based on a reaction stoichiometry that four moles of H⁺ are consumed per mole of O₂ by the cathodic reaction. The remaining part of the oxygen consumption (1−p) is attributed to aerobic respiration, which does not consume or release alkalinity. This simplified model assumes that no other biogeochemical reactions are releasing or consuming alkalinity within the oxic zone. In reality, some alkalinity consumption will take place, associated with nitrification, ferrous iron oxidation, calcium carbonate precipitation or direct sulphur oxidation by oxygen (Risgaard-Petersen et al., 2012 ). The rates of these reactions cannot be estimated from microsensor profiles. However, correcting for these reactions would increase the value of p. Accordingly, the p estimates derived here should be interpreted as conservative, lower bounds on the oxygen that is consumed electrochemically.