© ESO 2019

1 Introduction

The negative feedback provided by the long-term carbonate–silicate cycle is important for keeping the climate of Earth habitable in the long-term (Walker et al. 1981; Kasting & Catling 2003). Should plate tectonics on planets beyond our solar system work in a similar way to as they do on Earth, these planets could remain habitable over a wide range of orbital distances and long time-spans (Kasting et al. 1993; Franck et al. 2000). The negative feedback is mainly controlled by temperature-dependent continental weathering, which requires emerged land (Foley 2015). On the one hand, emerged land allows for a strong temperature dependence of weathering, and on the other hand, it allows for efficient erosion (e.g. Flament et al. 2008), exposing fresh uncarbonated surface for weathering. However, the continental volume on planets beyond our solar system may differ significantly (Höning et al. 2019). Furthermore, in the early evolution of Earth, the continents were largely flooded, resulting in a strongly reduced erosion rate (e.g. Flament et al. 2008), affecting the efficiency of continental weathering.

The fraction of emerged land required to keep the continental weathering feedback efficient may be small, and Abbot et al. (2012) argue for a strong continental weathering feedback for a continental fraction of more than 1%. However, without any emerged land, this particular weathering feedback mechanism does not work. A negative feedback mechanism that works for water-covered oceanic crust is provided by seafloor weathering (Brady & Gislason 1997; Sleep & Zahnle 2001). Seafloor weathering was more important in the history of Earth, when the surface was covered by water, than it is today (Krissansen-Totton & Catling 2017). Mills et al. (2014) argue for a shift in importance from seafloor to continental weathering 1.5–0.5 billion years ago and Coogan & Gillis (2013) argue that both feedbacks had a similar magnitude in the late Mesozoic. The major difference with respect to continental weathering is that seafloor weathering does not depend on the erosion rate to provide fresh surface that can be weathered, but rather on the supply rate of fresh basaltic crust by volcanism. On planets with plate tectonics, fresh basaltic crust is mainly generated by seafloor spreading. The way seafloor weathering behaves is debated. Seafloor weathering was commonly believed to directly depend on the atmospheric CO₂ partial pressure (Brady & Gislason 1997; Sleep & Zahnle 2001). However, more recent studies argue that it mainly depends on the pore-space temperature, which in turn is controlled by the planetary surface temperature (Coogan & Dosso 2015; Krissansen-Totton & Catling 2017).

Stagnant-lid planets lack mid-ocean ridges and subduction zones. The long-term carbonate–silicate cycle operates differently for these bodies than it does for planets with plate tectonics. Volcanic degassing occurs at hot spots, and the accumulation of CO₂ in the atmosphere affects the habitability of these planets (Noack et al. 2017; Tosi et al. 2017; Godolt et al. 2019). Tosi et al. (2017) showed that the oxidation state of the mantle crucially affects the amount of degassed CO₂, and the accumulated atmospheric CO₂ affects the outer boundary of the habitable zone. Large amounts of atmospheric CO₂ allow liquid water to exist on the planetary surface farther away from the star. The inner boundary of the habitable zone is also affected by the atmospheric CO₂ if the amount of degassed H₂O is small (Tosi et al. 2017). However, it remains unclear as to whether or not carbon recycling can occur on stagnant-lid planets and to what extent this can regulate their atmospheric CO₂ partial pressure.

The production of fresh crust by volcanism provides a potential carbon sink, which could establish a negative feedback similar to seafloor weathering (Foley & Smye 2018). Volcanic eruptions may bury carbonated crust, which then sinks deeper into the mantle. Increasing pressure and temperature with depth affect the stability of carbonates. The release of CO₂ back into the atmosphere caused by the increase of temperature with depth is known as metamorphic decarbonation (e.g. Bickle 1996). Accounting for burial and decarbonation of the crust, Foley & Smye (2018) showed that for stagnant-lid planets with an Earth-like total CO₂ budget, carbon cycling could provide negative feedback avoiding a supply-limited regime, which would occur through a complete carbonation of the crust. However, we still do not know whether this negative feedback is sufficiently strong to regulate the climate, or if CO₂ rather accumulates in the atmosphere with time, with weathering hardly influencing the atmospheric CO₂ budget.

On planets without plate tectonics, the production of fresh crust occurs via hot-spot volcanism. In the presence of water, the basaltic crust can be carbonated and subsequently buried by new volcanic eruptions. This also implies that even if subsequent erosion of the upper crust were to take place, it would not expose uncarbonated crust for weathering, in contrast to the present-day Earth. Instead, the crust exposed by erosion was already carbonated at the time it was formed. Therefore, the CO₂ sink on a stagnant-lid planet is directly coupled to the production rate of new crust and hence to the CO₂ source and not to the erosion rate. Valencia et al. (2018) studied the seafloor-weathering feedback for tidally heated planets. These authors calculated equilibrium states of the atmosphere and mantle and concluded that the negative feedback may keep these planets habitable if fresh rock is supplied at a sufficiently large rate and seafloor weathering is temperature dependent. However, the mantle temperature of a planet that is not tidally heated substantially evolves with time. This is particularly important for the stability of carbonates, which can enter the mantle in substantial quantities only if the temperature–depth gradient is shallow (Foley & Smye 2018). Foley (2019) explored the time-span over which a stagnant-lid planet can sustain active volcanism and therefore a temperate climate, depending on the planetary CO₂ budget and on the internal heating rate. In particular, he showed that the minimum CO₂ budget of a planet required to recover from a snowball state crucially depends on whether or not exchange of CO₂ between the ocean and the atmosphere is possible.

In this paper, we focus on the interplay between the mantle temperature and the atmospheric CO₂ during planetary evolution and explore factors controlling the accumulation of CO₂ in the atmosphere or climate regulation via seafloor weathering. In particular, we study the influence of different relationships for seafloor weathering: models where seafloor weathering is directly dependent on the CO₂ partial pressure (Sleep & Zahnle 2001) and models where seafloor weathering is dependent on the pore-space temperature, which is a function of the surface temperature (Krissansen-Totton & Catling 2017). We furthermore explore the effects of the tectonic regime (plate tectonics vs. stagnant lid) and other planet-specific parameters such as the initial mantle temperature, mantle rheology, and oxidation state of the mantle.

For a planet with plate tectonics, carbon recycling into the mantle occurs at subduction zones. The temperature profile of a subduction zone depends on several parameters such as subduction angle and plate speed and varies from one subduction zone to another (Penniston-Dorland et al. 2015). In hot subduction zones, the temperature can exceed the decarbonation temperature, causing the release of CO₂ (Johnston et al. 2011). Most of this CO₂ then finds its way through faults and volcanic units back to the surface and is not recycled into the mantle. In the present-day Earth, most subduction zones are sufficiently cold to avoid decarbonation, which enables carbon recycling into the mantle. We note that Ague & Nicolescu (2014) point out that the release of fluids in subduction zones can still drive decarbonation for present-day subduction zones, which we neglect in this paper, for simplicity. In the early Earth, most subduction zones were hot, decarbonation was efficient, and recycling of carbon into the mantle was difficult (Dasgupta et al. 2004; Dasgupta & Hirschmann 2010).

For a stagnant-lid planet, carbon recycling into the mantle could potentially occur through continuous volcanism, burial and sinking of carbonated crust, eventually followed by delamination. However, the crustal temperature increases with depth and the slow sinking of carbonated crust causes it to be close to thermal equilibrium with its surrounding. Carbonates are stable up to a certain pressure-dependent temperature above which metamorphic decarbonation releases CO₂ that may then find its way back into the atmosphere (Foley & Smye 2018). Whether or not decarbonation occurs therefore depends on the temperatureprofile in the carbonated crust. However, even if decarbonation occurs before the crust reaches the mantle, the slowly sinking crust can store carbon temporarily. After the planet has cooled sufficiently for mantle melting to cease, volcanism and mantle degassing will cease as well; no fresh crustal rock can be produced and weathering will greatly slow down. The atmospheric CO₂ partial pressure is then no longer controlled by interior-surface processes and will tend to remain constant. In the following, we derive a model of the atmospheric CO₂ partial pressure depending on the tectonic regime and thermal state of the planet.

