A warm or a cold early Earth? New insights from a 3-D climate-carbon model

TL;DR

This study uses a 3D climate-carbon model to explore early Earth's climate, suggesting it was likely cold or temperate with moderate CO2 levels, challenging previous hot Earth hypotheses based on isotope data.

Contribution

The paper introduces a 3D climate-carbon model to reassess early Earth's climate, incorporating cloud feedback and weathering effects, providing new insights into its temperature and CO2 stability.

Findings

Early Earth likely had cold or temperate climates.

Moderate CO2 levels (0.1-0.36 bar) could sustain these climates.

Weathering processes stabilized CO2 levels, preventing extreme greenhouse conditions.

Abstract

Oxygen isotopes in marine cherts have been used to infer hot oceans during the Archean with temperatures between 60{\deg}C (333 K) and 80{\deg}C (353 K). Such climates are challenging for the early Earth warmed by the faint young Sun. The interpretation of the data has therefore been controversial. 1D climate modeling inferred that such hot climates would require very high levels of CO2 (2-6 bars). Previous carbon cycle modeling concluded that such stable hot climates were impossible and that the carbon cycle should lead to cold climates during the Hadean and the Archean. Here, we revisit the climate and carbon cycle of the early Earth at 3.8 Ga using a 3D climate-carbon model. We find that CO2 partial pressures of around 1 bar could have produced hot climates given a low land fraction and cloud feedback effects. However, such high CO2 partial pressures should not have been stable because of the weathering of terrestrial and oceanic basalts, producing an efficient stabilizing feedback. Moreover, the weathering of impact ejecta during the Late Heavy Bombardment (LHB) would have strongly reduced the CO2 partial pressure leading to cold climates and potentially snowball Earth events after large impacts. Our results therefore favor cold or temperate climates with global mean temperatures between around 8{\deg}C (281 K) and 30{\deg}C (303 K) and with 0.1-0.36 bar of CO2 for the late Hadean and early Archean. Finally, our model suggests that the carbon cycle was efficient for preserving clement conditions on the early Earth without necessarily requiring any other greenhouse gas or warming process.

Tables

Table 1. Table 1: Global mean values from the GCM. Columns from left to rights are: pCO 2 (in bar), pCH 4 (in mbar), the mean surface temperature (in ∘ C), the mean planetary albedo, the mean precipitation rate (in mm/day), the shortwave cloud forcing (in W/m 2 ) and the longwave cloud forcing (in W/m 2 ).

pCO2 (bar)	pCH4 (mbar)	Tsurf (∘C)	Planetary albedo	Precipitation (mm/day)	SW forcing (W/m2)	LW forcing (W/m2)
0.01	0	-11.8	0.4	1.8	-38.3	+24.5
0.01	2	2.7	0.33	2.2	-42.0	+20.2
0.1	0	7.4	0.34	2.8	-47.1	+24.9
0.1	2	21.3	0.28	3.5	-39.8	+18.1
0.3	0	22.1	0.31	3.6	-42.4	+20.8
0.3	2	34.7	0.26	4.4	-32.6	+13.8
0.5	0	34.7	0.27	4.6	-30.4	+16.0
0.5	2	49.9	0.23	5.1	-22.8	+12.3
1.0	0	59.5	0.24	5.7	-19.4	+14.2
1.0	2	67.5	0.23	6.1	-20.3	+14.0

Table 2. Table 2: Equilibrium states from the carbon model coupled to the GCM. Columns from left to right are: the land fraction, the oceanic crust thickness (in km), the spreading rate, the recycling rate, the presence (2 mbar) or absence of methane, the outgassing rate, pCO 2 (in mbar), the mean surface temperature (in C ∘ superscript 𝐶 {}^{\circ}C ), the mean oceanic pH, the mean dissolved inorganic carbon (in mol/m 3 ), the CO 2 flux outgassed from mid-oceanic ridges and from arc volcanoes, and the CO 2 sink by continental weathering and seafloor weathering (all fluxes in 1e12 mol/yrs). The first case, with land fraction equal to 1, is the present-day Earth. The last case, with SR=0, is the heat-pipe scenario with no plate tectonics.

Boundary conditions						Results
Land fraction	Oceanic crust (km)	SR	Recyling	CH4	rdeg	pCO2 (mbar)	Tsurf (∘C)	pH	DIC	Fmor	Farc	Fcon	Fsfw
1	6	1	0.1	no	1	0.285	14.9	8.4	2.4	1.6	7.7	-10.5	-1.1
0.0	6	11.3	0.1	no	2.5	118	8.5	6.70	32.2	18	11.2	0.0	-29.2
0.0	6	11.3	0.3	no	2.5	225	16.5	6.60	39.8	18	19.5	0.0	-37.5
0.0	18	3.8	0.1	no	2.5	304	22.3	6.56	48.9	18	11.2	0.0	-29.2
0.0	18	3.8	0.1	yes	2.5	190	27.0	6.66	35.4	18	11.2	0.0	-29.2
0.0	18	3.8	0.3	no	2.5	359	25.8	6.53	52.4	18	19.5	0.0	-37.5
0.0	18	3.8	0.3	yes	2.5	243	30.6	6.61	39.6	18	19.5	0.0	-37.5
0.01	6	11.3	0.1	no	2.5	118	8.5	6.70	32.21	18	11.2	-0.05	-29.2
0.01	6	11.3	0.3	no	2.5	231	16.9	6.65	45.10	18
0.01	18	3.8	0.1	no	2.5	303	22.3	6.56	49.51	18
0.01	18	3.8	0.1	yes	2.5	186	26.8	6.67	35.39	18
0.01	18	3.8	0.3	no	2.5	358	25.8	6.53	52.73	18
0.01	18	3.8	0.3	yes	2.5	240	30.5	6.62	39.64	18
0.12	6	11.3	0.1	no	2.5	117	8.5	6.70	32.06	18	11.6	-0.4	-29.2
0.12	6	11.3	0.1	yes	2.5	97	20.2	6.80	25.52	18
0.12	6	11.3	0.3	no	2.5	235	17.3	6.66	46.69	18	23.7	-3.1	-38.7
0.12	6	11.3	0.3	yes	2.5	135	23.3	6.75	31.47	18
0.12	18	3.8	0.1	no	2.5	292	21.4	6.59	50.13	18	12.6	-4.1	-26.5
0.12	18	3.8	0.1	yes	2.5	163	25.2	6.71	34.38	18
0.12	18	3.8	0.3	no	2.5	356	25.6	6.57	56.31	18
0.12	18	3.8	0.3	yes	2.5	219	29.0	6.69	41.89	18
0.26	6	11.3	0.1	no	2.5	121	8.8	6.70	32.50	18	12.1	-0.64	-29.5
0.26	6	11.3	0.3	no	2.5	245	18.0	6.66	48.47	18
0.12	6	0	0	no	10	257	18.8	6.60	45.75	31.9	1.1	-0.9	-32.1

Equations

$$f(\ge m)=f_0\times {\left(\frac{m}{m_{max}}\right)}^{-b}$$

$$F_{mor}^0+F_{arc}^0+F_{sfw}^0=0$$

