Carbon cycling and habitability of massive Earth-like exoplanets

As the number of detected rocky extrasolar planets increases, the question of whether their surfaces could be habitable is becoming more pertinent. On Earth, the long-term carbonate silicate cycle is able to regulate surface temperatures over timescales larger than one million years. Elevated temperatures enhance weathering, removing CO$_2$ from the atmosphere, which is subducted into the mantle. At mid-ocean ridges, CO$_2$ is supplied to the atmosphere from the interior. The carbon degassing flux is controlled by the melting depth beneath mid-ocean ridges and the spreading rate, influenced by the pressure- and temperature-dependent mantle viscosity. The influences of temperature and pressure on mantle degassing become increasingly important for more massive planets. Here, we couple a thermal evolution model of Earth-like planets of different masses with a model of the long-term carbon cycle and assess their surface temperature evolution. We find that the spreading rate at 4.5 Gyr increases with planetary mass up to 3 $M_\oplus$, since the temperature-dependence of viscosity dominates over its pressure-dependency. For higher mass planets, pressure-dependence dominates and the plates slow down. In addition, the effective melting depth at 4.5 Gyr as a function of planetary mass has its maximum at 3 $M_\oplus$. Altogether, at 4.5 Gyr, the degassing rate and therefore surface temperature have their maximum at 3 $M_\oplus$. This work emphasizes that both age and mass should be considered when predicting the habitability of exoplanets. Despite these effects, the long-term carbon cycle remains an effective mechanism that regulates the surface temperature of massive Earth-like planets.


Introduction
Exoplanets are being discovered at a rapid rate. NASA's Kepler space telescope identified more than 4000 new planet candidates (Fulton et al., 2017), and advances in detection techniques are yielding observations of a vast array of distinct worlds (Madhusudhan et al., 2016). A significant fraction of currently known exoplanets has radii in between those of Earth and Neptune, R P = 1.0-3.9 R ⊕ (Batalha et al., 2013). The density of planets in this group with a radius smaller than 1.6 R ⊕ is generally consistent with a rocky composition, as opposed to the larger planets having lower densities, seemingly suggesting a low-density envelope (Fulton et al., 2017;Lozovsky et al., 2018). Although these small planets are strikingly common around Sun-like stars (Marcy et al., 2014), they have also been identified orbiting low mass stars. In the near future, the only properties of exoplanets that will be able to be explored beyond their ages, masses and radii are their atmospheres. For a few planets detailed information on their atmospheric compositions has been gathered and this will be extended to more super-Earths in the habitable zone in the near future by spectroscopic studies done by the ELT and the JWST .
The prospect of finding conditions at which liquid water is stable at a planetary surface is exciting. Conditions that allow for liquid water would need to be present for extended periods of time for life to develop and to be sustained. A critical factor in the regulation of the surface climate is the long-term carbon cycle (Franck et al., 2000;Kadoya & Tajika, 2014, 2015Haqq-Misra et al., 2016;B. J. W. Mills et al., 2019;Isson et al., 2020).
Volcanic systems have vented CO 2 from Earth's mantle into the atmosphere for billions of years (Kasting & Catling, 2003). Since silicate weathering and thereby the removal rate of CO 2 from the atmosphere depends on atmospheric CO 2 and surface temperature, a negative feedback is established. While these reactions have first been described by Urey (1952), the importance of the long-term carbon cycle as a feedback mechanism regulating the climate evolution of Earth and other planets has been described later (Walker et al., 1981;Kasting, 1989).
Carbon reaches the atmosphere from the deeper parts of the Earth through volcanism.
At mid-ocean ridges, lithospheric plates diverge and mantle rock ascends to replace it. In the process, decompression melting takes place and carbon is degassed. The rate of carbon degassing is of fundamental importance as it influences the width of the habitable zone both for stagnant-lid (Noack et al., 2017) and plate tectonics planets (Kadoya & Tajika, 2014).
At the outer boundary, inefficient carbon degassing could lower the surface temperature to below the freezing point of water, while at the inner boundary high degassing rates could enhance surface temperatures and cause water evaporation (Noack et al., 2017).
Degassing rates depend on the convective regime. While in the Solar System Earth is the only planet with active plate tectonics, the dominating tectonic regime on super-Earths is matter of debate (Kite et al., 2009). Noack et al. (2017) and Dorn et al. (2018) explore degassing rates for super-Earths in a stagnant-lid regime and find that degassing is most efficient on 2-3 M ⊕ planets (M ⊕ is the mass of Earth), depending on the initial mantle temperature of these planets. In this study, we explore the degassing flux and surface temperature for super-Earths with a mass between M ⊕ = 1 and 10 with active plate tectonics. To this end, we couple a parameterized model of mantle convection to a model of the long-term carbon cycle. In contrast to previous carbon-cycle models, we model the degassing flux by including temperature-and pressure-dependent mantle viscosity and by explicitly calculating the melting depth beneath mid-ocean ridges.

