A coupled carbon-silicon cycle model over Earth history: Reverse weathering as a possible explanation of a warm mid-Proterozoic climate

Article history: Received 13 August 2019 Received in revised form 20 February 2020 Accepted 24 February 2020 Available online xxxx Editor: L. Derry


Introduction
The carbon cycle controls Earth's climate evolution on long timescales. Silicate weathering is a temperature dependent net sink for atmospheric CO 2 , and thus provides a natural thermostat to buffer Earth's climate against the secular increase in solar luminosity (Walker et al., 1981). However, this thermostat may be modulated by "reverse weathering" reactions that produce clays on the seafloor, which consume ocean carbonate alkalinity generated by silicate weathering, but unlike carbonates do not consume atmospheric carbon (Garrels, 1965;Isson et al., 2019;Mackenzie and Kump, 1995). In this study, we investigate whether changes in the reverse weathering flux over The dissolved cation and alkalinity are ultimately deposited as carbonate on the seafloor according to the following reaction: The net result of equations (1) and (2)  Weathering of the continents and seafloor (equation (1)) is temperature dependent (Coogan and Dosso, 2015;Krissansen-Totton and Catling, 2017), and so acts as planetary thermostat. Numerous studies have applied carbon cycle models incorporating this continental weathering buffer-and sometimes with seafloor weathering also-to predict Earth's climate evolution over time Sleep and Zahnle, 2001;Tajika and Matsui, 1992). Earth's carbon cycle evolution also controls ocean pH evolution because pH is set by atmospheric CO 2 and the inventory of dissolved weathering products Tajika and Matsui, 1992). The dissolved silicic acid produced by weathering may ultimately precipitate as silica: This silica precipitation could be either biogenic or abiotic. Siever (1992) outlined the classic view of the Precambrian silicon cycle whereby all silica precipitation is assumed abiotic because biological silica precipitation had not yet evolved. Moreover, none of the abiotic silica precipitation was primary, but rather formed diagenetically via silicification of carbonate rocks, absorption by clay minerals, or evaporitic precipitation of silica. Ocean dissolved silica was assumed to be at silica saturation, or approximately 1500 μM (Siever, 1992).
Recent evidence, reviewed by , suggests that marine dissolved silica concentrations may have been below amorphous silica saturation for part of the Precambrian. Cyanobacteria such as Synechococcus precipitate substantial intracellular silica, possessing Si:P ratios comparable to that of diatoms (Baines et al., 2012). This could explain why Archean cherts formed from primary Si grains that precipitated in the water column (Stefurak et al., 2014). There is also compelling phylogenetic evidence for active Si transporters emerging in eukaryotes sometime between the Mesoproterozoic and the Ediacaran . With that said, the extent to which biology modulated the Precambrian silica cycle is unknown; there is no clear evidence from the fossil record to support significant Precambrian biogenic silicification. Moreover, chert is absent from organic rich shales in the mid-Proterozoic, contrary to what we would expect if silica precipitation from cyanobacteria-the ostensible primary producers-was large.
The silica cycle may affect the carbon cycle, climate, and ocean pH via reverse weathering reactions (Garrels, 1965;Isson et al., 2019;Mackenzie and Garrels, 1966;Mackenzie and Kump, 1995). Changes in clay-forming reactions could explain the evolution of seawater chemistry during the Cenozoic (Dunlea et al., 2017;Higgins and Schrag, 2015). Reverse weathering refers to a variety of clay-forming reactions that consume alkalinity without consuming carbon, thereby retaining CO 2 in the atmosphere-ocean system: The two alternative equations denote that there are a range of reverse weathering reactions, some of which alter cation-poor silicates, while others form silicates from exclusively aqueous reactants. The effect of reverse weathering can be understood as follows. In the absence of equations (5) or (6), cations and bicarbonate generated by continental weathering (equation (1)) will eventually form carbonate-bearing rock, thereby removing alkalinity and carbon from the ocean (2). However, if instead reverse weathering reactions consume alkalinity, then carbon remains in the atmosphere-ocean system. Intuitively, if the reverse weathering flux is large, then a greater flux of forward weathering (and therefore higher surface temperatures) are required to ensure carbonate burial matches carbon inputs from volcanism. Considerable uncertainty exists over the magnitude of the modern reverse weathering flux.  estimated the global reverse weathering sink to be 1.5 Tmol Si/yr, whereas Laruelle et al. (2009) and Tréguer et al. (1995) provide lower estimates of 1.0 Tmol Si/yr and 0.6 Tmol Si/yr, respectively. All these estimates are based, in part, on laboratory leaches of sediments to isolate biogenic silica and clays. However, it has recently been argued that traditional leaches underestimate authigenic silicates and that cosmogenic 32 Si abundances provide a more accurate measure of diagenetic Si deposition (Rahman et al., 2016;. In fact, Rahman et al. (2017) calculated 4.5-4.9 Tmol Si/yr for the modern reverse weathering sink based on 32 Si measurements in representative marine sediments. Rahman et al. (2019) revaluated Si inputs, namely submarine groundwater discharge, and concluded that a balanced Si cycle potentially requires an even larger reverse weathering sink. It should be noted, however, that these large reverse weathering estimates are based on extrapolations from a limited number of sites.
Even if reverse weathering is not a major silica sink today, in the Precambrian before the advent of ecologically significant biogenic silica precipitation, ocean silica concentrations were likely elevated Siever, 1992).  argued that these elevated silica levels may have increased the reverse weathering flux, maintained high levels of CO 2 , and warmed the Precambrian climate. A sediment diagenesis model coupled to a carbon cycle model was used to show high Precambrian pCO 2 , thereby helping to explain the apparent lack of glacial deposits in the mid Proterozoic (Kasting, 2005). Additionally,  proposed that reverse weathering is strongly pH dependent, and thus may have buffered ocean pH and CO 2 during the Precambrian.
Here, we investigate the role of reverse weathering on Earth's climate evolution by applying the first fully coupled carbon-silicon cycle model across all of Earth's history. We believe this is an improvement over imposing dissolved silica concentrations through time (e.g.  because it explicitly and selfconsistently balances carbon fluxes, alkalinity fluxes, and silicon fluxes. We also considered a range of initial (modern) conditions for modern reverse weathering fluxes, reflecting literature estimates. As noted above, recent papers have argued that the modern reverse weathering flux is considerably larger than traditional estimates. If true, this makes unlikely a large change in reverse weathering-and therefore a change in climate-across the Precambrian/Phanerozoic boundary.

Carbon-silica cycle model
We coupled the carbon cycle model from  to a new model of silicon cycling. The updated Python code is available on the lead author's Github upon publication. The biogenic silica sink flux, F bio (mol Si/yr), and authigenic clay reverse weathering sink flux, F RW (mol Si/yr), are parameterized using a sediment diagenesis model described below. The abiotic silica sink flux, F Abio (mol Si/yr), is related to dissolved silica concentrations in the ocean, whereas dissolved silica sources are dependent on continental weathering flux, F sil (mol C/yr), and seafloor weathering flux, F diss (mol C/yr). The evolution of dissolved silica in the ocean is determined by the following equation: Here, R Sil:Alk(FW) is the stoichiometric ratio of silicon to alkalinity release by continental and seafloor (forward) weathering reactions, which is assumed to be 0.5 (equation (1)). A factor of 2 is required to convert the equivalent carbon weathering fluxes to alkalinity fluxes. [Si] is the concentration of dissolved silica in the ocean (mol/kg) and M O is the mass of the ocean (kg).
The abiotic silica sink, F Abio , is parameterized using a power law dependent on silica saturation, as described in Appendix B. The only unknown parameter in the equation describing abiotic silica precipitation flux is the reaction coefficient, which was chosen to reproduce estimated Precambrian [Si] levels . However, our results are relatively insensitive to the choice of this coefficient (not shown). Note that estimated Precambrian [Si] levels  are potentially below amorphous silica saturation (Siever, 1992) because (i) Fe-Si gels, which likely precipitated in Precambrian oceans, have Si solubility below that of amorphous silica, and (ii) the Si homeostasis mechanisms of some cyanobacteria such as Synechococcus results in intracellular silica precipitation that may have decreased Precambrian [Si], although note the counterarguments to Precambrian biogenic silicification described above. In any case, whether the reaction coefficient is chosen to reproduce saturated or slightly undersaturated Precambrian dissolved silica does not strongly affect our results (the error envelope of model outputs tends to overlap with saturated [Si] levels).
The silica cycle is coupled to the carbon cycle model by adding a reverse weathering term to the ocean alkalinity budget: Here, R Alk:Si(R W ) is the ratio of alkalinity to silicon consumed by reverse weathering reactions. A range of 2 to 4  is uniformly sampled in our carbon cycle calculations. Other terms in the ocean alkalinity budget are the same as in : A O is ocean alkalinity (mol eq kg −1 ), A P is pore space alkalinity (mol eq kg −1 ), J is a mixing term with the pore space (kg yr −1 ), F carb is the carbonate weathering flux (mol C yr −1 ), and P ocean is carbonate precipitation in the ocean (mol C yr −1 ). Appendix B contains the complete set of equations governing carbon, alkalinity, and silica fluxes in the model. Unless stated otherwise, all nominal model outputs assumed negligible atmospheric methane throughout Earth history to highlight the dynamics of CO 2 cycling and reverse weathering. However, planetary redox budget considerations (Catling et al., 2001) and recent direct inferences of Archean rates of hydrogen escape to space from Xe isotopes  coupled with geologic evidence for low levels of tropospheric H 2 (Kadoya and Catling, 2019), all suggest high levels (∼0.5%) of methane in the Archean. We thus also present carbon-silica cycle outputs in supplementary materials with a modified climate model where 100 ppm and 1% methane are imposed in the Proterozoic and Archean, respectively .

Sediment diagenesis model
A sediment diagenesis model is used to calculate biogenic silica burial and reverse weathering fluxes. This model is not intended to be an exhaustive treatment of diagenetic processes in marine sediments. Instead, the model is presented as a plausible, self-consistent set of calculations that illustrates a wide range of outcomes due to uncertainties in key parameters.
We implemented a three-component reaction-diffusion model that tracks dissolved silica, biogenic silica, and authigenic clays as a function of depth in sediments. A full description of the sediment diagenesis model is given in the supplementary materials (Appendix C). Here, we provide the key equation governing authigenic clays formed via reverse weathering. The rate of reverse weathering (μmol/cm 3 /s) is described by: where A clay > 50 wt.% of sediments (9) Here, x (cm) is sediment depth, C Si (μmol/cm 3 ) is dissolved silica concentration, A clay (in wt.%) is the concentration of reverse weathering clays in sediments. The sediment model is illustrated in Fig. 1. The three variables that have the strongest influence on the reverse weathering flux over Earth history are: 1) k R W (μmol/cm 3 /s), the reaction coefficient for reverse weathering reactions, 2) z A (cm), the e-folding depth governing the decrease in reverse weathering with depth, 3) S A (μmol/cm 3 ), the solubility of authigenic silicates.
These three variables are uncertain and may vary spatially due to differences in cation availability, sedimentation rates, and sediment porewater chemistry.
The literature has various estimates for the rate coefficient for reverse weathering reactions.  fit a sediment diagenesis model to sediment core data from the Peruvian shelf to obtain a reverse weathering coefficient of k R W = 3.5 × 10 −6 μmol/cm 3 /s (110 μmol/cm 3 /yr).  assumed k R W is pH-dependent and ranged from approximately 5 × 10 −7 μmol/cm 3 /s (pH = 8.0) to 5 × 10 −9 μmol/cm 3 /s Table 1 Choices for key parameters controlling reverse weathering fluxes. Bold are fitted to modern conditions using inverse methods, unbolded were chosen heuristically to fit modern conditions. Here, T&D MCs refers to modern conditions matching modern silica cycle flux estimates from , whereas Rahman MCs refers to modern conditions matching modern silica cycle flux estimates from . Remaining parameters are shown in Table S1 in Appendix E.