$$F_{mor}+F_{arc}+F_{sfw}+F_{ej}=0$$

$$pC{O}_2=pC{O}_2^0\times {\left(\frac{S{R}_0}{SR}\right)}^{1/\alpha }{\left(\frac{1-R_0}{1-R}\right)}^{1/\alpha }{\left(1-\frac{(1-R)\times F_{ej}}{F_{mor}+F_{deep}}\right)}^{1/\alpha }$$

$$F_{\mathrm{mor}}=F_{\mathrm{mo}{\mathrm{r}}_0}\times \frac{SR}{S{R}_0}\times \frac{H}{H_0}$$

$$F_{\mathrm{arc}}=F_{\mathrm{sub}}\times R+F_{deep}$$

$$F_{\mathrm{con}}=f_{\mathrm{Ca}}^{bas}+f_{\mathrm{Ca}}^{gra}+f_{\mathrm{Mg}}^{bas}+f_{\mathrm{Mg}}^{gra}$$

$$f_{\mathrm{Ca}}^{bas}=20.94\times 10^{12}\times e^{fbas}\times runoff\times \frac{sr{f}_{\mathrm{bas}}}{sr{f}_{\mathrm{cont}}}\times \frac{sr{f}_{\mathrm{cont}}}{sr{f}_{cont0}}$$

$$f_{\mathrm{Ca}}^{gra}=3.49\times 10^{12}\times e^{fgra}\times runoff\times \left(1-\frac{sr{f}_{\mathrm{bas}}}{sr{f}_{\mathrm{cont}}}\right)\times \frac{sr{f}_{\mathrm{cont}}}{sr{f}_{cont0}}$$

$$f_{\mathrm{Mg}}^{bas}=19.86\times 10^{12}\times e^{fbas}\times runoff\times \frac{sr{f}_{\mathrm{bas}}}{sr{f}_{\mathrm{cont}}}\times \frac{sr{f}_{\mathrm{cont}}}{sr{f}_{cont0}}$$

$$f_{\mathrm{Mg}}^{gra}=3.31\times 10^{12}\times e^{fgra}\times runoff\times \left(1-\frac{sr{f}_{\mathrm{bas}}}{sr{f}_{\mathrm{cont}}}\right)\times \frac{sr{f}_{\mathrm{cont}}}{sr{f}_{cont0}}$$

$$sr{f}_{\mathrm{bas}}=\frac{F_{\mathrm{arc}}}{F_{arc0}}\times \frac{sr{f}_{\mathrm{cont}}}{sr{f}_{cont0}}\times sr{f}_{bas0}$$

$$f_{sfw}=\frac{SR}{S{R}_0}\times \sum _i{d}_{i}exp\left(\frac{-E_i}{R{T}_p}\right)C_i^n$$

$$f_{cd}=k_{cd}(\Omega -1{)}^{1.7}$$

$$f_{atm-ocean}=\alpha (pC{O}_2^{atm}-pC{O}_2^{ocean})$$

$$SR=S{R}_{\ast }\times {\left(\frac{T_p-T_0}{T_p^{\ast }-T_0}\frac{\eta (T_p^{\ast })}{\eta (T_p)}\right)}^{2/3}$$

$$m_{ej}=\frac{\pi D_{tc}^3{\rho }_c}{16\sqrt{2}}$$

$$D_{tc}=1.161{\left(\frac{{\rho }_i}{{\rho }_c}\right)}^{1/3}{L}^{0.78}{v}_i^{0.44}{g}_E^{-0.22}si{n}^{1/3}\theta$$

$$v_i=v_{0}exp\left(-\frac{3(C_{Datm}{\rho }_{air}{H}_{atm}+C_{Docean}{\rho }_{water}{H}_{ocean})}{4{\rho }_{i}Lsin\theta }\right)$$

$$V_m=8.9\times 10^{-12}Esin\theta$$

Full text

1 Introduction

Life likely emerged on Earth before 3.5 Ga, during the Hadean or the early Archean [Buick et al., 1981; Nisbet and Sleep, 2001; Bell et al., 2015]. Studying the atmosphere and the climates of the Earth during its first billion year is critical for understanding the environment in which life emerged and developed. Oxygen isotopic ratio of Archean marine cherts have been interpreted to suggest hot oceans with temperatures between 60∘C (333 K) and 80∘C (353 K) [Knauth and Lowe, 2003; Robert and Chaussidon, 2006], compatible with the thermophiles inferred from evolutionary models of the ancient life [Gaucher et al., 2008] and possibly with the indications for a low ocean water viscosity [Fralick and Carter, 2011]. However, such an interpretation has been strongly debated [Kasting and Howard, 2006; van den Boorn et al., 2007; Marin-Carbonne et al., 2014; Tartese et al., 2016] and other analyses suggest temperate oceans with temperatures lower than 40∘C [Hren et al., 2009; Blake et al., 2010]. Moreover, at 3.5 Ga, the presence of glacial rocks at 20-40∘ latitude implies temperatures below 20∘C [de Wit and Furnes, 2016], at least episodically.

In addition, previous modeling of the carbon cycle on the early Earth by Sleep and Zahnle [2001] and Zahnle and Sleep [2002] suggested that the Archean was cold unless another strong greenhouse gas was present. They also suggested that the Hadean was likely very cold because of the weathering of impact ejecta, particularly during the Late Heavy Bombardment (LHB). Therefore, the temperature of the early oceans remains an open question.

Here, we use a 3D climate-carbon model with pCO2 ranging from 0.01 to 1 bar in order to determine physically plausible climates of the late Hadean and early Archean Earth. In section 2, we describe the climate simulations. In section 3 and 4, we analyze the carbon cycle and the effect of the LHB with respect to carbon cycle responses. We finish with a conclusion in section 5.

2 Climate modeling

2.1 Description of the model

We simulated the atmosphere of the early Earth using the Generic LMDZ GCM (Global Climate Model). This model solves the primitive hydrostatic equations of meteorology using a finite difference dynamical core on an Arakawa C grid. It uses robust and general parameterizations in order to simulate planets very different from the present-day Earth. The Generic LMDZ GCM has already been used for studying the temperate and cold climates of the Archean Earth [Charnay et al., 2013], the climates of early Mars with high CO2 pressures [Forget et al., 2013], and the runaway greenhouse effect on terrestrial planets with high amounts of water vapor [Leconte et al., 2013]. It is thus adapted for simulating the early Earth with high amounts of CO2 and potentially hot and moist climates.

The radiative scheme is based on the correlated-k method. At a given pressure and temperature, correlated-k coefficients in the GCM are interpolated from a matrix of coefficients stored in a 7 $\times$ 9 temperature and log-pressure grid: $T$=100, 150, 200, 250, 300, 350, 400 K, $p$ = 10*-1*, 100, 101, …, 107 Pa. We used 36 spectral bands in the thermal infrared and 38 at solar wavelengths. Sixteen points were used for the g-space integration, where g is the cumulative distribution function of the absorption data for each band.