2 Carbon cycling and decarbonation on plate tectonics planets

For simplicity, we restrict our analysis to Earth-sized planets. We only consider planets with an entirely water-covered surface such that continental erosion can be neglected. We neglect the ocean as a carbon reservoir (as in Foley & Smye 2018) however, assuming that the climate is regulated by the long-term carbon cycle. For water-rich planets with an ocean thickness of several hundred kilometres, the climate will rather be controlled by partitioning of carbon between the atmosphere and the ocean (Kite & Ford 2018). Furthermore, we neglect atmospheric escape since Earth-sized planets are not expected to undergo Jeans escape for heavy molecules such as CO₂ (Catling & Zahnle 2009).

The weathering rate for water-covered planets F_w depends on the rate at which fresh crustal volume is produced, , and depending on the scaling used for seafloor weathering, either on the atmospheric CO₂ partial pressure P_{CO_2} or the pore-space temperature T_p. We test both models and scale the weathering rate using the present-day Earth seafloor weathering rate F_{sfw, E}. For the CO₂ -dependent weathering rate , we have (1)

where α ≈ 0.23 is a constant (Brady & Gislason 1997; Foley 2015). The index E denotes present-day Earth variables; variables scaled with their present-day Earth values are noted with an asterisk, that is and . The present-day Earth seafloor weathering rate F_sfw,E is determined by the present-day Earth ingassing rate F_ingas,E, assuming that only a fraction ξ_E of the total weathering rate is due to seafloor weathering, and that only fractions f_E and ϕ_E will not be removed by arc volcanism and decarbonation, respectively, during subduction: (2)

The total flux of carbon into subduction zones F_subd linearly depends on the crustal production rate and on the crustal carbon reservoir relative to the present-day Earth value, , and can be scaled using the present-day Earth seafloor-weathering rate: (3)

We note that Eq. (3) does not include the parameter ξ_E, since it describes the carbon flux into subduction zones of a water-covered planet for which seafloor weathering is the only weathering mechanism. The ingassing rate depends on the total flux of carbon into subduction zones multiplied by the fractions that are not removed by arc volcanism (f) or decarbonation (ϕ): (4)

We assume that the present-day Earth degassing and ingassing rates are in equilibrium and scale the degassing rate with the crustal production rate and with the mantle carbon reservoir. Neglecting minor amounts of carbon in the atmosphere and oceans, the mantle carbon reservoir R_mantle can be described as the difference between the total carbon reservoir R_tot and the crustal carbon reservoir. We have (5)

with R_mantle = R_tot − R_crust. The rate of change of the crustal carbon reservoir can be calculated as the difference between the weathering rate (Eq. (1)) and the rate at which carbon enters subduction zones (Eq. (3)): (6)

The rate of change of the atmospheric carbon is given by , that is (7)

We can now find steady-state values for the crustal and atmospheric reservoirs by setting Eqs. (6) and (7) equal to zero. We obtain (8)

and (9)

with (10)

We note that in steady state, the atmospheric carbon concentration does not depend on the crustal production rate. This isbecause the ingassing rate and the degassing rate both depend on the crustal production rate in the same manner. Equation (9) may be used to analytically calculate the atmospheric carbon concentration on water-covered planets with plate tectonics with variable total carbon reservoirs and carbon fractions f and ϕ that remain stable during subduction. We first restrict ourselves to an Earth-like total carbon reservoir () but also present results varying the total carbon reservoir. Keeping the arc volcanism fraction constant (f^* = 1) throughout our analysis, we calculate ϕ^* dependent on the temperature–depth profiles of subduction zones and a pressure-dependent decarbonation temperature. Upon scaling the model, we use a relative present-day Earth crustal carbon reservoir () of 10%, which is consistent with Sleep & Zahnle (2001).

In order to model a T_p-dependent seafloor-weathering rate , we use an Arrhenius expression following Krissansen-Totton & Catling (2017) (11)

where E_bas is the activation energy for basalt weathering, combined with an empirically derived linear relationship between the pore-space temperature and the surface temperature T_p = a_gradT_s + b_int + 9 with the constants a_grad = 1.02 and b_int = − 16.7 (Krissansen-Totton & Catling 2017). Equations (2)–(5) remain unchanged, while the CO₂ -dependent term in Eqs. (6)–(9) is replaced by the temperature-dependent term .

In both models, the carbon fraction ϕ that remains stable during subduction is a major parameter impacting the atmospheric CO₂ partial pressure and can be calculated using the decarbonation temperature. Foley & Smye (2018) combined a phase diagram for carbonate-bearing oceanic metabasalt (Kerrick & Connolly 2001) with the Holland and Powell DS622 thermodynamic database (Holland & Powell 2011). Foley & Smye (2018) found that the major part of metamorphic decarbonation occurs in a narrow temperature interval by the breakdown of dolomite. We follow their approach and assume that decarbonation occurs at a temperature T_decarb increasing linearly with depth and given by (12)

where z is the depth in m, A = 3.125 ⋅ 10⁻³ K m⁻¹, and B = 835.5 K (Foley & Smye 2018).

The fraction of stable carbon ϕ depends on how the temperature of the carbonated crust increases on its way into the mantle. However, temperature–depth profiles in subduction zones may differ significantly. While in hot subduction zones crustal decarbonation may occur, relatively cold subduction zones allow carbonated crust to reach the mantle. The global carbon fraction ϕ that remains stable during subduction depends on the distribution of the temperature–depth profiles of the actual subduction zones. Since the fraction ϕ may change with time, we derive a parameterization dependent on the thermal state of the planet.

The temperature gradient in subduction zones depends on various parameters, such as the mantle temperature, plate speed, subduction angle, and others. For most subduction zone parameters, it is unclear how they scale with planet evolution. For example, the plate speed increases with the convection strength as derived from boundary layer theory (e.g. Schubert et al. 2001), which would yield a larger plate speed when the Earth was hotter, but at the same time it decreases with the crustal thickness, which has also been larger in the past (Sleep & Windley 1982), possibly causing sluggish plate tectonics (Korenaga 2006). Altogether, the exact dependence of the plate speed on the thermal state of the planet is difficult to predict. In contrast, it is clear that the temperature gradient of a specific subduction zone directly depends on the mantle temperature (England & Wilkins 2004). Therefore, we set the temperature increase along the subducting slab from the surface to a specific depth z′, , to be proportional to the average temperature increase between the surface and the mantle, T_m − T_s: (13)

The constant c represents the specific subduction zone (i.e. hot or cold). We keep the depth z′ constant at 80 km, which is approximately the depth up to which a strong temperature gradient for the top of a subducting slab can be maintained (c.f. Penniston-Dorland et al. 2015 and references therein). We take a present-day distribution of subduction zone temperature gradients between 5 and 10 K km⁻¹. The temperature gradients in some subduction zones may be larger (Penniston-Dorland et al. 2015). However, such temperature gradients cannot be maintained to depths of 80 km. The temperature increase of the hottest subduction zones in Penniston-Dorland et al. (2015) is not larger than 800 K, which corresponds to an average temperature gradient of 10 K km⁻¹ throughout the upper 80 km. Using a present-day surface and upper mantle temperature of 285 and 1650 K, respectively, we obtain c values between 0.3 and 0.6. A small c value represents a rather cold subduction zone compared to the others while a large value represents a hot subduction zone. We assume that c is uniform and constant in time, and therefore the temperature of each subduction zone increases linearly with the mantle temperature. Depending on the thermal state of the planet, we can derive a critical value c_crit (t) that determines whether or not a particular subduction zone allows subducted carbonates to reach the mantle.

The fraction of c < c_crit(t) combines subduction zones that are sufficiently cold to avoid decarbonation while the fraction of c > c_crit(t) combines hot subduction zones in which decarbonation takes place (see Fig. 1). We note that c_crit increases as the mantle cools, such that an increasing fraction of subduction zones avoid decarbonation.