Scenario
Reverse weathering reaction coefficient, k RW (μmol/cm 3 /s) Plausible range: 4 × 10 −6 to 10 −10 (pH = 6.5). For comparison, silicon uptake rates in Amazon delta sediments are around 1.3 × 10 −6 μmol/cm 3 /s (Michalopoulos and Aller, 2004), which suggest an effective k R W around this order of magnitude for dissolved silica levels comparable to authigenic clay saturation. In contrast, silica consumption rates in sediment cores from the Sea of Okhotsk are around 10 −10 to 10 −11 μmol/cm 3 /s (Wallmann et al., 2008), therefore implying much smaller values for k R W . When fitting k R W we therefore assumed a log uniform prior from 3.8 × 10 −6 to 10 −10 μmol/cm 3 /s.  assumed that k R W is strongly pH dependent, with a pH 22 rate law. This is based on a polynomial fit to greenalite nucleation experiments in Tosca et al. (2016, their Fig. 2b), where a strong pH dependence is attributable to the exponential dependence of nucleation on saturation state (Rasmussen et al., 2019). Whether this saturation state dependence reflects the overall pH dependence of reverse weathering depends on whether reverse weathering is dominated by nucleation or crystal growth dominated. In any case, greenalite precipitation drops to zero at pH values <∼7.0. However, sediments in alkaline lakes such as Lake Turkana in Kenya (pH∼9.2) undergo extensive reverse weathering as evidenced by alkalinity consumption not associated with carbonate formation (Yuretich and Cerling, 1983). Thus, we explore the addition of a strong pH dependence of reverse weathering in our analysis but find that it is less important than changing biogenic silica precipitation and dissolved silica concentrations in the ocean (Appendix A). The solubility of authigenic silicates, S A , varies with clay composition and porewater chemistry. For authigenic clays formed by opal alteration, solubility estimates range from 0.22 μmol/cm 3 (Loucaides et al., 2010) to 0.33 μmol/cm 3 (Dixit et al., 2001). In Amazon delta sediments, dissolved silica asymptotes to around 0.5 μmol/cm 3 , which suggests a comparable value for authigenic clay solubility (Michalopoulos and Aller, 2004).  estimate S A using pH-dependent greenalite solubility derived from Tosca et al. (2016). For plausible Precambrian Fe 2+ con-centrations, S A (greenalite) ranges from 0.01 μmol/cm 3 (pH = 8) to 2.0 μmol/cm 3 (pH = 7.25) (Tosca et al., 2016, Fig. 13). Greenalite does not precipitate at pH < 7.25 because its solubility is greater than amorphous silica solubility. We adopt a uniform prior for S A (average for all authigenic clays) ranging from 0.2 μmol/cm 3 to 1.0 μmol/cm 3 (0.2-0.1 mM) when fitting S A .
There is also uncertainty in the e-folding depth of reverse weathering, z A . Reverse weathering reactions may decline with depth in sediments due to consumption of cations, especially aluminum , porosity decrease, or pH decrease.  derive a value of z A = 1.4 cm from fitting a sediment diagenesis model to sediment core data. However, Wallmann et al. (2008) reported reverse weathering increasing across the top 1-2 m of sediments, which would imply z A values that are very large (effectively no decay in reaction rate with depth). We adopt a log uniform prior for z A from 1 cm to 1000 cm when fitting z A .
The model ocean is partitioned into the proximal continental shelf (∼0.5% seafloor area), distal continental shelf (∼7.5% seafloor area) and open ocean (∼92% seafloor area) to allow for different sedimentation rates and kinetics in these different settings (see supplementary materials for details).
Sediment diagenesis parameters such as dissolution reaction coefficients, the exponential decay of biogenic silica dissolution due to opal maturation and porosity decline, the exponential decay of authigenic clay formation due to cation availability and porosity decline, and the mean biogenic flux are not known exactly, but instead have plausible ranges. Point estimates were chosen to reproduce modern conditions such as the correct biogenic silica burial fluxes, reverse weathering fluxes, approximate [Si] and biogenic silica concentrations at depth, and the ratio of biogenic sedimentation to permanent silica burial. This was either done heuristically or with a Markov Chain Monte Carlo (MCMC) retrieval depending on the case considered. Appendix E explains these estimates and procedures in detail, whereas Table 1 summarizes the Fig. 2. Carbon cycle and climate evolution over Earth history with imposed reverse weathering, which is an illustrative case only because the silica cycle is not selfconsistent. The solid blue and dashed blue lines show the median and 95% confidence intervals, respectively, for our nominal carbon cycle evolution assuming negligible reverse weathering throughout Earth history . The solid black line and gray shaded regions show the median and 95% confidence intervals, respectively, from the same carbon cycle model modified to include a reverse weathering flux as a fraction of the total dissolved silica sink (C), imposed by adopting the median f rw curve in . Elevated reverse weathering in the Proterozoic leads to somewhat higher atmospheric CO 2 (B) and temperature (D) envelopes. Subfigure (E) shows how the continental silicate weathering flux (black line and gray shaded) is higher in the Precambrian to balance the carbon cycle under elevated reverse weathering (red shaded). Individual data points represent proxy data, which are described in detail in  reverse weathering parameter choices for the analyses presented in the main text. The sediment diagenesis model produces (i) reverse weathering as a function of [Si] and biogenic silica precipitation as a function of [Si] for cases with biogenic silica deposition (Phanerozoic) and (ii) reverse weathering as a function of [Si] for no biogenic silica deposition. We fit 1D polynomials to these diagenesis model outputs to avoid calling the sediment diagenesis model in carbon cycle calculations (Appendix C). Since we are focusing on changes between the Proterozoic and Phanerozoic, we do not consider the effect of biological innovations during the Phanerozoic. Instead, we adopt parameterization (i) to determine reverse weathering and biogenic silica fluxes during the Phanerozoic, and smoothly transition to function (ii) over the preceding 0.5 Ga (Appendix C). Function (ii) governs reverse weathering for the remaining Precambrian. This transition was chosen to reflect recent estimates of the timing of the emergence of biogenic silica precipitation in the Neoproterozoic , although there is still considerable uncertainty in this timing (see Section 1).