Simulations were performed with a horizontal resolution of 64$\times$48 (corresponding to resolutions of 3.75∘ latitude by 5.625∘ longitude) and with 25 vertical layers with the lowest midlayer level at 5 m and the top level at 0.5 hPa. For the cloud microphysics, we either fixed the radii of water cloud particles (e.g. 12 $\mu$m for liquid droplets and 35 $\mu$m for icy particles for present-day Earth) or we fixed the density of cloud condensation nuclei (CCN) by mass of air (e.g. 5$\times$106 particles/kg for liquid droplets and 2$\times$104particles/kg for icy clouds for present-day Earth). As in Charnay et al. [2013], the oceanic transport and the sea ice formation were computed with the 2-layer oceanic model from Codron [2012].

To investigate the early Earth at 3.8 Ga, we used a solar constant of 1024 W/m2 corresponding to 75$\%$ of the present value (1361 W/m2). We ran simulations with no land. This hypothesis is valid for small continental surface fractions (e.g. $<$30$\%$ of the present-day fraction or $<$10$\%$ of the total Earth’s surface) as we assumed at this time [Flament et al., 2008] (see also section 3.1), although the land fraction is a matter of dispute [Viehmann et al., 2014]. Such small continental fraction should indeed have a negligible direct impact on the global climate. We also assumed that Earth’s rotation period was 14h. We used an atmospheric composition with 1 bar of nitrogen, a partial pressure of CO2 (pCO2) from 0.01 to 1 bar, and either no CH4 or 2 mbar of CH4, spanning the range of plausible methane concentrations estimated by a model for early ecosystems [Kharecha et al., 2005]. Wordsworth et al. [2017] showed that CO2-CH4 collision-induced absorption could have produced a strong warming on early Mars if the CH4 mixing ratio was higher than 1$\%$. With our CH4 mixing ratio of 0.1-0.2 $\%$, CO2-CH4 CIA produce a warming lower than 0.1 K with the 1D version of the model. We thus neglected CO2-CH4 CIA in our 3D simulations.

2.2 Results

Figure 1 shows the global mean and zonal mean surface temperatures obtained with the GCM for present-day properties of clouds (Global mean values are gathered in table 1). Using 3.8 Ga conditions, the Earth falls into a full glaciation for the present-day pCO2 (see Charnay et al. [2013]), but it is never fully ice-covered with at least 0.01 bar of CO2. For this case with no methane, the mean surface temperature is -11.8∘C, well below the freezing point, but there is still a cold ice-free equatorial belt between -37∘N and +37∘N. With 0.1 bar of CO2, the climate is temperate (global mean surface temperature around 7.4∘C without CH4 and 21.3∘C with CH4). Ice-free conditions are obtained for pCO2 higher than around 0.3 bar without methane or slightly higher than 0.1 bar with methane, corresponding to a mean surface temperature higher than around 23∘C. Our model has a climate sensitivity with pCO2 very similar to the GCM from Wolf and Toon [2013] going up to 0.3 bar. With 0.5 bar of CO2, the climate becomes warm (mean surface temperature around 34.7∘C without CH4 and 49.9∘C with CH4). With 1 bar of CO2, the mean surface temperature is around 59.5∘C without CH4 and around 67.5∘C with CH4.

Ozak et al. [2016] showed that CO2 collisional line mixing reduces the atmospheric opacity between 400 and 550 cm*-1* producing a significant cooling ($\sim$10∘C) for a 1 bar CO2 atmosphere for early Mars. Since we do not include collisional line mixing, our model might overestimate temperatures for high pCO2. However, for the early Earth with a 1 bar CO2 atmosphere, the peak of thermal emission from the surface would be around 1110 cm*-1* (compared to around 730 cm*-1* for early Mars), so quite far from the 400-550 cm*-1* spectral region. We therefore expect a weaker cooling effect by CO2 collisional line mixing for the early Earth than for early Mars.

Fundamental changes to cloud properties on early Earth have been proposed as alternative solutions to the faint young Sun paradox. In particular, it has been proposed that cloud particles were larger on the early Earth because of less cloud condensation nuclei from biological sources [Rosing et al., 2010]. In figure 2, we tested the impact of larger liquid cloud particles (17 $\mu$m versus 12 $\mu$m). Larger particles warm the planet by decreasing the planetary albedo [Charnay et al., 2013]. This effect is particularly strong for temperate climate with a warming of around 8∘C for pCO2=0.1 bar. However, it becomes weak for warmer climates. We also simulated the climates with a fixed amount of CCN and varying cloud radii (black line in left panel in figure 2). The mean surface temperature is almost identical to the case with 12 $\mu$m particles.

The right panel in figure 2 shows the shortwave (albedo cooling effect) and longwave (greenhouse warming effect) radiative forcing by clouds depending on pCO2 for the different particle size parametrizations. The trends are similar for all particles size parametrizations. The net cloud radiative forcing is higher for cold climates (e.g. pCO2=0.01 bar) than for temperate climates (e.g. pCO2=0.1 bar). This is due to the lower albedo cooling of clouds over sea ice and to a reduction of the cloud cover over the cold equatorial ocean, particularly close to the freezing line (see Charnay et al. [2013]). The net cloud radiative forcing increases with pCO2 when pCO2 exceeds 0.1 bar. The albedo cooling effect decreases because of a reduction of lower clouds (see right panel in figure 3). At the same time, the greenhouse warming effect decreases. This is mostly due to the gaseous infrared emission level (the altitude where most of the thermal emission is emitted) going up with more CO2, which reduces the greenhouse effect of clouds. The cloud albedo cooling decreases faster than the cloud greenhouse warming, which leads to a net forcing increasing with pCO2. Clouds therefore produce a destabilizing feedback on the climate. For hot climates, they tend to have a small net radiative effect (their cooling effect being reduced by around 30 W/m2), which depends little on the radii of cloud particles.

In response to clouds properties, our GCM predicts that hot climates can be obtained with 1 bar or less of CO2 at 3.8 Ga. This required amount of CO2 is significantly lower than previous 1D estimations by Kasting and Howard [2006]. They found that if there were no methane, around 3 bars of CO2 would be required for getting a mean surface temperature of 60∘C at 3.3 Ga. This discrepancy is due to cloud feedbacks, not simulated in the 1D model, and to a few other differences in the climate modeling. In particular, they used the present-day land cover while we assumed no continent, and we used more recent opacity data, likely increasing the visible and infrared gaseous absorption. From a purely climate modeling point of view (i.e., ignoring long-term geochemical feedbacks), a warm early Earth with temperature of 60-70∘C is possible with sufficient CO2 if there is no land or almost no land.

3 Carbon cycle modeling

3.1 Model set-up and boundary conditions

In order to quantify the carbon cycle equilibrium response to simulated climate, we used the model from Le Hir et al. [2008a,  b]. We generally assumed a classical plate tectonic regime as in Sleep and Zahnle [2001] (see Fig. 1 in Appendix A), although there is debate about when a modern style of plate tectonics began [Stern, 2007]. CO2 reservoirs are 1) the atmosphere, 2) the surface/deep ocean, 3) the mantle, 4) carbonates on continental platforms and 5) calcite precipitation in the oceanic crust resulting from seafloor weathering. The inventory of carbon between these six reservoirs evolves through CO2 outgassing by mid oceanic ridges and arc volcanoes, continental and seafloor weathering, carbonate formation/dissolution, subduction and recycling. The part of subducted carbon released by arc volcanoes depends on the recycling factor. Based on a steady state assumption, we considered an instantaneous carbon recycling in subduction zones. The model parametrizations for these processes are described in Appendix A.