Since the temperature–depth gradient for the top of a subducting slab is particularly large within the first 80 km, we assume that if decarbonation takes place, this occurs before the carbonated crust reaches this depth. If the temperature of a specific subduction zone at 80 km exceeds the decarbonation temperature, decarbonation occurs. On the contrary, if the temperature at 80 km does not exceed the decarbonation temperature, carbon can be recycled into the mantle. The threshold value c_crit, which determines whether or not decarbonation occurs in a specific subduction zone, can thus be obtained by setting the temperature of the top of a subducting slab at 80 km equal to the decarbonation temperature at this depth. Combining Eqs. (12) and (13), we have (14)

where the depth z = z′ = 80 km.

The evolution of the mantle temperature is calculated using a parameterized model of mantle convection. Details of this model can be found in Appendix A. We can now calculate the threshold value c_crit(T_m), which determines whether or not decarbonation occurs in a specific subduction zone. To obtain the global fraction ϕ of carbon that remains stable during subduction, we calculate the fraction of subduction zones for which c < c_crit. Assuming a uniform distribution of c between 0.3 and 0.6 for simplicity, we can calculate ϕ dependent on c_crit as (15)

with 0.3 ≤ c_crit ≤ 0.6. In order to remain consistent with our scaling to the present-day Earth, we also follow Eqs. 14 and 15 when calculating the present-day value ϕ_E of the stable global carbon fraction.

In Sect. 4.1, we first show examples obtained using fixed surface temperatures of 300 and 400 K, respectively. However, it is well known that CO₂ is a greenhouse gas, which exerts a feedback on the surface temperature and on the CO₂ partial pressure. To close the feedback cycle, we use a simple parameterized climate model based on Walker et al. (1981) (see Appendix C).

Fig. 1 Schematic cartoon depicting the distribution of the temperature in different subduction zones at depth z′ with respectto the decarbonation temperature in the early and late evolution.

3 Carbon cycling and decarbonation on stagnant-lid planets

On planets with plate tectonics, CO₂ degassing from the mantle mainly occurs via melting at mid-ocean ridges and is therefore independent of the temperaturegradient affecting carbon subduction, which is one of the reasons reasons why carbon can be recycled into the mantle despite volcanic activity. In contrast, degassing on stagnant-lid planets occurs at hot spots. Partial melting, which ultimately causes degassing, occurs at the transition from the lid to the upper mantle. On its way into the mantle, the buried carbonated crust will pass through the zone at which partial melting takes place. Comparing the decarbonation temperature with the solidus temperature of dry peridotite (e.g. Takahashi 1990) reveals that for depths of up to 250 km, the dry solidus temperature exceeds the decarbonation temperature by more than 500 K. As pointed out by Foley & Smye (2018), the main fraction of carbonates becomes unstable throughout this narrow temperature interval, and only a minor part of carbonates is stable at significantly larger temperatures. Although plumes may increase the temperature beneath the stagnant lid locally by the temperature difference across the bottom boundary layer, this increase is much smaller. The initial temperature drop of 200 K assumed in Tosi et al. (2017) decreases rapidly with ongoing evolution. Sinking carbonate in cold downwellings will have a lower temperature than the laterally averaged temperature, but the temperature difference should also be in the range. The presence of water reduces the solidus temperature. However, large amounts of water would be required to substantially reduce the solidus. A solidus reduction of 100 K would already require a bulk water concentration of more than 200 ppm (Katz et al. 2003). Maintaining sufficient amounts of water in the long-term without water recycling to significantly decrease the solidus would require a large initial mantle water concentration.

Altogether, stagnant-lid planets that are sufficiently hot for partial melting to take place should inevitably be too hot to allow carbon to enter the mantle in substantial quantities. In other words, as long as the mantle temperature of a stagnant lid planet is large enough to cause partial melting to take place regionally, metamorphic decarbonation will occur globally. On the other hand, if the mantle temperature is too low for partial melting to occur and if the crust is gravitationally stable, decarbonation will not occur. This is because the crust will only grow in thickness and carbonated crust will only reach the decarbonation depth if fresh crust can be produced by partial melting. It has been argued that gravitational instabilities could recycle the crust into the mantle (O’Rourke & Korenaga 2012; Johnson et al. 2014), which could cause sinking and decarbonation of part of the crust for planets in their late evolution even if partial melting no longer takes place. For simplicity, we neglect gravitational instabilities in the carbonated crust in our model.

The rates of change of atmospheric carbon reservoir R_atm and crustal carbon reservoir R_crust of a stagnant-lid planet can be calculated as (16)

and (17)

where F_w, F_decarb, and F_degas are the weathering, decarbonation, and degassing rates, respectively. If partial melting and decarbonation occur in stagnant-lid planets, Eqs. (16) and (17) imply that a combined steady state of the atmosphericand crustal carbon reservoirs can only be reached if the mantle degassing rate is zero, which would be the case if all of the carbon were stored in the atmosphere and crust. The reason for this is that a complete decarbonation implies zero ingassing, and therefore, the steady state would require zero degassing as well. However, it is unlikely that stagnant-lid planets can reach this steady state. Rather, the combined degassed CO₂ reservoir, which is distributed between the crust and the atmosphere, will increase with time as long as volcanism is active.

Figure 2 qualitatively illustrates carbonation, burial, and sinking of the carbonated crust, and decarbonation. For simplicity, our model assumes that all the carbon is initially stored in the mantle. Therefore, in early evolution (t₁), the CO₂ flux to the atmosphere only consists of mantle degassing. The weathering rate, and therefore the crustal carbon concentration, remain small. Later in the evolution (t₂), when carbonated crust reaches the decarbonation depth, the CO₂ flux to the atmosphere is given by the sum of the degassing rate and the decarbonation rate, the latter in turn is determined by the weathering rate at the time the crust was produced. As a result, the weathering rate and the carbon concentration of freshly produced crust increase. We note that an opposite effect is the increase of the decarbonation depth upon mantle cooling, which causes the decarbonation rate to diminish.

Weathering is a surface-related process and we assume that the growth of the crust is so slow that every finite layer can be carbonated before the next layer is added. Should massive volcanic eruptions produce large volumes of basalt so rapidly that the underlying crust is not carbonated, we may overestimate the weathering rate and the crustal carbon concentration with this approach. In assuming that the entire formed crust can be subject to weathering, we follow an approach that is typically made for planets with plate tectonics (Sleep & Zahnle 2001).

The evolution of the atmospheric CO₂ reservoir is given by Eq. (16). We follow Eqs. (1) and (2) and use , where X_E is the present-day Earth mid-ocean ridge CO₂ concentration in the melt. For the CO₂-dependent weathering scaling we obtain (18)

We use X_E = 125 ppm (e.g. Burley & Katz 2015) and P_{CO_2,E} = 4 × 10⁻⁴ bar. We again test the T_p-dependent weathering model by replacing the CO₂-dependent term in Eq. (18) with the temperature-dependent term .

For simplicity, we calculate the decarbonation depth assuming a linear temperature profile, thereby neglecting effects of spherical geometry and heat sources in the crust. The temperature throughout the lid with a thickness D_l and boundary layer with a thickness δ_m is then given by (19)

Combining Eqs. (12) and (19), the decarbonation depth is (20)

Parameter values for the carbon cycle model are given in Table 1.