Results
The first results (Fig. 2) are illustrative and did not use the self-consistent, coupled C-Si cycle model described above. Instead, we imposed the reverse weathering flux through time using the average reverse weathering fraction ( f r w ) curve derived in  and incorporated this into our carbon cycle model . Here, f r w is the fraction of total Si consumed through reverse weathering that  calculated from their sediment diagenesis model and an assumed seawater silica concentrations curve through time. The results in Fig. 2 show a comparison to carbon cycle outputs where reverse weathering is assumed to be negligible throughout Earth history. Consistent with the carbon cycle calculations of Isson and Planavsky (2018), we find increased reverse weathering in the Precambrian results in higher CO 2 and temperatures than in the negligible reverse weathering case. However, Earth's overall climate evolution is not dramatically changed by the inclusion of reverse weathering in this scenario. This is due to (i) the large uncertainties in other carbon cycle parameterizations drowning out the reverse weathering transition at 0.5-1.0 Ga, and (ii) our carbon cycle model being intrinsically less sensitive to reverse weathering changes than that of Isson and Planavsky (2018) (see Appendix F for detailed comparative calculations).
Next, we show the results of our self-consistent, coupled C-Si cycle model where reverse weathering and biogenic silica deposition are parameterized using our sediment diagenesis model. These model outputs are shown in Fig. 3. In this nominal scenario, reaction coefficients and boundary conditions in the sediment diagenesis model were chosen to correctly produce modern silica cycle fluxes as estimated by , and to approximately reproduce the fractional reverse weathering history in  to test its self-consistency (compare Fig. 2(C) to Fig. 3(I)). This model produces dissolved silica abundances through time ( Fig. 3(H)) consistent with literature estimates . In the Archean, reverse weather- Fig. 3. Carbon cycle and climate evolution over Earth history with self-consistent silica cycle evolution. The solid blue and dashed blue lines denote the negligible reverse weathering run as in Fig. 2. Grey, blue, and red shaded regions denote 95% confidence intervals for the carbon-silica cycle model with reverse weathering. Subfigure (G) shows biogenic silica (black), abiotic silica precipitation (blue) and reverse weathering (red) silica sink fluxes. Subfigure (H) shows the dissolved silica concentration in the ocean, and (I) shows the reverse weathering flux as a fraction of the total silica sink. Grey, blue, and red shaded regions denote 95% confidence intervals for carbon and silicon cycle variables. As in Fig. 2, the Proterozoic climate is warmer due to the elevated reverse weathering flux.
ing is lower than in the Proterozoic ( Fig. 3(E)) because reverse weathering fluxes are controlled by sedimentation rates, and these are assumed to be proportional to continental land fraction. Precambrian sedimentation rates and the precise functional relationship between continental land fraction and global sedimentation rates are unknown (Eriksson et al., 2013;Veizer and Mackenzie, 2003), but our assumed proportionality captures the expectation that siliciclastic deposition would be lessened with smaller subaerial continental area. Proterozoic temperature (Fig. 3(D)) and pCO 2 (Fig. 3(B)) are made higher by including reverse weathering. Earth's climate by ∼3-6 K in the Proterozoic. However, this result is sensitive to parameter choices and functional dependencies in the sediment diagenesis model, in addition to assumptions about initial (modern) conditions.
It is possible to create plausible model runs that show a different climate history. For example, Fig. 4 shows coupled C-Si cycle outputs where the modern Si cycle fluxes follow recent estimates described in . Here, the modern reverse weathering flux is considerably higher than early estimates by Tréguer and De La Rocha (2013) (4.5-4.9 Tmol Si/yr vs. 1.5 Fig. 4. Carbon cycle and climate evolution over Earth history with self-consistent silica cycle evolution and modern silica fluxes and parameters fitted to match . Grey, blue, and red shaded regions denote 95% confidence intervals for the carbon-silica cycle model with reverse weathering. Subfigure (G) shows biogenic silica (black), abiotic silica precipitation (blue) and reverse weathering (red) silica sink fluxes. Subfigure (H) shows the dissolved silica concentration in the ocean, and (I) shows the reverse weathering flux as a fraction of the total silica sink. Grey, blue, and red shaded regions denote 95% confidence intervals for carbon and silicon cycle variables. In this scenario, the modern reverse weathering flux is relatively large, and so there is only a modest change in reverse weathering between the Proterozoic and Phanerozoic. Consequently, surface temperature and CO 2 evolutions are comparable to the negligible reverse weathering case (solid blue and dashed blue lines).
Tmol/yr). The reaction coefficients were fitted to reproduce these modern conditions in Fig. 4 (see Appendix E for details). Similarly, Fig. 5 shows another example where reverse weathering has a minor effect on Earth's climate history. As in Fig. 3, we produce Tréguer and De La Rocha (2013) modern conditions, but in this case the fitted parameters ensured a negligible change in the reverse weathering flux between the Proterozoic and Phanerozoic.
It is also possible to select parameters from within the plausible ranges in Table 1 to deliberately engineer a large change in reverse weathering (Fig. 6). In this scenario, the Proterozoic is ∼15 K warmer than if reverse weathering were negligible. This would require some analyses of Proterozoic paleosols to have underestimated pCO 2 (Sheldon, 2006). The Proterozoic continental silicate weathering fluxes are 3-10 times higher than modern weathering fluxes in Fig. 6. This probably does not exceed the theoretical silicate weathering upper limit set by the supply of erodible rock, which is around 1000 Tmol/yr for modern Earth-like land fractions (Foley, 2015). However, if global weathering fluxes were as large as Fig. 6(E) suggests, then a larger fraction of the Earth's weatherable surface might become supply-limited, thereby decreasing the  Fig. 2. Grey, blue, and red shaded regions denote 95% confidence intervals for the carbon-silica cycle model with reverse weathering. Subfigure (G) shows biogenic silica (black), abiotic silica precipitation (blue) and reverse weathering (red) silica sink fluxes. Subfigure (H) shows the dissolved silica concentration in the ocean, and (I) shows the reverse weathering flux as a fraction of the total silica sink. Grey, blue, and red shaded regions denote 95% confidence intervals for carbon and silicon cycle variables. In this scenario, the reverse weathering flux is relatively insensitive to changes in the dissolved silica content of the ocean, and so there is only a modest change in reverse weathering between the Proterozoic and Phanerozoic. Consequently, surface temperature and CO 2 evolutions are nearly identical to the negligible reverse weathering case. sensitivity of global weathering to climate, and resulting in even higher temperatures and pCO 2 levels. For these reasons, we believe the scenarios presented in Fig. 3, 4, and 5 are more likely than the extreme Proterozoic warming in Fig. 6.
The sensitivity of these results to the pH-dependence of reverse weathering and Precambrian methane levels are discussed in Appendix A. In general, absolute changes in pH are small, and so the effect of a pH dependence proposed by Isson and Planavsky (2018) is minor (Fig. S1). This result is further validated by considering a more mechanistic diagenesis model that explicitly captures organic matter remineralization, dissolved inorganic carbon, carbonate alkalinity, and pH as a function of depth in anoxic sediments (Appendices A and C). Similarly, imposing high Precambrian methane levels results in somewhat higher temperatures and lower CO 2 levels, but does not dramatically change the climate evolution scenarios outlined above (Fig. S2).
Finally, it should be noted that the diminished effect of reverse weathering in the Archean is not guaranteed. If the model Fig. 6. Carbon cycle and climate evolution over Earth history with self-consistent silica cycle evolution and modern silica fluxes and parameters chosen to match  and to maximize the effect of reverse weathering in the Precambrian. The solid blue and dashed blue lines denote the negligible reverse weathering run as in Fig. 2. Grey, blue, and red shaded regions denote 95% confidence intervals for the carbon-silica cycle model with reverse weathering. Subfigure (G) shows biogenic silica (black), abiotic silica precipitation (blue) and reverse weathering (red) silica sink fluxes. Subfigure (H) shows the dissolved silica concentration in the ocean, and (I) shows the reverse weathering flux as a fraction of the total silica sink. Grey, blue, and red shaded regions denote 95% confidence intervals for carbon and silicon cycle variables. In this scenario, the reverse weathering flux is sensitive to changes in ocean dissolved silica content (k RW values are large and S A values are low). Consequently, Proterozoic pCO 2 and temperatures are very high.
is modified to remove any supply limit for and land fraction dependence of reverse weathering, then surface temperatures are warmer throughout the entire Precambrian, relative to the no reverse weathering case (Fig. 7). Whether Fig. 7 or Fig. 3-6 are better representations of Earth's climate evolution depends on which reverse weathering reactions dominated the global flux. Modern sediments with high silica concentrations are limited by terrigenous mineral supply Van Cappellen and Qiu, 1997), but it is unclear whether Precambrian sediments would be similarly supply limited.