We simulated the carbon cycle using the surface temperatures and precipitation-evaporation rates (P-E) from the GCM. In the absence of constraints about geography, we assumed that emerged continents are homogeneously spread over all the Earth as small continents (each grid cell has its own fraction of continents). That tends to maximize the continental weathering. We interpolated GCM temperatures and Precipitation-Evaporation (P-E) between pCO2 values at each grid point. We parametrized the seafloor weathering with a dependence on the deep ocean temperature (see Coogan and Gillis [2013] for geological evidence and Krissansen-Totton and Catling [2017] for application for the last 100 Ma), assumed to be equal to sea surface temperature at 60*∘* latitude. Most of the weathering occuring in the upper part of the oceanic crust ($<$1km), the seafloor weathering is explicitly described between 0 and 500 m by 6 layers with an ambient temperature increasing with depth (116∘C/km). We assumed that the potential temperature of the mantle was around 1550∘C at 3.8 Ga compared to 1330∘C today. Formulas from Flament et al. [2008] (with an activation energy for the mantle viscosity of $\sim$350-400 kJ/mol) lead to a oceanic crust mass production rate $\sim$11 times higher at 3.8 Ga than today. This value is consistent with the estimation by Ohta et al. [1996] of a 10 times higher oceanic crust production rate during the Archean.

We assumed that the CO2 flux outgassed by oceanic ridges is proportional to the oceanic crust mass production rate while the seafloor weathering, limited to the first km of the oceanic crust, is proportional to the spreading rate. For a fixed oceanic crust mass production rate, a thicker crust has no impact on the outgassed CO2 flux but reduces the seafloor weathering in our model, leading to an enhancement of pCO2. It is generally accepted that the oceanic crust was thicker during the Archean ($\sim$18 km versus 6 km in average) [Bickle, 1986; Ohta et al., 1996; Flament et al., 2008]. If we assume that the early oceanic crust was 3 times thicker at the mid-oceanic ridge, the spreading rate (SR) must have been 3.8 times higher than today, according to our estimation of the oceanic crust mass production rate. Since the oceanic crust thickness and the associated spreading rate during the Archean are highly uncertain [Bickle, 1986], we therefore tested cases with the present-day oceanic crust thickness (6 km) and SR=11.3.

For a warmer mantle and a higher spreading rate on the early Earth, the younger oceanic lithosphere was warmer and melted more easily when subducted. If arc volcanoes, associated with subduction zones, provide a quick return of materials to the surface, the direct melting of oceanic crust may increase the amount of fluid released, suggesting a recycling of suducted carbon higher than today. According to this assumption, the decarbonation efficiency may exceed Cenozoic values used as reference [Hoareau et al., 2015]. We thus explored the effect of the recycling (R) using the Cenozoic value as reference (R=0.1, corresponding to a recycling efficiency of 10$\%$) and values 3 times higher (R=0.3, a recycling efficiency of 30$\%$).

The fraction of emerged land during the Hadean/Archean is still debated [Roberts and Spencer, 2015]. We made the assumption that this fraction was low ($<$30$\%$ of the present-day fraction) during the early Archean as suggested by some geochemical data (e.g. [Dhuime et al., 2012]), some evidence for a more voluminous ocean during the Archean [Pope et al., 2012] and some models (e.g. [Flament et al., 2008]). This hypothesis allows us to use the outputs from the 3D land-free simulations for the carbon cycle model. We explored the effect of the emerged continental fraction considering four cases:

1. an aquaplanet (no land)

2. a very low emerged land fraction (1$\%$ of the present-day land fraction)

3. a moderate continental fraction (12$\%$ of the present-day land fraction)

4. a relatively high continental fraction (26$\%$ of the present-day fraction)

Although there is clear evidence of emerged land at 3.7-3.8 Ga [Nutman et al., 2015], an aquaplanet model is an idealized case, which is very useful for understanding the effect of seafloor weathering.

3.2 Results and discussion

We ran the climate-carbon model for the four land fractions (0, 0.01, 0.12 and 0.26 times present-day fraction), with a recycling factor R equal to 1 (present-day) or 3, with or without 2 mbars of atmospheric methane and with the present-day oceanic crust thickness (H=6 km, SR=11.3) or a 3 times thicker oceanic crust (H=18 km, SR=3.8). The equilibrium points span a temperature range from 8.5∘C to 31.5∘C and a pCO2 range from 0.12 bar to 0.36 bar (see figure 4 and table 1). Ice-free or partly ice-free oceans are obtained in any cases.

3.2.1 Cases with no land

If there is no emerged land, the carbon cycle is controlled by the seafloor weathering and volcanic/mid-oceanic ridge outgassing. At equilibrium, the seafloor weathering is simply given by $F_{sfw}=\frac{F_{mor}+F_{deep}}{1-R}$ ($F_{mor}$ and $F_{deep}$ are the CO2 fluxes from midocean ridges and mantle, respectively, see Appendix A) and depends on $R$. Because the carbon species partitioning between ocean and atmosphere is controlled by temperature, sea surface temperature (SST) implicitly determines the seawater pH and pCO2 for a fixed value of seafloor weathering rate. Light blue and purple curves in figure 4 reveal the CO2 partial pressure for a given SST with our model. They are computed assuming no CO2 climate feedback and thus represent the intrinsic response of the carbon model for Earth’s climate in the absence of continents. The light blue curves (6 km thick crust) reveal a maximum pCO2 of 0.14 bar for R=0.1 (left panel) and 0.26 bar for R=0.3 (right panel). The CO2 solubility in the ocean decreases with temperature. At temperatures lower than 15∘C, this feedback dominates and pCO2 increases with temperature. When the temperature becomes higher than 15∘C, the temperature dependance of the seafloor weathering dominates and pCO2 decreases with temperature. The combination of both effects (solubility and seafloor weathering) leads to an upper limit for pCO2 for the thin oceanic crust case. The reduced seafloor accretion rate associated to a thicker oceanic crust scenario combined with cold water strongly limits carbon consumption by seafloor weathering, hence the purple curves (18 km thick crust) do not show a maximum pCO2 below 1 bar. For these cases without land, the equilibrium point of the climate-carbon model corresponds to the crossing between the intrinsic response of the carbon model (i.e. the light blue and purple curves) and the GCM curves showing the global mean surface temperature as a function of pCO2 (i.e. the blue and red curves). For instance, the equilibrium point for a methane-free atmosphere with R=0.1 and a thick crust corresponds to a mean SST of 22.3∘C and pCO2=304 mbar.

3.2.2 Effects of land cover