Model Set-Up
Our model connects the long-term carbon cycle to the thermal evolution of plate tectonics planets of different masses. Carbon is degassed from the mantle reservoir R man into the atmospheric reservoir R atm by the degassing flux F deg . Exposed weatherable rock reacts with atmospheric CO 2 to form bicarbonate and calcium and magnesium ions, which are washed into the oceans where they convert to carbonate minerals. An increase of CO 2 in the atmosphere speeds up silicate weathering, composed of the continental weathering flux F w and the seafloor weathering flux F sf w , both of which draw CO 2 from the atmosphere.
Carbon is added to the plate reservoir R p and subsequently transported by the subduction flux F sub . A fraction of the subducted carbon is directly degassed back into the atmosphere at arc volcanoes (F arc ), while the rest enters the mantle, completing the long-term carbon cycle.
The basic carbon cycle model outlined above has been discussed extensively in the literature (Walker et al., 1981;Sleep & Zahnle, 2001;Kasting & Catling, 2003;Foley, 2015) and the main equations are given in Appendix A. Here, we focus on the degassing rate and its evolution through time for different planetary masses. Degassing at mid-ocean ridges is the dominant pathway for carbon to reach the surface from the mantle (Orcutt et al., 2019) and is therefore crucial for the long-term carbon cycle. In previous studies, the degassing flux has often been taken as a predefined input flux when considering Earth's evolution (Krissansen-Totton & Catling, 2017;Kadoya et al., 2020) or taken to be controlled by the mantle carbon reservoir when considering the evolution of Earth-sized planets (Sleep & Zahnle, 2001;Höning et al., 2019). Foley (2015) introduced a dependence of the plate velocity -and thereby degassing rate -on the surface temperature, but found that the effect on the degassing rate for planets in the habitable zone is small. Oosterloo et al. (2021) calculated the plate velocity as a function of the mantle heat flow but kept the melting depth constant and did not account for the pressure-dependence of mantle viscosity, making an application to more massive planets challenging. Kadoya & Tajika (2015) also neglected the pressure-dependence of viscosity when considering the carbon cycle for exoplanets.
Following Sleep & Zahnle (2001) and Foley (2015), the degassing flux F deg is given by where f d is the fraction of upwelling mantle that degasses, v p the plate velocity, L length of ridges and d melt the depth where melting begins. The carbon density of a planet is given by the ratio between the mantle carbon reservoir R man and the mantle volume V m . In the following, we will focus on the parameters v p and d melt , which differ substantially as a function of planet mass and thereby control the surface habitability.
With increasing planetary mass, the pressure increases, leading to an increase of the mantle viscosity. According to Tackley et al. (2013), the core-mantle boundary (CMB) region on a 10 M ⊕ planet with an Earth-like composition experiences a ten times pressure increase. This significant pressure increase substantially slows down mantle convection (Miyagoshi et al., 2013). Van den Berg et al. (2019) show that the effect of increasing mass on viscosity could result in a viscosity of two orders of magnitude higher for a 8 M ⊕ than for an Earth-sized planet. This increase in viscosity for higher mass planets is further supported by Schaefer & Sasselov (2015), who show that a M = 5 planet develops a mantle viscosity up to two orders of magnitude higher when a pressure-dependent viscosity case is considered as opposed to a pressure-independent viscosity, underscoring the importance of including the pressure effect on viscosity. However, higher-mass planets do not solely imply an increase in internal pressures, but also exhibit higher internal temperatures. For example, retention of primordial heat could be higher due to a relatively smaller ratio of surface area to planetary volume or due to ineffective convective heat transport caused by the higher pressure affecting viscosity.
We base out interior thermal evolution model on Schaefer & Sasselov (2015) and derive the plate velocity in Eq. 2 from boundary layer theory (Schubert et al., 2001) where κ is the mantle thermal diffusivity, R p and R c are the radius of the planet and core respectively, and δ u is the upper boundary layer thickness, given by where D is the mantle thickness, η( T m , P mm ) is the effective mantle viscosity with the midmantle pressure P mm and the spherically averaged mantle temperature T m , α the thermal expansivity, ∆T is the super-adiabatic temperature increase throughout the mantle, and Ra and Ra crit are the Rayleigh and critical Rayleigh number, respectively.
Since in contrast to Schaefer & Sasselov (2015) we do not model the deep water cycle, we calculate the viscosity for dry silicate material following Stamenković et al. (2011Stamenković et al. ( , 2012.
The effective viscosity is then given by where the values for the reference temperature and pressure are taken to be 1600 K and 0 Pa, respectively, and the pressure-dependent activation volume is given by the scaling law derived by Stamenković et al. (2011) V a = 1.38 + 2.15 · exp −0.065 P mm 10 9 + 10 0.485 · 10 −6 .
The average temperature follows from an adiabatic temperature profile extrapolated to the surface (Schaefer & Sasselov, 2015) with the adiabatic temperature profile where T p is the potential temperature of the mantle.
We calculate the thickness of the lower thermal boundary layer similar to the upper thermal boundary layer (Eq. 3) but use a local instability criterion, since the high pressure at the core-mantle boundary substantially influences the heat flow here. This calculation differs from the calculation of the upper boundary layer thickness (Eq. 3), which depends on the convection rate of the whole mantle, and therefore requires the use of an effective viscosity.
The thickness of the lower thermal boundary layer δ l can then be calculated as where T c and P c are the temperature and pressure at the core-mantle boundary and T l is the temperature at the top of the lower thermal boundary layer.
The initial core temperature T c,i is estimated using a scaling law derived by Stixrude where T c,i⊕ = 4180 K (Fiquet et al., 2010). Following Schaefer & Sasselov (2015), we assume an initial mantle potential temperature of 2520 K for all planetary masses in our reference model, which translates into an initial average mantle temperature of 3000 K for the Earth, increasing with planetary mass. The effect of assuming different initial mantle potential temperatures is assessed separately (Section 3.3). Calculating the heat fluxes out of the core and mantle, we apply a standard parameterized thermal evolution model from Schubert et al. (2001), which is described in Appendix B.
In addition to the plate velocity, the mantle melting depth must be calculated in order to determine the degassing flux. Since we consider melting at mid-ocean ridges, we follow Kite et al. (2009) and assume that melting occurs as long as the mantle potential temperature is higher than the zero-pressure solidus. The melting depth is considered to be the region between the surface and the point where the solidus crosses the adiabatic temperature gradient (T sol > T grad ), as seen in figure 1. Since with increasing planetary mass planets are expected to cool more slowly or even heat up, the melting depth increases with time relative to 1 Earth mass planets (Papuc & Davies, 2008).
For pressures higher than 12 GPa, the solidus temperature becomes irrelevant, since any produced melt is likely unable to rise to the surface: The P-T location of crystal-liquid density inversions is determined by the compressibility, phase equilibria and element partitioning (Agee, 1998). At high pressures, density inversions are most likely if the liquid has a high compressibility relative to the crystalline phase, when the crystalline phase is stable over a large pressure range or when the liquid is 'chemically dense' (Agee, 1998). According to Agee (1998), crystal-liquid density inversions occur in the Earth in the transition zone and lower mantle, due to the extensive stability range of phases such as garnet and perovskite. For simplicity, we follow Noack et al. (2017) and use a fixed value of 12 GPa as estimated by Agee (1998) and Ohtani et al. (1995). In our model, we impose a limit to the effective melting depth at this pressure.
As noted above, the remaining carbon fluxes are calculated following standard carbon cycle models (Sleep & Zahnle, 2001;Foley, 2015) and given in Appendix A. All parameters that directly depend on the size and/or mass of the planet (such as subduction zone length, land area, and gravity) are adjusted accordingly. For Earth-sized planets, the initial carbon distribution between the atmosphere, crust, and mantle has been shown to become unimportant after 1 Gyr (Foley, 2015;Höning et al., 2019). For simplicity, we assume for our reference model that all carbon is initially stored in the mantle, but will elaborate on the effect of the initial carbon distribution in Section 3.3. We also do not explicitly model ocean chemistry, which only affects the climate on considerably shorter timescales, since the carbonate-silicate cycle takes over control of the system on timescales longer than ≈1 Myr (Sundquist, 1991;Colbourn et al., 2015;Höning, 2020). Valencia et al. (2006Valencia et al. ( , 2007 calculate the interior structure of super-Earths up to 10 Earth masses M ⊕ . The authors find that the dependence of several parameters on mass can be adequately described by a power-law relationship of the form: where Z is the considered parameter and ξ is a scaling exponent. In this way parameters such as gravitational acceleration, mantle density and planetary radius can be related to planet mass. For example, we use the following power law for the mantle density ρ m : Values for ξ for the other scaled parameters can be found in Table 2.

Results
We first focus on the effective melting depth at 4.5 Gyr for planets with masses between 1 and 10 Earth masses (Section 3.1). In Section 3.2, we elaborate on the interior and surface temperature evolution, and Section 3.3 addresses the effect of different initial mantle temperatures and carbon concentrations.

Effective Melting Depth and Plate Thickness at 4.5 Gyr
The melting depth is a main parameter controlling the surface temperature in our model as it directly affects the degassing rate. With increasing planetary mass, both the mantle temperature at 4.5 Gyr and the solidus temperature increase, which are competing factors determining the melting depth. In addition, the 12 GPa limit to positively buoyant melt moves to shallower depth. This limit is imposed to account for the densification of magmatic liquids up to the point where the melt will become neutrally or even negatively buoyant with respect to co-existing minerals (Agee, 1998;Ohtani et al., 1995). All three parameters for different planetary masses are illustrated in Fig. 2, where the effective melting depth is marked. It becomes apparent that as long as the melting depth is reached for pressures below the 12 GPa limit, the melting depth at 4.5 Gyr steadily increases with planetary mass. In contrast, for planets whose effective melting depth is determined by the 12 GPa   Gyr as a function of planetary mass. For the effective melting depth, we find a minimum at 3 Earth masses: On the one hand, the mantle temperature generally increases stronger with planetary mass than the solidus temperature does, resulting in an increase of the melting depth. On the other hand, the 12 GPa limit moves shallower to shallower depths, which becomes dominant for planets from 3 Earth masses. As a result, the effective melting depth, i.e. the depth that controls degassing, has its peak at 3 Earth masses (see Fig. 3a).
For the plate thickness (Fig. 3b), we also find a minimum at 3 Earth masses. This is a result of the pressure-dependence of viscosity, which dominates over its temperaturedependence for planets more massive than 3 M ⊕ . Since the plate velocity directly follows from the boundary layer thickness (see Eq. 2), this effect also promotes high degassing rates particularly for planets of 3 Earth masses. We note that this minimum does not exist for a model that neglects the pressure-dependence of viscosity (dashed curve, following Valencia et al. 2007). Altogether, accounting for pressure-dependent viscosity, both the plate thickness and the melting depth promote a peak of the degassing rate at 4.5 Gyr for

Climate Evolution for Different Planetary Masses
Increasing planetary mass generally implies a delayed cooling of the mantle or even an increasing temperature with time due the higher mantle viscosity. As a consequence, more massive planets remain hotter for a longer period of time (Fig. 4a). However, the plate velocity ( Fig. 4c) depends on the mantle viscosity and therefore not only on temperature but also on pressure. From ≈1 to 4.5 Gyr, this combination results in a particularly high plate velocity for planets with 3 Earth masses. This result, however, depends on the initial mantle temperature (discussed later, Section 3.3). Since the degassing rate (Fig. 4d) depends on both the plate velocity and the melting depth, its evolution is similar to that of the plate velocity, with an even stronger decrease of the degassing rate with time for 1 Earth mass planets compared to more massive planets. This is a consequence of the more rapid cooling of 1 Earth mass planets, which strongly reduces the melting depth with time. least, does not significantly increase). We also find that the most massive planets have a surface temperature below the freezing point of water in their early evolution, (discussed in more detail in Section 4).
Considering a planet with only 1% emerged land area (Fig. 5c, d), the differences between planets of different masses increase, since climate regulation does not work as efficiently in the absence of a substantial contribution of temperature-dependent continental weathering. However, the trend found in Fig. 5 (a and b)  neglected, the surface temperature in our model with a degassing rate dependent on plate velocity and melting depth depends on the thermal evolution more strongly. We also find that the high initial degassing rate for 1 Earth mass planets in combination with the limited weathering rate owed to the small land fraction imply hot, potentially uninhabitable surface conditions during the first 0.5 Gyr.

Initial Mantle Temperature and Carbon Budget
For Earth-sized planets the temperature-dependence of viscosity ensures that the late evolution hardly depends on the initial mantle temperature (Schubert et al., 2001). This does not necessarily apply to more massive planets if the pressure-dependence of viscosity is considered (Schaefer & Sasselov, 2015). For massive planets, the degassing rate at 4.5 Gyr and therefore the atmospheric CO 2 abundance increase with the initial mantle temperature, resulting in higher surface temperatures (see Fig. 6a). However, the planet mass that yields the highest surface temperature at 4.5 Gyr is not very sensitive to the initial mantle temperature: In most cases, the peak remains at 3 M ⊕ , even though the initial mantle temperature is varied over a large range. However, for very high initial mantle temperatures, the peak is slightly shifted towards higher planet masses: For example, for a very high initial potential temperature of 2750 K, we find the peak at 4 M ⊕ (Fig. 6a).
Interestingly, numerical models of stagnant-lid planets  show a similar behavior, even though in these models the degassed CO 2 accumulates in the atmosphere over the entire evolution of the planet and is not only determined by the degassing rate at 4.5 Gyr at 4.5 Gyr for planets between 2 and 4 M ⊕ with a shift of the peak with higher initial mantle temperatures towards higher planet masses. We note that the effect of the initial mantle temperature on the surface temperature at 4.5 Gyr is stronger for stagnant-lid than for plate tectonics planets, since the initial mantle temperature particularly controls early degassing, which makes up a significant part of the accumulated degassed CO 2 abundance on stagnant-lid planets in their late evolution.
The chemical composition of rocky exoplanets is thought to vary greatly (Bond et al., 2010;Madhusudhan et al., 2016;Moriarty et al., 2014). Even for Earth it is still debated which combination of meteorites best represents Earth's bulk chemical composition (Moynier et al., 2015). In order to assess the effect of the planetary carbon density on the habitability of super-Earths, we ran models varying the carbon concentration and present the results in Fig. 6b, keeping the land fraction constant at present-day Earth's value. Since the atmospheric CO 2 is controlled by the degassing rate, which in turn depends linearly on the mantle carbon concentration, our results do not change qualitatively. Similar to Fig. 6a, we find a maximum surface temperature at 4.5 Gyr for planets of 3 Earth masses. The effect of planetary mass becomes more pertinent for higher bulk carbon concentrations, which has important implications for assessing the habitability of exoplanets: Kopparapu et al. (2013) argues that around 340 K the moist-greenhouse effect is reached, where the stratosphere becomes water-dominated and hydrogen loss to space occurs. This state can be considered as the inner edge of the habitable zone. For carbon-rich planets near the inner edge of the habitable zone, it is therefore crucial to consider the planetary mass when assessing the potential for life.

Discussion
Coupling a model of the long-term carbon cycle to thermal evolution models with temperature-and pressure-dependent viscosity and to a model of the melting depth, we assessed the influence of planetary mass on the surface temperature-evolution of Earth-like planets. In this section, we compare our model to previous studies, discuss simplifications and limitations of our model, and consequences for the habitability of exoplanets.
Starting with Walker et al. (1981), models of the long-term carbon cycle for Earth and other planets have been used to assess the climate evolution on Earth and other plate tectonics planets (Sleep & Zahnle, 2001;Kasting & Catling, 2003;Foley, 2015). A key parameter that influences the climate on long timescales is the mantle degassing rate, which is balanced by the weathering rate for specific values of atmospheric CO 2 and surface temperature. For Earth's evolution, it is suitable to approximate the degassing rate over time with scaling laws (Krissansen-Totton & Catling, 2017;Krissansen-Totton et al., 2018;Kadoya et al., 2020), but in particular an extrapolation of the model to planets of different masses requires the consideration of the interior evolution. Rushby et al. (2018) modelled the carbon cycle for planets of different masses by setting the degassing rate proportional to the internal heat production rate, thereby neglecting effects of the planet mass on the mantle viscosity and melting depth. As a result, they find only a little influence of planet mass on surface temperature, with a small trend towards higher surface temperature with increasing mass (<3 K difference between 1 and 10 Earth masses). However, this simplification neglects the fact that the melting depth below midocean ridges crucially depends on the pressure-temperature profile and that the mantle heat flow (which controls the plate velocity) is determined by pressure-and temperaturedependent viscosity. In contrast, we find a peak in the surface temperature for 3 Earth mass planets, which is ≈8 K above the surface temperature for a 1 Earth mass planet. At higher masses, the surface temperature at 4.5 Gyr decreases again and may even fall below the surface temperature that is found for a 1 Earth mass planet. Whereas planets of 1 Earth mass have their highest mantle temperature and therefore degassing rate throughout the first 1 Gyr, after which the initial mantle temperature becomes unimportant (Schubert et al., 2001), the mantle temperature of more massive planets is much longer influenced by its initial value. Regardless of that, at 4.5 Gyr we find a peak in the surface temperature that is hardly influenced by initial conditions. Kadoya & Tajika (2015) include a thermal evolution model of the planetary interior and calculate the degassing rate as a function of melting depth and plate velocity for planets up to 5 Earth masses. These authors use a viscosity law that neglects a pressure-dependence, which results (a) in a steadily increasing degassing rate with planetary mass and (b) in a steadily decreasing degassing rate with time for all planetary masses. These results are in accordance with Oosterloo et al. (2021), who calculate the plate velocity from a thermal evolution model of the mantle and derived the surface temperature evolution for different compositions and masses. In contrast, our results indicate that these trends do not hold if a pressure-dependence of viscosity is considered. In addition, we find that including a pressure limit to the melting depth substantially reduces the melting depth for massive planets.
The peak in surface temperature at 4.5 Gyr for planets of ≈ 3 Earth masses is similar to what has been found for stagnant-lid planets using numerical models of mantle convection in combination with a parameterization of the degassing rate (Noack et al., 2017;. The reason for this is similar: First, the pressure-dependence of viscosity inhibits efficient heat flow and thereby high convection rates (or plate velocities in our model) for massive planets, whereas 1 Earth mass planets have already cooled down more strongly at 4.5 Gyr. Second, the increasing mantle temperature with planetary mass generally increases the melting depth up to the imposed pressure-limit to degassing, which reduces the effective melting depth for planets larger than 3 Earth masses.
Whether or not plate tectonics operate on exoplanets is a matter of debate and may depend on the plate thickness (Valencia et al., 2007), plate buoyancy (Kite et al., 2009), and temperature (Foley et al., 2012;Noack & Breuer, 2014). Even for Earth, a transition from a stagnant-lid tectonic regime to modern plate tectonics may have occurred ≈ 3 billion years ago (Gerya, 2014;Naeraa et al., 2012). We note, however, that the surface temperature at 4.5 Gyr is hardly influenced by the nature of the early tectonic regime, as the atmospheric CO 2 is controlled by an equilibrium between degassing and weathering, with a residence time of has a large influence on the occurrence of plate tectonics. Altogether, plate tectonics on planets somewhat more massive than Earth seems to be likely, although further studies are certainly needed.
An additional simplification of our model is the assumption of a constant land area.
However, Abbot et al. (2012) argue that lower land fractions should be expected for the more massive terrestrial planets. A planet's mass increases stronger than its surface area, which results in deeper oceans. The higher gravity furthermore reduces topography and creates shallower ocean basins (Cowan & Abbot, 2014), increasing the surface area covered with water.
Considering stellar evolution (Fig. 5b), the temperatures found in the early evolution for massive planets are slightly below the freezing point of water. However, this does not necessarily imply a snowball state, since these temperatures should be regarded as globally averaged values and do not rule out an equatorial section with temperatures above freezing.
According to Warren et al. (2002), at surface temperatures above 261 K ice layers are too thin (<1 m) to form coherent sheets. This point is supported by Feulner & Kienert (2014) who find that the lowest steady-state configuration for a partially ice-covered state in the Sturnian glaciation generates a global mean surface air temperature of 266 K. For the Marinoan glaciation the coldest temperature at which ice-free ocean regions still exist is 268 K. In fact, Ye et al. (2015) found preserved benthic macroalgae in the Marinoan-age Nantuo Formation in South China. According to the authors, this suggests that during glaciation muddy substrates existed in the photic zones, suggesting open-water areas in coastal environments. These areas might have acted as a shelter for benthic macro algae.
Together with photoautotrophs, macro algae could have been part of coastal food webs and the carbon cycle, just as they are today in glacial environments. In addition, autotrophes that do not rely on solar energy could exist below an icy crust. Altogether, life as we know it on Earth could potentially still evolve under these conditions. Nevertheless, as discussed above, massive planets are expected to keep their internal heat longer than planets of 1 Earth mass (see Fig 4a). As a consequence, the degassing rate of massive planets remains high for a longer period of time, or even increases with time (see Fig 4d). In combination with increasing stellar luminosity, the surface temperature of a massive planet that is in the conservative habitable range in the early evolution would increase stronger with time than that of a 1 Earth mass planet. We note that the initial mantle temperature affects the early degassing rate and could therefore also affect the onset time of habitability, but such a shift in time would not significantly affect the total duration of the habitable period. Altogether, with increasing mass of the planet, the time window for life to evolve on it appears to be become smaller.

Conclusions
In this paper, we have coupled a model of the long-term carbon cycle to a thermal evolution model of the mantle with temperature-and pressure-dependent viscosity and applied it to different planet masses. In particular, we calculated the degassing rate dependent on the melting depth and plate velocity. Our conclusions are summarized below: • The planet mass plays an important role in the evolution of the degassing rate, which causes different points in time throughout their evolution where the atmospheric CO 2 reaches its maximum. While low-mass planets have a high atmospheric CO 2 content particularly in their early evolution, high-mass planets are expected to build up CO 2rich atmospheres later.
• At 4.5 Gyr, surface temperature increases with planetary mass up to ≈3 Earth masses, since (i) the temperature-dependence of mantle viscosity dominates over its pressuredependence, which yields higher plate velocities with increasing planetary mass, and (ii) the melting depth increases with the mantle temperature. However, above ≈3 Earth masses, the surface temperature decreases with planetary mass, since (i) the pressure-dependence of viscosity dominates and (ii) a pressure-limit above which melt is not positively buoyant reduces the effective melting depth.
• Despite these effects of planetary mass, the long-term carbon cycle remains an important stabilizing feedback mechanism for massive Earth-like planets for a wide range of bulk carbon densities.

A1 Seafloor Weathering Flux
Seafloor weathering adds carbon to the oceanic plates, resulting in the withdrawal of CO 2 from the combined atmosphere and ocean reservoir. The CO 2 dissolved in seawater hydrothermally alters ocean floor basalt. Exactly how and to what extent the seafloor weathering flux operates is not fully understood at present. Brady & Gíslason (1997) showed that the CO 2 uptake during hydrothermal alteration directly depends on the partial pressure of CO 2 in the atmosphere, P CO2 . In addition, seafloor weathering demands the supply of fresh rock, which is created at mid-ocean ridges. The plate velocity therefore dictates the supply of weatherable material. These parameters are combined in the following equation as given by Foley (2015): where v ⊕ is the modern day plate speed on Earth. The plate speed v p is derived in equation 2, and as ratio to v ⊕ depicts the effect of spreading rate. This makes the seafloor weathering flux mass dependent, where increasing planetary size results in a higher F sf w . P * CO2 is the partial pressure for Earth's conditions and α = 0.25. The present-day value for the seafloor weathering flux is 1.75 · 10 12 mol yr -1 as determined by B. Mills et al. (2014).
Other scalings of seafloor weathering exist: Krissansen-Totton & Catling (2017) argue for a direct temperature dependence of basalt dissolution and therefore of seafloor weathering.
The effect of the applied scaling on our results is expected to be small, since continental weathering (section A2) is the dominating flux controlling the climate, and directly depends on both, temperature and atmospheric CO 2 .

A2 Continental Weathering Flux
The continental weathering flux F w is calculated following Foley (2015) as follows: with the asterisk indicating values at present, f l denoting the land fraction, E a the activation energy, and R g the universal gas constant. The exponents a and β are constants. A breakdown of the remaining parameters is provided in the following.
A kinetically-limited continental weathering regime is assumed in this model. A kinetically controlled regime, although ill-constrained, is suggested for the Earth (West et al., 2005), indicating that the supply of weatherable bedrock exceeds the physical erosion rates necessary for complete denudation. In this regime, the parameter F ws calculates the supply limit to weathering and is retrieved by Foley (2015) as follows: where f l is the fraction of land above sea level that is subject to surface weathering, A P the surface area of the planet, E r is the erosion rate by physical processes, k cc the fraction of the cations Mg, Ca, K and Na in the continental crust, ρ r the regolith density and m cc the molar mass of the aforementioned cations.
With the supply limit to continental weathering established, three physical parameters remain to be determined for the weathering flux: the partial pressure of atmospheric CO 2 , the temperature evolution of the atmosphere over time, and the saturation vapor pressure.
The partial pressure of atmospheric CO 2 can be related to the amount of carbon present in the atmosphere, R atm . The initial value for R atm is assumed to be zero and CO 2 is added to the atmosphere by the degassing flux. The partial pressure affects the weathering rate by increasing the reaction rate with increasing PCO 2 , and can be calculated by (Foley, 2015): where mCO 2 is the molar mass of CO 2 .
The surface temperature can be derived by using the effective temperature T e . The effective temperature depends on the amount of solar irradiation and the albedo: with S the solar irradiance, A the albedo of the planet, and σ the Stefan Boltzmann constant.
For the model that considers stellar evolution, S is assumed to increase linearly by a factor of 1/3 throughout 4.5 Gyr (Ribas, 2009). The surface temperature follows from the effective temperature (Foley, 2015): where T * , T * e , and P * CO2 are the surface temperature, effective temperature, and CO 2 partial pressure of present-day Earth. Finally, the saturation vapor pressure P sat represents the fluctuation in runoff with changing surface temperature and can be calculated as follows (Foley, 2015) with m w denoting the molar mass of water, L w the latent heat of water, P sat0 the reference saturation vapor pressure and T sat0 the reference saturation vapor temperature.

A3 Subduction Flux
The concluding part of the long-term carbon cycle in this model is the subduction flux.
Unlike other atmospheric gases, CO 2 reacts with water to form carbonic acid rain, dissolving silicate rock through weathering. Carbon is then removed from the atmospheric system by chemical reaction of the dissolved CO 2 with Ca and Mg containing silicate minerals. This can be in the form of sediment as limestone (CaCO 3 ), or as organic matter (CH 2 O), and is emplaced on the oceanic plates. The fate of 99% of the oceanic plates is to be subducted.
Part of the carbon content on the plates will be accreted to the overriding plate, but it is estimated that this term is approximately balanced by the addition of carbon to the subducting plate, through tectonic erosion (Sleep & Zahnle, 2001). At any time, the amount of carbon R p on the oceanic plates can be calculated by the mass balance equation: The weathering flux F w as obtained by equation A2 is divided by two to account for the carbon that is re-released to the atmosphere when the carbonates are formed. Then, the total amount of subducted carbon can be obtained by (Foley, 2015): where A plate represents the area of the ocean plates. The area of the plates is retrieved by multiplying (1 − l f ) with the surface area of a planet. The parameters v p and L are the velocity of the plates and length of trenches, respectively.
It needs to be accounted for that not all of the subducted carbon will reach the deeper parts of the mantle, but instead that part of it will be re-emitted into the atmosphere by arc volcanism. This is achieved by establishing the arc volcanism flux F arc : where f is the fraction that degasses by arc volcanism. As of today, f is not well-constrained and estimates vary as much as 25-70% (Foley, 2015). Hence, an average value of 50% will be assumed since it does not affect the outcomes in significant ways.
Finally, the reservoirs need to be connected through the fluxes. The amount of carbon on oceanic plates was earlier established in equation A8, so the mantle and atmospheric reservoirs are left to be determined. The quantity of carbon in the atmosphere and oceans combined is given by the mass balance equation (Foley, 2015): and the amount of carbon in the mantle at any given time is:

Appendix B Thermal Evolution Model
The spherically averaged mantle temperature T m evolves as where ρ m , C p , V m are the mantle density, heat capacity, and volume, A s and A c are the surface area of the planet and the core, and q s and q c are the surface and core heat flux, respectively. For the heat production Q m (t) we assume Earth's relative mantle abundances of the main radiogenic heat producing elements K, U and Th, as given by (Korenaga, 2008): where t is the time in Gyr, Q 0 is the total present-day heat production rate, and q 1 ...q 4 and h 1 ...h 4 are the relative present-day heat production rates and half-life times of 238 U, 235 U, 232 Th, and 40 K as given in Korenaga (2008).
The mantle heat flux is given by where T u is the temperature at the base of the upper boundary layer (the crust-mantle boundary), T s the surface temperature.
The core temperature evolution is calculated based on the assumption that there are negligible concentrations of K, U and Th in the core, and is given by