Discussion
Why is there so much variability in our model outputs? Equation (9) for the sediment diagenesis model has reverse weathering linearly proportional to dissolved silica content. Ocean dissolved silica is estimated to have decreased by a factor of ∼20 during the transition from the Precambrian to Phanerozoic , and so we might naively expect the Precambrian reverse weathering flux to be much higher than the modern flux. However, global compilations reveal that the dissolved silica content in  6)). The solid blue and dashed blue lines denote the negligible reverse weathering run as in Fig. 2. Grey, blue, and red shaded regions denote 95% confidence intervals for the carbon-silica cycle model with reverse weathering. Subfigure (G) shows biogenic silica (black), abiotic silica precipitation (blue) and reverse weathering (red) silica sink fluxes. Subfigure (H) shows the dissolved silica concentration in the ocean, and (I) shows the reverse weathering flux as a fraction of the total silica sink. Grey, blue, and red shaded regions denote 95% confidence intervals for carbon and silicon cycle variables. Here, the absence of supply limit means that even in the Archean when the subaerial land fraction was smaller than today, the effect of reverse weathering on climate is potentially significant. modern ocean sediments-which is what controls reverse weathering rates-is typically a few hundred μM (Frings, 2017;. These elevated abundances relative to the ocean concentrations are due to the dissolution of biogenic silica . Additionally, Precambrian ocean dissolved silica concentration may have been below amorphous silica saturation . These factors combined mean that dissolved silica concentrations in Precambrian sediments were, at most, only a factor of a few higher than in modern sediments. Furthermore, sedimentation rates across most of the open ocean are low, and so reverse weathering reactions in the Precambrian often become supply limited across much of the seafloor (see supplementary figures), although even if this supply limit is omitted the reverse weathering response to changed [Si] is less than linear (Fig. 7). A linear relationship between reverse weathering and dissolved silica content likely only held under high sedimentation rates on the continental shelf.
Considerable uncertainty in the reaction coefficient, k R W , and the solubility of authigenic clays, S A , also contribute to variability in model results. If the reaction coefficient, k R W , falls at the lower end of the nominal range, then the change in reverse weathering that results from increasing ocean dissolved silica concentration is modest (Table 1). Similarly, high authigenic silicate saturation, S A , ensures a less than linear response in reverse weathering to an increase in dissolved ocean silica concentration (Table 1).
The model outputs presented in this study suggest enhanced reverse weathering could have warmed the Proterozoic by ∼5 K (Fig. 3). Alternatively, negligible Proterozoic warming ( Fig. 4 and 5) and dramatic warming of ∼15 K (Fig. 6) are also possible, although the latter scenario is probably the least likely because it contradicts Proterozoic CO 2 proxies (Sheldon, 2006) and predicts dramatically enhanced continental weathering fluxes. Model uncertainty is caused, in part, by the wide range of estimates of the modern reverse weathering flux. If the modern flux is large, as suggested by  and illustrated in our Fig. 4, then there is less capacity for a dramatic increase in reverse weathering fraction at the Phanerozoic-Proterozoic boundary. However, when considering different estimates for the modern reverse weathering flux, the broader implications for the modern carbon and alkalinity budgets must be considered. One potential problem with the large reverse weathering fluxes proposed in  is that for both the modern silica and carbon cycles to be balanced, the modern carbon outgassing flux must be 3-7 Tmol C/yr (Fig. 4). If a larger portion of the alkalinity sourced from the weathering of silicates and carbonates is consumed by reverse weathering, then the resulting lowered carbonate burial results implies lower outgassing to balance the carbon budget. This required outgassing flux is low compared to most estimates, which converge around 6-10 Tmol C/yr (Berner, 2004;Lee and Lackey, 2015). Another important caveat here is that our modeling framework implicitly assumes that all forms of alkalinity are interchangeable. In practice, an elevated modern reverse weathering sinks for K + cations may not affect the carbon budget because K + does not contribute to carbonate sinks, and therefore elevated K-rich clays may only require a larger riverine supply from continental weathering. Better empirical constraints are needed not only on the total modern reverse weathering flux, but also on the individual cation contributions, as well as the possible effect of cation-exchange reactions in both sedimentary and hydrothermal systems.
Kanzaki and Murakami (2018a) developed a detailed silicate weathering model to calculate how the temperature dependence of silicate weathering changes as a function of atmospheric composition. Specifically, higher pCO 2 results in a lower effective activation energy for silicate weathering and therefore a weaker weathering feedback (Kanzaki and Murakami, 2018b). We do not consider this feedback in this study to isolate the effect of reverse weathering. For plausible Precambrian pCO 2 ranges (10 −3 to 1 atm), the weathering model presented in Kanzaki and Murakami (2018b) predicts effective activation energies ranging from ∼140 kJ/mol to ∼14 kJ/mol, which corresponds to e-folding temperatures ranging from 5 K to 45 K. The combined uncertainties in the CO 2 -dependence and temperature dependence of continental weathering in our model  effectively encompass this range of e-folding temperatures.
In contrast to evidence for snowball episodes in the Neoproterozoic and early Paleoproterozoic, the apparent absence of glaciation between 1.0-1.8 Ga has typically been interpreted as implying a persistently warm mid-Proterozoic climate (Kasting, 2005). This is difficult to reconcile with lower solar luminosity and Proterozoic pCO 2 estimates from paleosols (Sheldon, 2006), and so hypotheses for warming the Proterozoic have been proposed such as elevated atmospheric methane (Catling et al., 2002;Pavlov et al., 2003) or nitrous oxide (Buick, 2007;Roberson et al., 2011;Stanton et al., 2018). The extent to which marine sulfate and/or inefficient degradation of organic matter may have prevented Proterozoic atmospheric CH 4 accumulation Laakso and Schrag, 2019;, and the possible contribution from terrestrial CH 4 fluxes (Zhao et al., 2017) are subjects of ongoing debate.
This study, and that of , suggest that elevated reverse weathering could easily provide ∼5 K warming in the Proterozoic, thereby helping explain mid-Proterozoic warmth. The decline in ocean dissolved silica concentrations in the Neoproterozoic and Paleozoic caused by the rise of ecologically significant biogenic silica precipitation results in cooling (e.g. Fig. 3, 6, S2) approximately coincident with Neoproterozoic snowball episodes, although note that the emergence of ecologically significant biogenic silica precipitation may postdate snowball episodes (see section 1). Our nominal model outputs predict a diminished role for reverse weathering-and therefore no additional warming-in the Archean due to lower continental land fraction (Korenaga et al., 2017) and increasing total carbon throughput (Marty et al., 2019). Alternatively, if reverse weathering is not limited by the supply of silicates from the continents, then reverse weathering may also contribute to warmer temperatures in the Archean (Fig. 7). Finally, it should be noted that elevated Proterozoic temperatures are not necessarily required by our model (e.g. Fig. 4, 5), and may not even be necessary given proposed evidence for mid-Proterozoic glaciation (Geboy et al., 2013;Kuipers et al., 2013), although some purported mid Proterozoic glaciations may in fact be Neoproterozoic in age (Alvarenga et al., 2019).
Our results are highly sensitive to unconstrained parameters in sedimentary diagenesis models, and so additional empirical constraints are required to determine the importance of Precambrian reverse weathering for Earth's carbon cycle evolution and possible contribution to Proterozoic warmth. Possible constraints could come from Si isotopes, better quantification of clay minerals in the Precambrian, or an improved mechanistic understanding of diagenetic processes, which we discuss subsequently, in turn.  performed an inverse analysis of Si isotope data over Earth history to predict the silica clay sink as a fraction of total silica deposition over Earth history. This cannot be used to constrain reverse weathering through time because this total clay deposition term does not differentiate between authigenic clay formation that consumes alkalinity (i.e. reverse weathering), and clay formation that produces alkalinity (kaolinite-forming weathering reactions). Since the Si isotopic abundance of the authigenic clay reservoir through time is unknown, the isotope system cannot be solved to isolate the reverse weathering fraction . Instead, Trower and Fischer (2019) adopted the f r w curve in  to predict the Si isotopic ratio of the authigenic clay reverse weathering reservoir through time. The f r w model outputs presented in the various scenarios in this paper can similarly be used to predict the Si isotopic ratio of reverse-weathering clays through time.
Given current available data, the Si isotope system cannot be used to constrain the reverse weathering fraction through time. Fig. S9 in Appendix D shows a statistical analysis of the Si isotope record presented in  where the reverse weathering fraction is solved for rather than imposed. This analysis demonstrates that more Si isotope data, preferably data that separates alkalinity-consuming clays and alkalinity-producing clays, are required to constrain the reverse weathering fraction.
Isson and Planavsky (2018) presented a literature compilation of greenalite, minnesotaite, and stilpnomelane through time as suggestive evidence for elevated clay formation in the Precambrian. These minerals were chosen because of their probable marine origin. But since all three minerals are iron-bearing and/or commonly found in banded iron formations, their occurrence through time may reflect banded iron formation deposition in the Precambrian. In particular, the minnesotaite occurrence data closely resembles the record of banded iron formation deposition. The Si isotope analysis presented in  demonstrated that banded iron formations were probably a negligible contributor to the Precambrian silica budget. Additionally, the Proterozoic continental shelf waters were likely predominantly euxinic or oxic, which would predict K, Na, or Mg-bearing reverse weathering products rather than Fe-bearing clays such as greenalite and minnesotaite. Note also that the literature compilation in  is based on "number of papers published" rather than distinct greenalite/minnesotaite occurrences, and thus does not distinguish early diagenesis from later alteration.
Progress could also be made with an improved mechanistic understanding of authigenic clay kinetics within anoxic sediments. Better experimental constraints on clay formation kinetics, including their pH and redox sensitivities, the saturation states of authigenic phases, and the kinetics of cation cycling in different marine settings could enable precise reverse weathering model predictions. Although our analysis suggests the pH dependence of reverse weathering is of secondary importance compared to changes in ocean dissolved silica content (Appendix A), an enhanced pHdependence cannot be completely ruled out without a more detailed mechanistic understanding of clay forming feedbacks. Such an approach could also help shrink the uncertainties in Precambrian climate calculations (Fig. 3, 4, 5, 6, and 7).
A final, related caveat is that our reverse weathering formulation effectively assumes that silica abundances are the primary limitation on global reverse weathering fluxes. We accounted for possible cation limitations empirically by sampling a range of efolding depths, authigenic saturation concentrations, and reaction coefficients in our sediment diagenesis model, and by exploring a rudimentary mechanistic model in the supplementary materials with an explicit treatment of pH and organic matter oxidation. However, the empirical ranges we adopt are derived primarily from modern settings, including river deltas with sizeable terrigenous inputs of weathered cations and organic matter (e.g. Michalopoulos and Aller, 2004). Precambrian sediments may have been much more severely cation-limited, and so it is possible that an improved mechanistic treatment of these other limitations to authigenic reactions might disfavor strongly elevated Precambrian reverse weathering scenarios.

Conclusions
The geological carbon cycle controls Earth's climate and ocean chemistry evolution. The carbon cycle may be modulated by the silica cycle because reverse weathering reactions remove ocean alkalinity whilst retaining atmospheric CO 2 . A decline in reverse weathering caused by the proliferation of biogenic silica precipitation at the end of the Proterozoic has been hypothesized to have triggered dramatic climate cooling and could help explain a warm Proterozoic climate.
By applying a coupled C-Si cycle model, we find that a broad range of climate evolutions is possible. Under most plausible scenarios, elevated reverse weathering in the Precambrian results in a Proterozoic climate that is ∼3-6 K warmer than if reverse weathering were ignored. However, self-consistent scenarios with negligible changes in reverse weathering throughout Earth historyimplying that the silica cycle does not affect climate evolution-are also consistent with modern constraints. Scenarios with extremely large changes in reverse weathering, implying ∼15 K warming in the Proterozoic, are less likely but cannot be ruled out.
The evolution of biogenic silica precipitation and subsequent decline in dissolved ocean silica content has the strongest reverse weathering effect on Earth's carbon cycle evolution. In contrast, a strong pH-dependence for reverse weathering reactions has little effect on Earth's climate or ocean pH evolution.
Key uncertainties are the modern reverse weathering flux, and parameters in a sediment diagenesis model that affect the global reverse weathering rate. In particular, the rate coefficient for reverse weathering reactions has large uncertainty, and if it is small, reverse weathering could be relatively unimportant in the Precambrian. Also, whether authigenic clays are highly soluble or not can affect the efficacy of reverse weathering.
Additional empirical constraints such as those from Si or Li isotopes, better quantification of reverse weathering reaction kinetics, and refined estimates of Phanerozoic reverse weathering fluxes are required to fully assess the magnitude of the effect of reverse weathering on Earth's climate evolution. An improved mechanistic treatment of reverse weathering kinetics including cation-and terrigenous-matter limitations on authigenic reactions would be especially helpful for evaluating the likelihood of large Precambrian RW fluxes.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Supplementary material Appendix A: Additional sensitivity tests
The effect of pH Fig. S1 is identical to Fig. 3 except that the equation governing the reverse weathering flux has an ocean pH to the power of 22 dependence when the biogenic silica precipitation flux drops to zero, following an assumed equation in . This has a relatively minor effect on carbon cycle evolution: the Precambrian reverse weathering flux is lower than in the no pH-dependence case (Fig.   S1E), and consequently Precambrian temperatures and pCO2 levels are slightly lower than the no pHdependence case ( Fig. S1B and D). This effect is small because cooling introduced by the pH-dependent reverse weathering decline is offset by the silicate weathering thermostat (the carbon cycle must still be balanced). Consequently, fractional reverse weathering (Fig. S1I) is lower than in the no pH-dependence case and surface temperature and pCO2 are only slightly lower than the no pH-dependence case. It is also worth noting that even with a pH 22 dependence, absolute changes in pH are small, and so the effect of this pH dependence is relatively minor. For example, a change from pH=7.4 to pH=7.2 across the Proterozoic results in a 1.8-fold change in reverse weathering. This is much smaller than the 8-fold change in reverse weathering between the Phanerozoic and Proterozoic (Fig. S1E).
One potential criticism of this approach to incorporating a pH-dependence is that an overall ocean pH 22 dependence doesn't account for possible feedbacks in sediments between organic carbon oxidation, alkalinity consumption, and pH. To explore these effects, we modified our sediment diagenesis model to track organic carbon, alkalinity, dissolved inorganic carbon, and pH (see Appendix C for details). The overall sensitivity to or reverse weathering to pH was less sensitive than pH 22 , suggesting that the conclusions drawn from Fig. S1 are reasonable. However, it is possible that a more complex model that 2 better captures the pH and cation-concentration dependence of different amorphous clay saturation states may possess a stronger overall pH dependence.

High Precambrian methane
Fig. S2 is identical to Fig. 3 except that we have imposed CH4 abundances of 100 ppm during the Proterozoic and 1% in the Archean Kasting & Brown 1998;) as representative of a 'high CH4' scenario. Unsurprisingly, this results in a somewhat warmer Precambrian than in Fig. 3, given that imposed CH4 has no negative feedback on CH4 itself.
Carbon dioxide levels are lower than in Fig. 3 (and pH is higher) to ensure a balanced carbon cycle against methane induced warming.   Fig. 2. Grey, blue, and red shaded regions denote 95% confidence intervals for the carbon-silica cycle model with reverse weathering (identical to Fig. 3 except for addition of high methane). Subfigure (G) shows biogenic silica (gray), abiotic silica precipitation (blue) and reverse weathering (red) silica sink fluxes. Subfigure (H) shows the dissolved silica concentration in the ocean, and (I) shows the reverse weathering flux as a fraction of the total silica sink. Grey, blue, and red shaded regions denote 95% confidence intervals for carbon and silicon cycle variables.