To obtain the evolution of the mantle temperature, degassing rate and crustal production rate of a stagnant-lid planet we employ a 1D thermal-evolution model based on a parameterization of heat transport by convection (e.g. Grott et al. 2011; Tosi et al. 2017). Details of these models can be found in Appendix B (see also Tosi et al. 2017). The numerical model solves the energy-conservation equations of the core, mantle, and stagnant lid. Crustal growth and extraction of incompatible elements occurs via partial melting, and the melt fraction is calculated by comparing the temperature at each depth with the solidus and liquidus temperatures. The mantle viscosity is a function of temperature and water concentration. Since we do not model the water cycle in this paper, we keep the mantle water concentration constant. For the reference case, we use a water concentration of 500 ppm, which is the initial mantle water concentration used in Tosi et al. (2017), but also test a smaller water concentration of 125 ppm. We model CO₂ degassing based on a model of redox melting following Grott et al. (2011) calculating the solubility of CO₂ in the melt dependent on the oxidation state of the mantle. For the reference case, we use an oxygen fugacity of the iron-wüstite buffer (IW), but also test a larger oxygen fugacity of IW+1.

Our numerical model tracks each layer of carbonated crust on its way to decarbonation depth individually and does not assume mixing of carbon in the crust. Whenever a certain layer exceeds decarbonation depth, the carbon content of this layer is set to zero and the CO₂ is supplied to the atmospheric reservoir. However, for an increasing decarbonation depth with time, the lowermost carbonated layer can alternately be subject to decarbonation or remain stable. If the atmosphere-crust equilibrium timescale is fast compared to the temporal resolution of the numerical model, this implies an alternating carbon concentration of the newly produced crust. When this crustal layer in turn is subject to decarbonation, the oscillations can amplify and cause numerical instabilities. To avoid these instabilities, we include a delay of the decarbonation flux such that the decarbonated layer gradually releases its CO₂ to the atmosphere rather than instantaneously. For this purpose, we introduce a temporary volatile CO₂ reservoir within a crustal matrix, which is fed by the decarbonation flux and which releases at each time-step a fixed proportion ω of its CO₂ reservoir to the atmosphere. Implementing such a delay is reasonable to first order: The CO₂ flux throughout the crustal matrix depends on continuous formation and closure of pores and cracks in the crust and on 3D inhomogeneties not captured by our parameterized model. Small values of ω imply on the one hand efficient damping of oscillations and a smooth atmospheric CO₂ curve, but on the other hand a delay in the supply of CO₂ to the atmosphere. For our models, we choose a value of ω = 0.02, which yields both a smooth curve and a relatively small delay. After approximately 80 time-steps (or 0.8 Myr using a temporal resolution of 0.01 Myr), 50% of the CO₂ released by decarbonation is already supplied to the atmosphere. In presenting our results, we illustrate this delay.

Fig. 2 Schematic cartoon depicting carbonated crustal burial and sinking for a stagnant-lid planet. The dots qualitatively represent the crustal carbon concentration.

4 Results

4.1 Planets with plate tectonics

As the mantle cools and an increasing fraction of subduction zones avoid decarbonation, the equilibrium atmospheric CO₂ partial pressure decreases. To illustrate this effect, we plot in Fig. 3 the evolution of (a) the mantle temperature and (b) atmospheric CO₂ partial pressure for initial mantle temperatures of 1800 K (blue), 1900 K (yellow), and 2000 K (red) keeping the surface temperature constant at 300 K (solid) and 400 K (dashed). We use a CO₂-dependent weathering model with α = 0.23. Differences due to the different initial mantle temperature vanish at ≈1 Gyr. The fixed surface temperature is a boundary condition to the thermal evolution of the mantle and to the decarbonation depths (compare the solid and dashed lines). Accordingly, the atmospheric CO₂ approaches different states depending on the surface temperature.

Carbon dioxide is a greenhouse gas affecting the surface temperature. The fixed surface temperature in Fig. 3 does not account for this feedback. In Fig. 4, we apply a simple climate model with the surface temperature dependent on the atmospheric CO₂ partial pressure for a planet orbiting a Sun-like star (Appendix C). We illustrate the influence of different mantle temperatures, which we keep fixed (a and d: 1650 K, b and e: 1800 K, c and f: 1950 K), on the atmospheric CO₂ partial pressure using (a–c) a CO₂-dependent weathering model with α = 0.23 and (d–f) a T_p -dependent weathering model. In addition, we show that the atmospheric CO₂ obtained from the steady-state models (Eq. (9)), plotted as solid black lines, are stable fixed points towards which the atmospheric CO₂ partial pressure evolves. The time-dependent models (dashed lines) start with relative crustal carbon reservoirs of 0, 0.2, 0.4, 0.6, 0.8, and 1 of the total carbon budget (the rest being in the mantle). The initial atmospheric CO₂ concentration is zero. However, this only affects the very early evolution (<0.1 Gyr). The time required to reach a steady state regarding different initial crustal carbon reservoirs is much longer. The models that use the same fixed mantle temperature start to converge at approximately 0.1 Gyr. Neglecting the extreme initial conditions, i.e. all carbon in the mantle or all carbon in the crust, the differences of initial conditions after 1 Gyr are small, and differences of initial conditions at the age of the solar system have vanished. We note that the convergence rate is determined by the crustal production rate, which we keep constant. If the carbon cycle operated faster during early evolution, equilibrium would be reached faster accordingly.

In Fig. 5, we use the equilibrium models of CO₂-dependent weathering with α = 0.23 (solid) and α = 0.25 (dashed) and of T_p-dependent weathering (dashed-dotted) and let the mantle and surface temperatures evolve self-consistently. While we keep the solar luminosity constant in Figs. 5a and b, we use a simple stellar-evolution model to calculate increasing solar luminosity (Appendix C) in plotting Figs. 5c and d. It becomes apparent that the scaling used for seafloor weathering largely determines the evolution of the atmospheric CO₂ and thereby the surface temperature.

Using a CO₂-dependent weathering model (solid and dashed lines in Fig. 5), the influence of the solar luminosity on the atmospheric CO₂ partial pressure is small (compare Figs. 5a and c). Here, the time-dependence is mainly controlled by the cooling of the mantle. However, the surface temperature is greatly affected by the solar luminosity. While in the early evolution (<1 Gyr) the surface temperature decreases with mantle cooling, increasing solar luminosity becomes dominant after 1 Gyr and the surface temperatureincreases accordingly.

Using the T_p-dependent weathering model (dashed-dotted lines in Fig. 5), the surface temperature changes very little with time, and is neither strongly affected by the mantle temperature nor by the solar luminosity. However, the equilibrium CO₂ partial pressure strongly decreases as the luminosity increases.

In Fig. 6, we plot the equilibrium surface temperature as a function of the solar luminosity and the mantle temperature for CO₂-dependent weathering using (a) α = 0.23 and (b) α = 0.25. The surface temperature increases with both luminosity and mantle temperature. Therefore, a cooling mantle can partly compensate for the increasing luminosity, thereby limiting the net change of the surface temperature. We note that we use a thermal evolution model that is based on boundary-layer theory, which results in a rapid early mantle cooling. If thermal evolution models with a slower mantle cooling were applied instead, mantle cooling would compensate for increasing solar luminosity over a longer time span accordingly.

As discussed above, the solar luminosity affects the steady-state surface temperature very little if a T_p -dependent weathering model is used. However, the equilibrium surface temperature is still a function of the degassing rate, which in turn is sensitive to the total carbon reservoir. In Fig. 7, we plot the equilibrium surface temperatureas a function of the total carbon reservoir (relative to the Earth value) and of the mantle temperature for T_p -dependent weathering. We note that the model does not account for a limit of the carbon uptake of the crust due to complete crustal carbonation. However, this limit requires a total carbon budget of more than approximately two times the value for Earth (Foley 2015), which is not shown here. As long as this limit is not reached, the surface temperature does not vary significantly for the T_p -dependent weathering model.

Fig. 3 Evolution of the mantle temperature (a) and the atmospheric CO2 partial pressure (b) keeping the surface temperature constant at 300 K (solid lines) and 400 K (dashed lines). We use the CO2 -dependent weathering steady-state model with α = 0.23 and initial mantle temperatures of 2000 K (red), 1900 K (yellow), and 1800 K (blue). The partial pressure decreases with time as more and more subduction zones become sufficiently cold so as to avoid decarbonation.