According to our model, land cover has only a limited impact on the carbon cycle and the global climate for the early Earth (temperature change of $\sim$1∘C between the four land distributions). Indeed, while the carbon cycle was mostly controlled by continental weathering during the Proterozoic and the Phanerozoic (e.g. Krissansen-Totton and Catling [2017]), it was mostly controlled by seafloor weathering on the early Earth. That difference is mainly due to the higher spreading rate (high accretion rate provides a sufficient quantity of weatherable materials which enhances the seafloor weathering) and to the assumed lower land fraction on the early Earth. Under such circumstances, our assumption about the land distribution (i.e. small continents distributed homogeneously) is therefore of minor importance for the results. We also notice that the feedback produced by land does not have the same sign with a 6 or 18 km thick oceanic crust. With a 6 km oceanic crust, the seafloor weathering is slightly stronger than the continental weathering per unit of surface. In that case, the global mean surface temperature and pCO2 are higher with a higher land fraction. In contrast, with a 18 km thick oceanic crust, the continental weathering is stronger than the seafloor weathering per unit of surface (as today) and a higher land fraction leads to lower pCO2. The coldest case with our assumptions (mean SST of 8.5∘C) is thus obtained with no land, a methane-free atmosphere, R=0.1 and a thin crust. The warmest case (mean SST of 31.5∘C) is obtained with no land, a methane-rich atmosphere, R=0.3 and a thick oceanic crust, as shown by the blue point in Figure 4b.

3.2.3 Effects of methane, oceanic crust thickness and recycling rate

A thicker oceanic crust or a higher recycling efficiency always leads to an increase of pCO2 and so to an increase of temperature. In contrast, adding methane in the atmosphere reduces pCO2 but still increases the global mean surface temperature. If we use as a reference, the case with a moderate land cover, a 6 km thick crust, no methane and R=0.1 (mean surface temperature=8.5∘C), adding 2 mbar of methane produces a 11.7∘C warming, a 18 km crust produces a 12.9∘C warming, and a 3 times higher recycling rate produces a 8.8∘C warming. These three parameters therefore produce significant temperature and pCO2 changes. In particular, the thicker oceanic crust is generally expected for the early Earth, with an effect on pCO2. A thicker oceanic crust would have had a large impact, allowing warmer climate for the early Earth than today. The combination of a thicker crust, a high recycling rate, and the presence of methane leads to global mean surface temperatures around 30∘C for the three land distributions. That temperature change ($\sim$20∘C) is much lower than the sum of individual warming ($\sim$33∘C), mostly because of the non-linear dependence of the seafloor weathering with the deep ocean temperature.

3.2.4 Absence of plate tectonic: the heat-pipe scenario

The onset time of plate tectonic is unknown but most geological and isotopic data suggest that subduction started during the Hadean or early Archean [Sleep, 2007; Tang et al., 2016]. Moore and Webb [2013] suggested that before the onset of plate tectonic, Earth was in a heat-pipe regime in which volcanism dominated surface heat transport. We tested the consequences of this possible early regime for the carbon cycle with our model. We considered that all the heat flux is released by volcanism, that the spreading rate is null and that submarine basalt flows are weathered as seafloor basalts. In our simulation of the heat-pipe Earth, we assumed a very low land fraction (1$\%$ of the present-day land cover). The volcanic outgassing becomes around 4 times as high as with the plate tectonic regime, or 10 times as high as the present volcanic outgassing. At equilibrium, pCO2 is 0.26 bar and the mean surface temperature is around 18.8∘C. A heat-pipe regime would thus have led to a temperate early Earth.

3.2.5 Discussion

Our climate-carbon model never reaches warm climates for the early Earth because of the strong temperature feedback by the seafloor weathering. It predicts global mean surface temperatures around 8.5-30.5∘C, with pCO2 around 0.1-0.36 bar. This temperature range is compatible with the most recent measurements by Hren et al. [2009] and Blake et al. [2010], and also with the evolutionary model from Bousseau et al. [2008] for a mesophilic last universal common ancestor (LUCA). Thermophily near the roots of the tree of life may reflect local warm environments (like hydrothermal vents) or survival after hot climates produced just after large impacts [Bousseau et al., 2008; Abramov and Mojzsis, 2009]. The pCO2 range of 0.1-0.35 bar is generally higher than the upper limits initially derived from the stability of different minerals from Archean samples [Rye et al., 1995; Sheldon, 2006; Rosing et al., 2010; Driese et al., 2011]. These constraints globally suggest a pCO2 lower than 30 mbars. However, they are not all compatible together and were essentially obtained for the late Archean. Recent work by Kanzaki and Murakami [2015] gives an upper limit for pCO2 of 0.14-0.7 bar at 2.77-2.46 Ga, compatible with our results.

4 Climate and carbon cycle during the Late Heavy Bombardment

4.1 Impact modeling

The Earth was likely struck by numerous of impactors during the Late Heavy Bombardment (LHB), which would have occurred between around 4.1 and 3.8 Ga [Gomes et al., 2005; Bottke and Norman, 2017]. These impactors may have been a sink of CO2 by producing large amounts of ejecta, which can be easily weathered [Sleep and Zahnle, 2001]. They may also have been a source of CO2 by vaporizing carbonates from the oceanic crust. For instance, Chixulub impact may have outgassed 350-3500 Gt of CO2 [Pierazzo et al., 1998], or the equivalent of 0.06-0.6 mbar of CO2. We investigated the long term effects of impacts on the carbon cycle to understand their consequences for the mean climate during the LHB, in particular if they had a global net warming or cooling effect.

We assumed that the frequency of impact follows a scaling law given in Anbar et al. [2001] and Sleep and Zahnle [2001]:

[TABLE]

where f is the frequency of impactors having a mass higher than $m$. For the LHB, we chose $b$ = 0.7 [Anbar et al., 2001]. We considered that the maximal mass was $m_{max}=10^{19}$ kg (i.e. diameter around 190 km). This upper limit is constrained by the largest crater observed on the Moon. We neglected impacts with a mass lower than 1013kg (i.e. diameter lower than around 1.8 km), which produce little ejecta because of the strong deceleration by the ocean. We fixed $f_{0}$ = 2.9$\times 10^{-8}$ impacts/year in order to have a cumulative impactor mass equal to 2$\times 10^{20}$ kg for the 300 millions of years of the LHB. That corresponds to an impact flux of 6.7$\times 10^{11}$ kg/year.

As Sleep and Zahnle [2001] and Zahnle and Sleep [2002], we computed the weathering of ejecta, which provides a source of cations, precipitating as carbonates in the ocean. We assumed that the mass distribution of ejecta is given by $\it N(\geq m)\propto m^{-\gamma}$, where $\gamma$=0.87 (derived from scaling laws of Collins et al. [2005]). Zahnle and Sleep [2002] used a similar law with $\gamma$=0.9-0.95. We considered that the diameter of each particle decreases of 1 mm in 167 kyrs [Crovisier et al., 1987] and that cations are present at 0.005 mol/g. With the oceanic crust lifetime for the early Earth, ejecta with a diameter larger than 100 m contribute little to the weathering compared to smaller ejecta.