Appendix B: Carbon-silica cycle model
The full set of equations governing the time-evolution of carbon, alkalinity, and silica reservoirs are as follows: Here, The power law exponent is taken from Carroll et al. (1998). A sediment diagenesis model (see section 2.2) cannot easily be used to predict abiotic silica precipitation because abiotic silica precipitation probably occurs sporadically in evaporative settings (Maliva et al. 2005). Our parameterization of abiotic silica precipitation is analogous to standard parameterizations of carbonate precipitation using a power law related to supersaturation by an empirical exponent (Opdyke & Wilkinson 1993). The reaction coefficient , As noted in the main text, our results are relatively insensitive to choices of Abio k .

Appendix C: Sediment diagenesis model
This section outlines the sediment diagenesis model used to create the parameterizations for biotic silica precipitation and reverse weathering fluxes in the carbon cycle model. The precipitation of authigenic clays and dissolution of biogenic silica in ocean sediments is modeled using a reactiondiffusion model similar to those described in Schink et al. (1975), Aller (2014), and .
The governing equations are: Here, Here, B K =10 -9 s -1 is the biogenic silica dissolution coefficient,  (1/cm) governs the decay of silica dissolution with depth due to opal maturity, and f C =1.0 μmol/cm 3 is the saturation concentration for dissolved silica (DeMaster 2003;Rickert et al. 2002). Note that the dissolution of biogenic silica ceases when dissolved concentrations exceed the saturation concentration. Similarly, the reverse weathering reactions cease when dissolved silica is undersaturated with respect to authigenic phases, or when >50% of sediments (by mass) have already undergone clay forming reactions, except where stated otherwise. The remaining variables in equation (8) are defined in the main text.
The system of equations described above are solved using the following boundary conditions:      open ocean curves were calculated using the sediment diagenesis model described above.
The polynomial fit to the total burial (in Tmol Si/yr) in Fig S5 is given by the following equation: Here, t is time in Ga. Biogenic silica precipitation is assumed to gradually fall to zero by around 1 Ga There is a gradual transition between the polynomial fit with biogenic precipitation, 1 R f − , (Fig. S6) and the polynomial fit before biogenic precipitation, 2 R f − (Fig. S7). The polynomial fits to the sediment diagenesis model using different parameters (i.e. those used to generate Fig. 4, Fig. 5, Fig. 6, and Fig. 7) are provided in the published code, which is available on the website of the lead author. 16

