Modeling the Geological Carbon Cycle

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.