We included CO2 delivered by the impactors as done by de Niem et al. [2012]. We assumed that the mass fraction of CO2 is 1$\%$ in asteroids and 20$\%$ in comets. For comparison, Rosetta measurements of Comet 67P / Churyumov-Gerasimenko indicate a mass ratio of CO2 to H2O around 40$\%$ [Marty et al., 2016] for a mass fraction of H2O to dust around 20$\%$, implying a mass fraction of CO2 around 8$\%$. If all impactors were comets, the LHB would have provided the equivalent of 3.8 bar of CO2 (taking into account the escaping mass of impactors). The fraction of comets during the LHB is not known and Bottke et al. [2012] suggest that comets were a minor player. We did simulations using fractions of comets of 0 and 100$\%$. We also included the CO2 outgassed by oceanic carbonates and lithosphere melting after impact. We considered that the oceanic crust is covered by a 40 m layer of carbonate, which outgassed CO2 when shocked at a high pressure (above 60 GPa). We assumed that the lithosphere below the carbonate layer (considered as pure calcite) contains 100 ppmv of CO2 and outgasses all CO2 when melted. The complete assumptions for the ejecta weathering and CO2 outgassing are described in Appendix B.

Figure 5 shows the amount of CO2 outgassed from impactor/crust and consumed by ejecta as a function of impactor mass for asteroids (density=3000 kg/m2 and initial velocity=21 km/s) and comets (density=1000 kg/m2 and initial velocity=30 km/s). The solid red lines show the maximal amount of CO2 that can be consumed by ejecta. The solid red line show the amount of CO2 consumed by ejecta in 15 Ma (the lifetime of the oceanic crust in our model). With our size distribution for ejecta, most of ejecta are weathered during the lifetime of the oceanic crust. In all cases, impacts produce a net sink of CO2 by ejecta weathering.

4.2 Results and discussion for the LHB

To simulate the effect of impacts during the LHB, we computed sequences of 10 millions of years with stochastic impacts following the statistical distribution described above. For simplicity, we assumed that impacts never overlap and that they all fall into the ocean. This assumption is valid because we consider small continental fractions and because the total surface of transient craters (diameter $D_{tc}$ in Appendix B) accumulated in 10 million years only represents 0.6 $\%$ of the Earth’s surface. These sequences provide a time-dependent CO2 outgassing rate and alkaline source, which are used in the carbon cycle model to computed the climate and pCO2 evolution. The top panel in figure 6 shows the cumulated CO2 outgassed by impacts and consumed by ejecta during a typical 10 million year sequence with stochastic impacts. In this sequence, one 143 km impactor strikes the Earth at 7 Ma and dominates CO2 fluxes. On average, the weathering of ejecta produces a cations flux of 16 $\times 10^{12}$ mol/yr (compared to 300 $\times 10^{12}$ mol/yr in Sleep and Zahnle [2001]), around half the CO2 flux from mid oceanic ridges and arc volcanoes (around 30$\times$1012mol/year). The vaporization of carbonates produces a CO2 flux of 0.05$\times$1012mol/year, negligible on the long term compared to the other fluxes. The bottom panel in figure 6 shows the evolution of pCO2 during this sequence for the 0.01 land fraction, no methane, R=0.1, a thin or thick oceanic crust and with asteroids or comets. The giant impact at 7 Ma produces a strong decrease of pCO2 by ejecta weathering triggering a full glaciation for thin oceanic crust cases and a very cold climate with a global mean surface temperature of -10∘C for the thick oceanic crust case. Before this impact, the carbon cycle is stabilized with pCO2 around 50 mbar for the thin crust and around 250 mbar for the thick crust. As in Sleep and Zahnle [2001] and Zahnle and Sleep [2002], we thus find a strong cooling by ejecta weathering leading to cold climates, but not always completely ice covered. In particular, cases with a thick crust, methane and R=0.3 should remain temperate or at least partially ice-free most of the time during the LHB.

To better understand the role of each parameter, we derived a simple analytical expression of pCO2. If we assume a case with no land, a recycling R0 and an oceanic crust thickness H0, the equilibrium of the carbon cycle can be expressed for a reference state with no impact as:

[TABLE]

where $F_{mor}^{0}$, $F_{arc}^{0}$ and $F_{sfw}^{0}$ are the CO2 fluxes from mid-oceanic ridges, arc volcanoes and seafloor weathering from the (crust+mantle) to the (ocean+atmosphere). For a new state with other values of R and spreading rate and with impacts, the equilibrium is given by:

[TABLE]

with $F_{mor}=F_{mor}^{0}$ and $F_{arc}=F_{deep}-R\times(F_{sfw}+F_{ej})$ (see appendix A). We make the assumption that the seafloor weathering does not depend on oceanic temperature, what is valid for cold climates. By using the simple expression: $F_{sfw}=F_{sfw}^{0}\left(\frac{pCO_{2}}{pCO_{2}^{0}}\right)^{\alpha}\left(\frac{SR}{SR_{0}}\right)$ with $\alpha$=0.5 [Haqq-Misra et al., 2016], we obtain an expression for pCO2:

[TABLE]

The term $\frac{(1-R)\times F_{ej}}{F_{mor}+F_{deep}}$ simply corresponds to the ratio of the amount of CO2 lost in the mantle (due to ejecta weathering) by the amount of CO2 outgassed directly from the mantle (e.g. without taking into account recycling). When this ratio is higher than 1, there is a net sink of CO2 toward the mantle and the Earth necessary falls into a full snowball Earth. That occurs for $F_{ej}\leqslant$ -29$\times$12 mol/yr for R=0.1 and $F_{ej}\leqslant$ -37$\times$12 mol/yr for R= 0.3 with our model values. With our impact distribution, the ejecta cations flux (16 $\times$ 1012 mol/yr) is below these limits, and the Earth should have avoided a permanent full glaciation during the LHB. Because of the large uncertainties on the carbon cycle parameters and the duration and intensity of the LHB, it is not possible to be definitive on that point. We conclude that the LHB likely had a strong impact on the carbon cycle and climate with potentially snowball Earth events after large impacts.

Large impactors were also frequent during all the Archean [Bottke et al., 2012; Johnson and Melosh, 2012; Bottke and Norman, 2017] but much less than during the LHB. According to the dynamical model from Bottke et al. [2012], the impact frequency was an order of magnitude lower at 3.7 Ga than at the peak of the LHB. Therefore, impacts probably had only a limited effect on the long term carbon cycle and climate during most of the Archean. However, our LHB impact sequence reveals that individual giant impacts could trigger episodic glaciations. Johnson and Melosh [2012] evaluated the diameters of ancient large impactors from spherule layers during the Archean and the Proterozoic. A few of them with diameters higher than around 50 km could have been large enough to trigger glaciations.

The existence of the Late Heavy Bombardment is still debated by some authors. For instance, Boehnke and Harrison [2016] suggest that the LHB could be an artifact of monotonically decreasing impact flux combined with episodic crust formation. If it did not occur, the impact flux at the end of the Hadean would have been too small to affect the climate on the long term, but strong cooling after episodic large impacts may still have been possible.

5 Conclusions

Our 3D climate simulations (taken in isolation) reveal that warm climates, with mean surface temperatures around 60∘C, could have been obtained during the late Hadean/early Archean with around 1 bar of CO2. However, our carbon model suggests that such states would not have been stable. Indeed, we found that warm climates can be maintained only if the seafloor weathering is assumed independent of temperature. Such an assumption is inconsistent with recent measurements indicating a higher weathering during the late Mesozoic, when the deep ocean was around 10∘C warmer [Coogan and Gillis, 2013]. Stable very warm climates with a N2-CO2-CH4-H2O early atmosphere are therefore excluded.