Mechanistic model modifications
This section describes the modifications made to the sediment diagenesis model to explore the effect of pH-dependent reverse weathering reactions. Ocean alkalinity, dissolved inorganic carbon, and organic carbon are explicitly tracked in the model: [ALK] [DIC] 1 , , , ,  Figure S8 shows model outputs from these calculations. We see that this change in pH results in an approximately 3-fold change in the reverse weathering flux. This is considerably less that what we would expect from a pH 22 dependence, (7.4/6.9) 22.4 ~ 5-fold, therefore implying that the overall ocean pH dependence adopted in Fig. S1 if anything, overestimates the effect of pH on reverse weathering. 19 Appendix D: Silicon isotope calculations Fig. S9 was calculated using the isotopic data for chert and modern authigenic clays in . Mass balance between bulk silicate Earth inputs and chert and clay outputs was assumed, the iron formation component was assumed to be negligible, and the fractionation between chert and authigenic clays formed by reverse weathering is assumed to be -2‰ ).
This figure shows fractional reverse weathering, rw f , over Earth history. Clearly, the reverse weathering fraction is not well constrained by the existing isotope record. This is because the isotope system is degenerate. Fig. S9 (21) for the reverse weathering sink, discarding all rw f values greater than unity or less than zero. The degeneracy in this isotope system produces the large spread in Fig. S9. Progress requires separate constraints on 30 Si  in clays from reverse weathering (alkalinity consuming reactions) and kaolinite clays (alkalinity producing reactions). 21