Fig. 4 Black solid lines: equilibrium atmospheric CO2 partial pressure. Dashed lines: evolution curves of the atmospheric CO2 using initial crustal reservoirs of 0, 0.2, 0.4, 0.6, 0.8, and 1 times the total carbon reservoir. The mantle temperature is kept fixed at 1650K (a and d), 1800 K (b and e), and 1950 K (c and f), respectively. (a−c) CO2 -dependent weathering model with α = 0.23, (d− f) Tp -dependent weathering model.

Fig. 5 Evolution of (a and c) the atmospheric CO2 partial pressure and (b and d) the surface temperature using initial mantle temperatures of 2000 K (red), 1900 K (yellow), and 1800 K (blue). In a and b we keep the solar luminosity constant while in c and d the solar luminosity is assumed to increase according to stellar-evolution models by a factor of 1.4 in 4.5 Gyr. We test equilibrium models of CO2-dependent weathering with α = 0.23 (solid) and α = 0.25 (dashed) and equilibrium models of Tp-dependent weathering (dashed-dotted).

Fig. 6 Surface temperature as a function of the relative solar luminosity and mantle temperature for CO2 -pressure dependent weathering using α = 0.23 (a) and α = 0.25 (b).

Fig. 7 Surface temperature for temperature-dependent weathering as a function of the total CO2 reservoir relative to the Earth and the mantle temperature. Although the surface temperature increases with both parameters, temperature-dependent weathering is very efficient in keeping the surface temperature relatively low.

4.2 Stagnant-lid planets

In Fig. 8, we show thermal-evolution models of the mantle of a stagnant-lid planet using different initial mantle temperatures of 1800 K (blue), 1900 K (yellow), and 2000 K (red), and the corresponding degassing rates. Neglecting a minor adjustment owing to initial conditions during the first few million years, the degassing rate closely follows the mantle temperature evolution. In particular, an initially hot mantle of 1900 K (yellow) or 2000 (red) results in extensive degassing during the first gigayear. In contrast, an initial mantle temperature of 1800 K results in a slight increase of the mantle temperature and of the degassing rate during the first 0.5 Gyr, followed by a cooling of the mantle and a decreasing degassing rate. After ~1 Gyr, the effect of the initial mantle temperature vanishes. Following that, the degassing rate strongly decreases, and between 4 and 5 Gyrs, volcanism ceases completely.

Modelling the evolution of the crustal and atmospheric CO₂ reservoirs of a stagnant-lid planet (Fig. 9), we assume that all the CO₂ is initially stored in the mantle. We use an initial mantle temperature of 1900 K and compare (a) a CO₂-dependent weathering model with α = 0.23 with (b) a T_p -dependent weathering model. The accumulated degassed CO₂ (magenta) is distributed between three reservoirs: the carbonated crust (red), the atmosphere (blue), and the volatile CO₂ reservoir in the crustal matrix (yellow), which is fed by decarbonation and gradually releases its CO₂ to the atmosphere (implemented to avoid numerical instabilities; see Sect. 3).

Using the CO₂-dependent weathering model (Fig. 9a), the degassed CO₂ in the earlyevolution is mainly stored in the crust such that the atmospheric CO₂ reservoir remains small. With ongoing degassing and an increase of the crustal carbon reservoir, the CO₂ supply rate to the atmosphere increases and with it the atmospheric CO₂ reservoir. After 2–3 Gyrs, the atmospheric CO₂ reservoir decreases despite an increase of the crustal CO₂ reservoir. This is caused by an increase of the decarbonation depth, which results in a reduced decarbonation rate. Since the degassing rate at this time is also small and therefore contributes very little CO₂ to the atmosphere, a net predominance of weathering results in a decrease of the atmospheric CO₂. When volcanism ceases completely between 4 and 5 Gyrs, the atmospheric CO₂ remains approximately constant.

The T_p-dependent weathering model (Fig. 9b) implies an extremely strong weathering feedback and hardly allows any CO₂ to remain in the atmosphere. Initial rapid degassing goes along with an increase of the crustal carbon reservoir and the atmospheric CO₂ remains small. We note that independent of the weathering scaling, the crustal CO₂ reservoir steadily increases as long as partial melting takes place, whereas the atmospheric reservoir decreases inthe late evolution as the decarbonation depth increases. This qualitatively differs from stagnant-lid models that do not account for weathering and decarbonation, in which the atmospheric CO₂ partial pressure steadily increases (e.g. Tosi et al. 2017).

The effect of different initial mantle temperatures becomes apparent from Fig. 10. We again test initial mantle temperatures of 1800 K (blue), 1900 K (yellow), and 2000 K (red), and CO₂-dependent weathering models with α = 0.23 (solid) and α = 0.25 (dashed) and a T_p -dependent weathering model (dashed-dotted). For the CO₂-dependent weathering model, the atmospheric CO₂ is greatly affected by the initial mantle temperature. This is because any CO₂ that is degassed from the mantle cannot be recycled back. A large initial mantle temperature implies rapid CO₂ degassing during early evolution, which in turn causes a rapid increase of the crustal carbon reservoir. Since the crustal reservoir is continuously supplied to the atmosphere via decarbonation, the atmospheric CO₂ partial pressure throughout the entire evolution crucially depends on the initial mantle temperature. Since the decarbonation depth depends on the surface temperature, the difference in the atmospheric CO₂ partial pressure grows larger. Altogether, differences in the initial mantle temperature of 100 K can cause large differences in the surface temperature at 4.5 Gyr of up to 50 K if a CO₂-dependent weathering model is used. In contrast, the T_p-dependent weathering model hardly allows any CO₂ to remain in the atmosphere, independently of the initial mantle temperature. In the late evolution, the remnant CO₂ in the atmosphere is depleted and approaches zero.

Particularly for the CO₂-dependent weathering model the atmospheric CO₂ of a stagnant-lid planet is strongly affected by the accumulated degassed CO₂, which in turn depends on the mantle rheology and the oxidation state (Appendix B), parameters that are poorly known for extrasolar planets. Figure 11 illustrates the sensitivity of the model to these parameters using a CO₂-dependent weathering model with α = 0.23 and initial mantle temperatures of 1800 K (blue), 1900 K (yellow), and 2000 K (red). Figures 11a–c compare the reference case (mantle water concentration of 500 ppm and oxygen fugacity of IW+0, solid) with a case of a reduced mantle water concentration of 125 ppm affecting the viscosity (dashed) and Figs. 11d–f compare the reference case (solid) with a model of a greater oxygen fugacity of IW+1 (dashed).

Although the mantle viscosity impacts the evolution of the mantle temperature (Fig. 11a), the effect on the atmospheric CO₂ partial pressure (Fig. 11b) or surface temperature (Fig. 11c) is less clear. On the one hand, a large mantle viscosity delays mantle cooling, but on the other hand it results in a more sluggish convection, which reduces the rate of melt production. Whether the combined effect results in a higher or lower atmospheric CO₂ partial pressure depends on the initial mantle temperature and is not significant regarding other model uncertainties. In contrast, the oxygen fugacity impacts the evolution of the atmospheric CO₂ partial pressure (Fig. 11e) and the surface temperature (Fig. 11f) even more strongly than the choice of the initial mantle temperature. This is because a greater oxygen fugacity enhances the rate at which CO₂ is degassed without enhancing the melt production rate.

Fig. 8 Panel a: mantle temperature evolution. Panel b: degassing rate for a stagnant-lid planet with an initial mantle temperature of 1800 K (blue), 1900 K (yellow), and 2000 K (red) using a CO2 -dependent weathering model with α = 0.23.

Fig. 9 Evolution of carbon reservoirs for a stagnant-lid planet with an initial mantle temperature of 1900 K and (a) CO2 -dependent weathering model with α = 0.23 and (b) Tp -dependent weathering model. The accumulated degassed CO2 (magenta) is distributed between the carbonated crust (red), the atmosphere (blue), and the volatile CO2 reservoir in the crustal matrix (yellow).