Our results favor instead temperate climates with global mean surface temperatures around 8.5-30.5∘C and with pCO2 around 0.1-0.36 bar. With the increasing solar luminosity over the Archean, the carbon cycle would have led to a reduction of pCO2. In parallel, with the emergence of land and the slowing seafloor spreading rate during the Archean, the carbon cycle would have progressively become controlled by continental weathering while initially controlled by seafloor weathering. However, a better understanding of the data from O isotopes in ancient cherts together with better constraints on pCO2 are required to validate our conclusions about the long-debated climate of early Earth.

During the Late Heavy Bombardment, the weathering of ejecta would have dramatically decreased the partial pressure of CO2 leading to cold climates but not necessary to a snowball Earth during all that period. We predict that large impactors (i.e. with diameters higher than around 50 km) could have trigger extreme glaciations during the Hadean and the Archean.

Finally, our modeling work suggests that the CO2-dependent carbon cycle could have maintained temperate climates on the early Earth, without requiring any additional greenhouse gas or warming process. This work also highlights the importance of the seafloor weathering, which can dominate the continental weathering. Under such circumstances, the climate becomes relatively insensitive to the presence of land. The strong temperature dependence expected for the seafloor weathering also has implications for the habitability of exoplanets, allowing in particular a stabilizing climate-weathering feedback on waterworlds. Waterworlds might therefore not be much more subject to water loss by moist greenhouse than partially ocean-covered Earth-like planets, as suggested by Abbot et al. [2012].

Acknowledgments:

B.C. acknowledges support from Paris Sciences et Lettres and from an appointment to the NASA Postdoctoral Program, administered by USRA. This work was performed as part of the NASA Astrobiology Institute’s Virtual Planetary Laboratory, supported by NASA under Cooperative Agreement No. NNA13AA93A. DC also acknowledges support from NASA Exobiology Program grant NNX15AL23G. We thank J. Kasting and an anonymous reviewer for helpful comments.

Appendix A Carbon cycle modeling

The carbon cycle model has been previously used for the Proterozoic and the Phanerozoic [Goddéris et al., 2006; Le Hir et al., 2008a,  b]. Figure 7 shows the different processes included in the model. The parametrizations for the different CO2 flux are briefly described below.

CO2 outgassing by oceanic ridges and arc volcanoes

We considered that the CO2 flux from mid oceanic ridges $F_{mor}$ is proportional to the production of oceanic crust:

[TABLE]

where $F_{mor_{0}}$,= 1.8$\times$1012 mol/yr, $SR$ is the spreading rate and $H$ is the oceanic crust thickness ($H_{0}$ = 6 km).

We assumed that the flux from arc volcanoes (Farc) is the total degassing flux due to carbonate recycling at subduction zones and mantle outgassing (with the exception of mid-oceanic ridge). Degassing over geologic times being assumed proportional to the heat flux from the interior of the Earth, the mantle outgassing rate is adjusted, as appropriate, using a forcing factor, rdeg=2.5 at 3.8 Ga [Goddéris and Veizer, 2000]):

[TABLE]

where $R$ is the recycling rate and $F_{deep}=C_{\rm deep}^{0}\times rdeg$ with $C_{deep}^{0}$ = 3.3 $\times$ 1012 mol/yr. For the heat-pipe case, we considered that $SR$ = 0, $F_{\rm sub}$ = 0 and $rdeg$ =10. Indeed, according to Moore and Webb [2013], the heat flux had to be around 10 times higher than today to be in the heat-pipe regime.

Continental weathering

Weathering depends on rocks available for reactions (basalt vs. granite) and fraction of continents deduced from crustal growth models. The weathering rate of continental granites and basalts is given by:

[TABLE]

with:

[TABLE]

with $fgra=\frac{-48200}{8.314}\left(\frac{1}{T}-\frac{1}{T_{0}}\right)$, $fbas=\frac{-42300}{8.314}\left(\frac{1}{T}-\frac{1}{T_{0}}\right)$ and where $\frac{srf_{\rm cont}}{srf_{\rm cont0}}$ is the fraction of emerged continental area.

We used $pCO_{2}^{0}$ = 285 ppm. These values give a present-day flux of 3.49$\times$1012 mol/yr for Ca2+ and 3.31$\times$1012 mol/yr for Mg2+.

To represent the presence of very active volcanoes and potential extent of basaltic flows over small continents, rocks available for weathering is assumed proportional to volcanism ( $F_{\rm arc}$). Based on the extent of basalt outcrops of the Cenozoic era, the equilibrium between basaltic flows / volcanism is prescribed as follow:

[TABLE]

with $F_{\rm arc0}=6.8\times 10^{12}$ mol/yr [Gaillardet et al., 1999], $srf_{cont0}=146\times 10^{6}$ km2 and $srf_{bas0}=7.3\times 10^{6}$ km2 [Dessert et al., 2003].

Using this assumption, our model predicts that early Archean continents were mostly mafic, composed   65-70$\%$ of basalts and 35-30$\%$ of granite, consistent with analysis of Sb/Sr, Ni/Co and Cr/Zn ratios by Dhuime et al. [2015]; Tang et al. [2016].

The runoff and the surface temperature depend on the latitude and are properly computed by the GCM. We assumed that lands are homogenously distributed over all the Earth. Then, we computed the global continental weathering flux by summing the contribution of each latitudinal band from the GCM.

Seafloor weathering

Contrary to Caldeira [1995], we assume that the seafloor weathering does not occur in solutions at steady state with carbonate minerals, but when oceanic water enters in contact with fresh basaltic rock, it evolves along its downward path towards saturation with calcite, which is the ultimate state of the reaction. Our percolating reservoir stands for water during the early stages of seafloor weathering, when all cations are produced. The seafloor weathering rate is given by [Goddéris et al., 2006; Le Hir et al., 2008a,  b]:

[TABLE]

where $SR$ is the spreading rate. The sum extends up to all the dissolving species (H+, OH– and H2O), while the index $i$ stands for the dissolving species. $E_{i}$ and $d_{i}$ are, respectively, the activation energy and dissolution constant for each mineral depending on the species $i$ promoted dissolution. $T_{p}$ is the temperature of the percolating waters into the oceanic crust at which dissolution occurs. We divided the oceanic crust in 6 layers from the seafloor to a depth of 500 m in the crust. We assumed that the percolating water has the temperature of surrounding crust fixed with the present-day temperature gradient (around 116 K/km). For the first layer, the deep ocean temperature is assumed to be equal to the sea surface temperature at 60*∘* N/S (see idealized simulations by Manabe and Bryan [1985]; Enderton and Marshall [2009]). We scaled the seafloor weathering to its present-day value. $C_{i}$ are the concentration in the species $i$, while $n_{i}$ is the order of the dissolution reaction. We considered that the seafloor weathering occurs only in the first 500 m of the oceanic crust. A change of oceanic crust thickness thus does not impact the seafloor weathering, which remains proportional to the spreading rate.