Appendix E: Method for fitting unknown parameters
There are six unknown parameters in the sediment diagenesis model: 1) RW k (μmol/cm3/s), the reaction coefficient for reverse weathering reactions,

2)
A z (cm), the e-folding depth governing the decrease in reverse weathering with depth

4)
B K (s -1 ), the reaction coefficient for biogenic silica dissolution 5)  (1/cm), the coefficient governing the decrease in biogenic silica dissolution with depth 6) B F (μmol/cm 2 /s), the diffusion flux of biogenic silica at the sediment-ocean interface.
Parameters (1), (2), and (3) are discussed in the main text and shown in Table 1. For parameters (4), (5), and (6), plausible ranges and chosen values are shown in Table S1. Parameters (4), (5), and (6) are less important for controlling the Precambrian reverse weathering flux since they relate to the rate of biogenic silica burial. However, their values are important for fitting the modern reverse weathering flux as the dissolved silica available for reverse weathering is mostly sourced from biogenic silica dissolution (as opposed to in the Precambrian where it is exclusively sourced from seawater). As noted in the main text, parameter values were either chosen heuristically from plausible ranges (unbolded numbers in Table 1 and Table S1) or fitted using Markov Chain Monte Carlo (MCMC) inverse methods (bolded numbers in Table 1 and Table S1). Inverse method for fitting parameters MCMC methods were used to fit model parameters using approximate observational constraints. This is not intended to be a rigorous solution to the inverse problem, but rather a 'ballpark approach' to ensure the parameters used calculate carbon-silica cycle behavior in the Precambrian can at least approximately reproduce modern conditions.