Fig. 10 Evolution of (a) atmospheric CO2 and (b) surface temperature for a stagnant-lid planet with initial mantle temperatures of 2000 K (red), 1900 K (yellow), and 1800 K (blue). We test CO2-dependent weathering models using α = 0.23 (solid) and α = 0.25 (dashed) and a Tp dependent weathering model (dashed-dotted).

Fig. 11 Evolution of mantle temperature, atmospheric CO2, and surface temperature for a stagnant-lid planet with an initial mantle temperature of 1800 K (blue), 1900 K (yellow), and 2000 K (red)using a CO2-dependent weathering model with α = 0.23. The solid lines refer to the reference model (mantle water concentration of 500 ppm and oxygen fugacity of IW+0). Plotting the dashed lines in a–c we test a smaller mantle water concentration (125 ppm) and plotting the dashed lines in d–e we test a larger oxygen fugacity (IW+1).

5 Discussion

We present a model of CO₂ recycling for planets with a water-covered surface with and without plate tectonics. In order to focus on the main mechanisms and feedbacks, various simplifications have been made.

The complexity of individual subduction-zone processes determining their temperature–depth gradients could not be covered within our simple analytical model. In particular, it is difficult to account for parameters that vary from one subduction zone to another, such as the angle and speed of subduction. Instead, we used the present-day distribution of the subduction zone temperature–depth profiles and scaled them according to the evolving mantle temperature. As a result, the fraction of subduction zones avoiding decarbonation increases as the mantle cools. However, this simplification neglects indirect effects of the mantle temperature on subduction zone decarbonation. For example, if the plate speed increases with the mantle heat flux, this may result in limited slab heating during early evolution and therefore in a smaller decarbonation rate than obtained from our model. Altogether, a more detailed study approaching subduction zone temperatures during the thermal evolution of terrestrial planets is required to obtain a more robust connection between the decarbonation rate and the thermal state of the planet.

Modelling decarbonation, we used the scaling for the decarbonation temperature from Foley & Smye (2018), based on Kerrick & Connolly (2001), which was derived for dry carbonates. Therefore, we did not account for CO₂-release by fluid-mediated reactions (Ague & Nicolescu 2014), thereby possibly underestimating the decarbonation flux in subduction zones. In particular, a functional dependence of water release on the temperature–depth profile of a subduction zone, which could then trigger decarbonation, is not captured by our model and may result in a stronger effect of the temperature on the decarbonation rate. Furthermore, we kept the parameter f, describing the fraction of subducted CO₂ that is not released by arc volcanism, constant. Following Eqs. (9) and (10), the equilibrium CO₂ partial pressure directly depends on the stability of carbonates with respect to decarbonation (ϕ) and arc volcanism (f). Future work should consider the effects of temperature and water release on these parameters in order to improve our understanding of the long-term carbon cycle with mantle cooling.

Although we studied planets with a water-covered surface, we restricted our analysis to planets whose habitability remains controlled by the long-term carbon cycle. For planets with a water content of more than 100 times that of the Earth and with a water-layer of several hundred kilometres, the long-term carbon cycle becomes unimportant in stabilizing the climate (Kite & Ford 2018). The habitability of these planets is rather determined by the partitioning of carbon between the ocean and the atmosphere, by the pressure of the water layer affecting volcanism, and by a possible formation of high-pressure ice phases (Noack et al. 2016; Kite & Ford 2018).

A quantitative comparison of the results between the plate-tectonics model and the stagnant-lid model is not straightforward owing to different formulations and parameters of both models. On the one hand, scaling the plate-tectonics model to the present-day Earth is advantageous as it reduces the number of free parameters to a minimum (Appendix A). We assumed an Earth-like CO₂ degassing rate per melt volume for the carbon reservoir of the present-day mantle and scaled it with the latter. Using a thermal evolution model then directly allows us to derive conclusions for the early Earth. On the other hand, the stagnant-lid model requires further assumptions, in particular regarding the rates of crustal production and volcanic degassing (Appendix B), which in turn greatly depends on the oxidation state of the mantle (e.g. Hirschmann & Withers 2008; Grott et al. 2011; Tosi et al. 2017). However, despite the differences in the model formulations, the positive feedback connecting surface temperature, decarbonation, and atmospheric CO₂ affects the habitability of both model planets. The cooling mantle in the late evolution reduces the atmospheric CO₂, regardless of the tectonic mode. For the case of plate tectonics, this is caused by an increasing fraction of subduction zones that avoid decarbonation. For the case of a stagnant-lid planet, this is caused by an increasing decarbonation depth.

The dependence of the seafloor weathering rate on CO₂ partial pressure and/or temperature used in the model is crucial. In contrast to continental weathering, which is known to be a strong function of the surface temperature (Walker et al. 1981), a temperature dependence of seafloor weathering is less certain. Large uncertainties are associated with scaling hydrothermal carbonation with the atmospheric CO₂ concentration, as discussed in Sleep & Zahnle (2001). The effect of using a different value for α becomes apparent from Eq. (9) since the equilibrium partial pressure of CO₂ directly scales with α. Coogan & Dosso (2015) and Krissansen-Totton & Catling (2017) argued that the water temperature at the bottom of the water layer could control the basalt dissolution rate, thereby implying a strong temperature dependence of seafloor weathering. As illustrated in Fig. 5 (for planets with plate tectonics) and Fig. 10 (for stagnant-lid planets), a larger value of α substantially reduces the atmospheric CO₂ by increasing the weathering rate, and a temperature-dependent weathering rate stabilizes the surface temperature even more strongly.

Kadoya & Tajika (2014) argued that continuous CO₂ degassing is a necessary requirement for an Earth-like planet to maintain surface habitability. This is because with continuous erosion and weathering on a planet with plate tectonics and emerged land, the carbon sink is always present and continuous degassing is required to keep the balance. This requirement has been applied to stagnant-lid planets without continental weathering (Foley & Smye 2018). However, on these planets the carbon sink is directly coupled to the carbon source. Stopping mantle degassing will inevitably imply stopping the production of new crust, which will extensively slow down weathering. Whether or not hydrothermal circulation and seafloor weathering can continue in the absence of volcanic activity remains unclear. In any case, a complete carbonation of the upper part of the crust would imply a halt to weathering. If other sinks such as atmospheric escape can be neglected, the change of the atmospheric CO₂ partial pressure will be small. In our stagnant-lid models (Fig. 10), a decline of the rate of volcanism to a minimum (or even a complete stop) does not imply a rapid decrease of the atmospheric CO₂ if a CO₂-dependent weathering model is used (solid and dashed lines in Fig. 10). However, the T_p -dependent weathering model (dashed-dotted lines in Fig. 10) hardly allows any CO₂ to remain in the atmosphere, and therefore the small atmospheric reservoir can be entirely consumed by weathering when the degassing and decarbonation rates become low (at ≈3.5 Gyrs).

Our model is not suited to exploring a transition to a potential snowball state for very small atmospheric CO₂ partial pressures. Using a simple climate model following Walker et al. (1981), we focus on exploring the ability of planets to maintain a sufficiently low atmospheric CO₂ partial pressure in order not to limit their surface habitability by overly high temperatures. However, a complete freeze of the surface of an ocean planet at a low CO₂ partial pressure can also limit the surface habitability. Foley (2019) argued that fluctuations in the weathering rate could cause stagnant-lid planets withvery small total carbon contents to evolve into a snowball state. This latter author showed that in this case, a continuous exchange of CO₂ between the ocean and the atmosphere may be crucial for the planet to recover from this state. Since we found that the initial mantle temperature impacts the atmospheric CO₂ of stagnant-lid planets, in particular if a weak weathering feedback (i.e. a CO₂-dependent weathering model) is used, it could also play a role in the ability of a planet to avoid a permanent snowball state.