Experimental studies and measurements of carbonate mineral in the oceanic crust reveal that the seafloor weathering strongly depend on the bottom ocean temperature and was around 5 times higher during the late Mesozoic [Brady and Gíslason, 1997; Coogan and Gillis, 2013], when the deep ocean temperature was around $\sim$10 K higher. It is therefore expected to play a significant role on the early Earth with a lower land cover and a faster spreading rate. Our model predict a 2-3 times higher seafloor weathering during the late Mesozoic for 4 times the present-day pCO2. However, the volcanic outgassing may have been higher and continental weathering may have been less efficient at that time, what would have increased the seafloor weathering closer to the measured rates. Because of these uncertainties, we consider that our values remain acceptable but we may underestimate the temperature sensitivity of the seafloor weathering.

Precipitation of carbonates

In the absence of pelagic producers, carbonates are assumed to be precipitated on continental margins (at a depth of less than 100 m) through:

[TABLE]

where $k_{cd}$ = 0.172 mol/yr/m2 and $\Omega$ is the saturation ratio of calcite: $\Omega$ = $\frac{[Ca^{2+}][CO_{3}^{2-}]}{K_{s}}$ where $K_{s}$ is the solubility constant depending on the temperature and the salinity (expression given in Sacramento and Gruber [2006]).

Ocean-atmosphere interface

The CO2 exchange flux between the atmosphere and the oceanis given by:

[TABLE]

where $\alpha$ = 16 mol/m2/PAL and $pCO_{2}^{ocean}$ is the CO2 partial pressure computed with Henry’s law for the concentration of dissolved CO2.

Seafloor spreading rate

We computed the seafloor spreading rate $SR$ considering it evolves with the potential temperature of the mantle $T_{p}$ as Rayleigh-Bénard convection [Flament et al., 2008]:

[TABLE]

where $T_{0}$ is the deep ocean temperature (we assumed $T_{0}$ = 0∘C, but the value has a very small impact on $SR$), $\eta(T_{p})=\eta_{0}\left(\frac{E}{RT_{p}}\right)$ is the viscosity of the mantle with $E$ = 380 kJ and $R$ = 8.31 J/K/mol. Assuming $T_{p}^{*}$=1330*∘*C for the present-day Earth and $T_{p}$=1550∘C for the early Earth then gives $u/u_{*}$ = 11. If the oceanic crust thickness is increased by a factor 3, then the spreading rate becomes 3 times lower while the crust production rate and the heat flux remain constant. A spreading rate 11 times higher than today gives a lifetime of the oceanic crust of around 15 Ma.

Appendix B Effects of impacts during the Late Heavy Bombardment

Ejecta weathering

We computed the effect of ejecta weathering following Sleep and Zahnle [2001]. Impacts produce a global layer of particles with diameters typically from 1 mm to 100 m. Ejecta constitute a source of cations that precipitate with carbonate ions. The mass of ejecta was determined using the scaling laws from Collins et al. [2005]. We took into account the deceleration caused by the atmosphere and the ocean but we neglected the breakup in the atmosphere and the mass loss by sublimation. These assumptions are valid for large impactors. We also included the escaping mass of impactor after the impact according to the hydrodynamic simulations from Shuvalov [2009].

We assumed that impactor have a velocity $v_{0}$ of 21 km/s for asteroids and 30 km/s for comets, a diameter $L$, an impact angle $\theta$ of 45∘ and a density $\rho_{i}$ = 3000 kg /m3 for asteroids and $\rho_{i}$ = 1000 kg /m3 for comets. The Earth is supposed covered by a 2700 m deep global ocean. We chose a density $\rho_{c}$ = 2700 kg/m3 for the crust.

The mass of ejecta is computed according to Collins et al. [2005]:

[TABLE]

where the transient crater $D_{tc}$ (in m) is given by [Collins et al., 2005]:

[TABLE]

where $g_{E}$ is Earth gravity and $v_{i}$ is the impactor velocity at the seafloor.

Taking into account the deceleration by the atmosphere and the ocean, the impactor velocity at the seafloor is [Collins et al., 2005]:

[TABLE]

where $C_{D}$ is the drag coefficient, taken equal to 2 for the atmosphere and 0.877 for the ocean, $H_{atm}$ is the atmospheric scale height, and $H_{ocean}$ is the ocean depth. For a 10 km diameter impactor, the atmosphere reduces the velocity by just 0.1$\%$ but the ocean reduces it by around 8$\%$.

Collins et al. [2005] suggest that the mean diameter of ejecta decreases with distance from the crater center $r$ as $d\propto r^{-\alpha}$, where $\alpha$=2.65 and that the thickness of the ejecta layer decreases as $h\propto r^{-3}$. By integrating these relations, the mass distribution of ejecta is given by $\it N(\geq m)\propto m^{-\gamma}$, where $\gamma$=0.87. Zahnle and Sleep [2002] use the same law with $\gamma$=0.9-0.95. We assumed that the diameter of each particle decreases of 1 mm in 167 kyrs [Crovisier et al., 1987] and that cations are present at 0.005 mol/g.

With the impact statistics that we use, the impact flux is around $6.7\times 10^{11}$ kg/yr and the ejecta production rate is around $3.2\times 10^{12}$ kg/yr, corresponding to a cation flux of $1.6\times 10^{13}$ mol/yr (compared to 3 $\times 10^{14}$ mol/yr in Sleep and Zahnle [2001]). In our carbon cycle model, these cations precipitate if the saturation ratio of calcite is higher than 1 (see Appendix A).

CO2 outgassing by impacts

We assumed that the early oceanic crust is covered by a 40 m thick carbonate layer ($\rho_{carb}$ = 1500 kg/m3). We therefore assumed the equivalent of 4 bars of CO2 ($4.8\times 10^{20}$ moles of CO2) are in this layer and can be released during impacts. This amount is lower than the total quantity of CO2 currently stored in the oceanic crust (around $1.2\times 10^{21}$ moles of CO2 [Sleep and Zahnle, 2001]). The shock wave produced by an impact can lead to decarbonation of carbonates. For calcite, the reaction is CaCO3 $\Rightarrow$ CaO + CO2. This process starts when the pressure is higher than 20 GPa and is total when the pressure is higher than 60 GPa [Agrinier et al., 2001]. However, the amount of outgassed CO2 is reduced by recombination through the back reaction. The percentage of CO2 back reacted is between 37 and 66 $\%$ [Agrinier et al., 2001]. For simplicity, we assumed that CO2 outgassing occurs only at pressures higher than 60 GPa with an efficiency of 50 $\%$ due to the reverse reaction. For each impact, we computed the area where the shock pressure is higher than 60 GPa, assuming that it evolves with the distance $r$ to the impact center as $\propto r^{-1.5}$.

We also took into account CO2 outgassed by melting lithosphere, assuming that its volume $V_{m}$ is [Collins et al., 2005]:

[TABLE]

where $E$ is the impactor kinetic energy and $\theta$=45∘ is the impact angle. We assumed that the lithosphere below the carbonate layer contains 100 ppmv of CO2 and outgasses all CO2 when melted.