Appendix F: Comparison to carbon cycle model of Isson and Planavsky (2018)
The carbon cycle model used by  is intrinsically more sensitive to changes in reverse weathering than the carbon cycle model used in this study. This appendix presents comparisons between steady state calculations from both models, as well as analytic flux-balanced calculations.
In these calculations, we adapted our carbon cycle model to match that of  as best as possible. For example, we removed seafloor weathering and adopted the following overall function for silicate weathering with the same weathering feedback exponent: Solar luminosity, continental land fraction, and the biological enhancement of weathering were held constant for these steady state calculations. The sediment diagenesis model is not used. Instead, an increase in fractional reverse weathering is imposed and the steady state response of the carbon cycle parameters is calculated. Parameter values were chosen to match those in  wherever possible, although parameter choices were sometimes unclear since the source code is unavailable.
In addition to comparing the two models, we also present analytic calculations of the relationship between reverse weathering flux and atmospheric CO2. In steady state, carbon fluxes and alkalinity fluxes are balanced: By combining these two equations and substituting equation (22) we obtain the following expression:   Fig. S10 and S11). In both Fig. S10 and S11 analytic calculations of pCO2 as a function of reverse weathering match our numerical model exactly. The quantitative differences between  and our model are not explained. represent different ALK/Si ratios for reverse weathering (4, 3.5, 3, 2.5, 2). The qualitative behavior of the two models is similar, but Isson and Planavsky (2018) is more sensitive to changes in reverse weathering than our model for unknown reasons. reverse weathering sink, i.e., the fraction of total Si consumed through reverse weathering. The five contours in each subfigure represent different ALK/Si ratios for reverse weathering (4, 3.5, 3, 2.5, 2).