The strict separation between the plate tectonics and stagnant-lid model throughout the entire evolution of a planet is certainly a simplification. The early Earth may have been in a stagnant-lid regime (Debaille et al. 2013). Changes in the tectonic mode may be caused by changes of the surface temperature (Gillmann & Tackley 2014), which would establish a feedback to the climate (Lenardic et al. 2016). We did not account for temporal changes in the tectonic mode and applied simple parameterized thermal-evolution models for both tectonic regimes.

We only considered Earth-sized planets with an Earth-like composition in our model. In particular, the fractions f and ϕ that remainstable during subduction are affected by the interior temperature and pressure gradients. If both gradients increase with planetary mass, the effect of varying the planetary mass on f and ϕ will presumably be small. In contrast, if the pressure-dependent viscosity substantially reduces the heat flow for massive Earth-like planets (Stamenkovic et al. 2012), a large fraction of subducting carbonates will remain stable, thereby reducing the steady-state partial pressure of CO₂ in the atmosphere. Furthermore, volcanic degassing may be restricted for massive planets (Kite et al. 2009). For smaller planets with smaller surface gravitational acceleration such as for example Mars, atmospheric escape may provide an additional important sink to CO₂ (e.g. Tian et al. 2009), which is not considered in our model. Certainly, the planetary mass also has an effect on the occurrence of plate tectonics. However, it is not clear whether an increasing planetary mass would make plate tectonics more likely due to a larger shear stress (Valencia et al. 2007) or whether a more buoyant crust would resist subduction (Kite et al. 2009).

6 Conclusions

Since surface erosion on planets entirely covered with water is negligible, volcanism producing fresh uncarbonated crust is the only mechanism to maintain a carbon sink. Therefore, the weathering rate linearly depends on the rate of volcanism, and is therefore directly coupled to the carbon source. This circumstance becomes particularly important if the stability of carbonated crust with depth is taken into account. Accounting for degassing, seafloor weathering, decarbonation, and thermal evolution we explored the habitability of planets with a water-covered surface in the plate-tectonics and stagnant-lid regimes. Our findings can be summarised as follows.

• Coupling seafloor weathering, decarbonation, and interior evolution strengthens the conclusions of previous studies finding that neither emerged land (Abbot et al. 2012; Foley 2015) nor plate tectonics (Tosi et al. 2017; Foley & Smye 2018; Foley 2019) are required for planets to remain habitable in the long-term. For both tectonic regimes, the CO₂ concentration decreases during late evolution. For planets with plate tectonics this is caused by an increasing fraction of subduction zones that avoid crustal decarbonation and for stagnant-lid planets this is caused by an increasing decarbonation depth with mantle cooling.

• For planets with plate tectonics, the initial mantle temperature is of minor importance. Differences in the initial mantle temperature are negligible after ≈1 Gyr in our model, although this time depends on the parameterization of the convective heat transport. With continuous carbon recycling into the mantle, the atmosphere-crust-mantle system approaches an equilibrium state, which does not depend on the initial conditions. At the age of the solar system, initial conditions are negligible. However, linking the stability of carbonates in subduction zones to the mantle temperature has implications for the atmospheric CO₂ and the surface temperature in the early evolution. Whereas previous studies argue for other greenhouse gases than CO₂ in the atmosphere (e.g. Haqq-Misra et al. 2008) or for a strongly temperature-dependent weathering rate (Krissansen-Totton & Catling 2017; Krissansen-Totton et al. 2018) in order to explain temperate conditions in the early evolution of the Earth, our model results indicate that a hotter mantle may partly compensate for the smaller solar luminosity (see Fig. 6).

• For stagnant-lid planets, the initial mantle temperature has an effect on the evolution of the atmospheric CO₂. This is because a large initial mantle temperature results in a large melt fraction and CO₂ degassing rate during early evolution. Once CO₂ has been degassed, it cannot be recycled back into the mantle. Carbonation of the crust can act as a negative feedback to the atmospheric CO₂, but since the decarbonation rate depends on the crustal carbon concentration, and because the entire crustal carbon reservoir will be supplied to the atmosphere repeatedly, this negative feedback is weaker than for planets with plate tectonics in the long-term. Differences in the initial mantle temperature of 100 K can result in differences in the surface temperature at 4.5 Gyr of up to 50 K if a CO₂-dependent weathering model is used, thereby greatly impacting the habitability of the planet. In contrast, a temperature-dependent weathering model is very effective in keeping the surface temperature low, independent of the initial mantle temperature. Our results should also hold for stagnant-lid planets with an emerged surface: even if erosion on such a planet takes place, the thereby exposed crust would already have been carbonated at the time it was produced, such that the weathering rate still linearly depends on the crustal production rate.

• Since the stability of carbonated crust with depth is temperature-dependent, a positive feedback loop between the surface temperature, decarbonation rate, and atmospheric CO₂ is established. Future work studying the difference in the atmospheric CO₂ and the runaway greenhouse for planets with different received solar heat fluxes (c.f. Driscoll & Bercovici 2013; Foley & Driscoll 2016) should account for this positive feedback, as it reinforces differences in the atmospheric CO₂ partial pressure.

Acknowledgements

We thank Brad Foley for his constructive review, which greatly helped to improve this paper. D.H. has been supported through the NWA StartImpuls and N.T. has been supported by the Helmholtz Association (project VH−NG−1017) and by the German research foundation (DFG) through the priority programs 1922 “Exploring the diversity of extrasolar planets” (TO 704/3-1) and 1883 “Building a habitable Earth” (TO 704/2-1).

Appendix A Thermal evolution model for planets with plate tectonics

Calculating the thermal evolution of a planet with plate tectonics, we follow boundary layer theory and a parameterization of the heat transport by convection (Schubert et al. 2001; Davies 2007; Höning & Spohn 2016). The evolution of the mantle temperature is calculated by integrating the energy conservation equation (A.1)

where ρ_m is the mantle density, c_m the heat capacity, V _m the volume, A_m the surface area, and A_c the core surface area. The mantle heat production rate is (A.2)

where t represents the time in gigayears, t_E = 4.5 Gyr is the present-day age of the Earth, Q₀ is the total present-day heat production rate, are the relative present-day heat production rates, and h₁–h₄ the halt-life times of ²³⁸U, ²³⁵ U, ²³² Th, and ⁴⁰ K as given in Korenaga (2008).

The heat fluxes q_u and q_b through the upper and lower thermal boundary layers are given by (A.3)

and (A.4)

where k is the thermal conductivity, δ_u and δ_b are the upper and lower boundary layer thicknesses, and T_s, T_c, and T_b are the surface, core, and bottom mantle temperatures. We use a local instability criterion following Stevenson et al. (1983) and set the Rayleigh number Ra at the upper boundary layer equal to the critical Rayleigh number Ra_crit to calculate the thicknesses of the boundary layer (A.5)

where g is the gravitational acceleration,η the viscosity, and the thermal expansion coefficient. For simplicity, we assume a symmetric temperature profile throughout the mantle and set the thickness of the lower thermal boundary layer equal to the thickness of the upper boundary layer, that is δ_b = δ_u. The viscosity at the upper thermal boundary layer is temperature-dependent: (A.6)

where η_E and T_m,E are the present-day Earth mantle viscosity and mantle temperature, E^* the activation energy, and R the gas constant. The adiabatic temperature profile in the mantle is calculated using the linearized form (A.7)

where d_m is the mantle thickness. The evolution of the core temperature T_c is calculated by (A.8)

where A_c, V _c, ρ_c, and c_c are the surface area, volume, density, and specific heat capacity of the core. Parameter values are given in Table A.1.

Appendix B Thermal evolution model for stagnant-lid planets

We use a parameterized model of mantle convection in the stagnant lid regime and CO₂ degassing that has been explained in detail in Tosi et al. (2017). Therefore, we only give a brief overview of the model here. The conservationof energy is given by (B.1)

where St is the Stefan number controlling the release and consumption of latent heat upon melting and solidification, q_l, A_l, and V _l are the heat flux through the lid and its area and volume, ρ_cr, c_cr, and D_cr are the crustal density, heat capacity, and thickness, and L is the latent heat of melting. The evolution of the stagnant lid is given by (B.2)

where D_l is the lid thicknessand is the radial temperature gradient at the base of the lid (i.e. at the radius r = R_l), which is calculated using the steady-state heat conduction equation (B.3)

where r is the radial coordinate, and Q_l and k are the heat production rate and the thermal conductivity of the lid. The heat fluxes q_u and q_l are obtained following boundary layer theory (Appendix A). We account for the water-dependence of the mantle viscosity as explained in Tosi et al. (2017) (B.4)

where A′ is the pre-exponential factor and the water concentration in parts per million (for the reference models, we use = 500 ppm), E^* is the activation energy, V ^* is the activation volume, and p_m is the lithostatic pressure at the upper mantle. We follow Grasset & Parmentier (1998) and calculate the lid temperature by (B.5)

where Θ = 2.9 accounts forspherical geometry.

Calculating the local melt fraction Φ(r), we compare the mantle temperature profile T(r) with the solidus T_sol(r) and liquidusT_liq(r) temperatures as given in Katz et al. (2003): (B.6)

The volume-averaged melt fraction is calculated as (B.7)

where V _Φ is the melt volume. The crustal growth rate is calculated as (B.8)

where R_p is the planet radius, f_p is a constant describing the fraction of the surface covered by hot plumes in which partial melting takes place, and u is the convective velocity, (B.9)

where u₀ is the characteristic mantle velocity scale and β = 1∕3 is a constant. If the crustal thickness exceeds the lid thickness, we set D_cr = D_l, causing crustal recycling into the mantle.

We account for partitioning of incompatible elements between the crust and the mantle. The concentration of a trace element in the liquid phase X_liq relative to its bulk mantle concentration X_m is calculated from (B.10)

where δ = 10⁻³ is the partition coefficient for heat producing elements (Blundy & Wood 2003). Combining Eqs. (B.6) with (B.10) yields the average concentration of incompatible elements in the crust, and the total rate at which incompatible elements are extracted is set proportional to the crustal production rate (Eq. (B.8)). The volumetric heating rate of the mantle as needed in Eq. (B.1) is then calculated dependent on the respective abundance of heat producing elements in the mantle, and the heat production in the lid (Eq. (B.3)) is derived from the heating rates in the crust and mantle using their relative volume fractions. A more thorough explanation of the thermal evolution model is given in Tosi et al. (2017).

To calculate the CO₂ degassing rate, we follow the redox melting and degassing model of Grott et al. (2011) and assume that the mantle oxygen fugacity is sufficiently low for carbon to be available in its reduced form. Upon partial melting, some graphite is dissolved in the melt and subsequently degassed. The concentration of CO₂ in the melt depends on the carbonate concentration, which in turn is a function of the oxygen fugacity as given in Grott et al. (2011) and Tosi et al. (2017), using scaling constants appropriate for Hawaiian basalts (Holloway 1998). We keep the oxygen fugacity fixed at the iron-wüstite buffer IW but also test an oxygen fugacity of IW+1. The CO₂ degassing rate directly follows from the concentration of CO₂ extracted from the mantle and the melt fraction as given in Eq. (B.7), using a fixed extrusive-to-intrusive volcanism ratio of 2.5 (Grott et al. 2011). More details on the CO₂ degassing model used for stagnant-lid planets can be found in Tosi et al. (2017).

Appendix C Climate model

To test the feedback cycle between mantle temperature, atmospheric CO₂, and surface temperature, we apply a simple climate model following Walker et al. (1981). We relate the atmospheric CO₂ partial pressure to the surface temperature as (C.1)

where T_s,E is the present-day Earth surface temperature, P_{CO_2} the atmospheric CO₂ partial pressure (and P_{CO_2,E} the present-day Earth value), and the difference of the effective temperature is given by (C.2)

where σ is the Stefan-Boltzmann constant, A the albedo, and the solar flux following Gough (1981) is given by: (C.3)

where t is the time in Gyr (t_E = 4.5 Gyr is the age of the present-day solar system) and S_E the present-day Earth solar flux. Parameter values are given in Table C.1.

References

All Tables

All Figures

	Fig. 1 Schematic cartoon depicting the distribution of the temperature in different subduction zones at depth z′ with respectto the decarbonation temperature in the early and late evolution.
In the text

	Fig. 2 Schematic cartoon depicting carbonated crustal burial and sinking for a stagnant-lid planet. The dots qualitatively represent the crustal carbon concentration.
In the text

Fig. 3 Evolution of the mantle temperature (a) and the atmospheric CO2 partial pressure (b) keeping the surface temperature constant at 300 K (solid lines) and 400 K (dashed lines). We use the CO2 -dependent weathering steady-state model with α = 0.23 and initial mantle temperatures of 2000 K (red), 1900 K (yellow), and 1800 K (blue). The partial pressure decreases with time as more and more subduction zones become sufficiently cold so as to avoid decarbonation.

In the text

Fig. 4 Black solid lines: equilibrium atmospheric CO2 partial pressure. Dashed lines: evolution curves of the atmospheric CO2 using initial crustal reservoirs of 0, 0.2, 0.4, 0.6, 0.8, and 1 times the total carbon reservoir. The mantle temperature is kept fixed at 1650K (a and d), 1800 K (b and e), and 1950 K (c and f), respectively. (a−c) CO2 -dependent weathering model with α = 0.23, (d− f) Tp -dependent weathering model.

In the text

Fig. 5 Evolution of (a and c) the atmospheric CO2 partial pressure and (b and d) the surface temperature using initial mantle temperatures of 2000 K (red), 1900 K (yellow), and 1800 K (blue). In a and b we keep the solar luminosity constant while in c and d the solar luminosity is assumed to increase according to stellar-evolution models by a factor of 1.4 in 4.5 Gyr. We test equilibrium models of CO2-dependent weathering with α = 0.23 (solid) and α = 0.25 (dashed) and equilibrium models of Tp-dependent weathering (dashed-dotted).

In the text

	Fig. 6 Surface temperature as a function of the relative solar luminosity and mantle temperature for CO2 -pressure dependent weathering using α = 0.23 (a) and α = 0.25 (b).
In the text

	Fig. 7 Surface temperature for temperature-dependent weathering as a function of the total CO2 reservoir relative to the Earth and the mantle temperature. Although the surface temperature increases with both parameters, temperature-dependent weathering is very efficient in keeping the surface temperature relatively low.
In the text

	Fig. 8 Panel a: mantle temperature evolution. Panel b: degassing rate for a stagnant-lid planet with an initial mantle temperature of 1800 K (blue), 1900 K (yellow), and 2000 K (red) using a CO2 -dependent weathering model with α = 0.23.
In the text

	Fig. 9 Evolution of carbon reservoirs for a stagnant-lid planet with an initial mantle temperature of 1900 K and (a) CO2 -dependent weathering model with α = 0.23 and (b) Tp -dependent weathering model. The accumulated degassed CO2 (magenta) is distributed between the carbonated crust (red), the atmosphere (blue), and the volatile CO2 reservoir in the crustal matrix (yellow).
In the text

	Fig. 10 Evolution of (a) atmospheric CO2 and (b) surface temperature for a stagnant-lid planet with initial mantle temperatures of 2000 K (red), 1900 K (yellow), and 1800 K (blue). We test CO2-dependent weathering models using α = 0.23 (solid) and α = 0.25 (dashed) and a Tp dependent weathering model (dashed-dotted).
In the text

Fig. 11 Evolution of mantle temperature, atmospheric CO2, and surface temperature for a stagnant-lid planet with an initial mantle temperature of 1800 K (blue), 1900 K (yellow), and 2000 K (red)using a CO2-dependent weathering model with α = 0.23. The solid lines refer to the reference model (mantle water concentration of 500 ppm and oxygen fugacity of IW+0). Plotting the dashed lines in a–c we test a smaller mantle water concentration (125 ppm) and plotting the dashed lines in d–e we test a larger oxygen fugacity (IW+1).

In the text