© ESO 2021

1 Introduction

In the past two decades, astronomy has seen the emergence of the field of exoplanets, starting with the discovery of the first exoplanet around a solar-type star (Mayor & Queloz 1995). Since the discovery of this hot Jupiter, the number of discovered exoplanets has increased to over 4200, including planets which are Earth-like both in terms of mass and in terms of effective temperature (Schneider et al. 2011; van Hoolst et al. 2019). A key question with respect to these Earth-like exoplanets is whether these planets are habitable and whether they are able to remain habitable over timescales long enough for life to develop and evolve (e.g., Kasting & Catling 2003).

Notable examples of potentially habitable exoplanets are Proxima Centauri b (Anglada-Escudé et al. 2016), Trappist-1 e, f & g (Gillon et al. 2017), and Teegarden’s Star c (Zechmeister et al. 2019). However, these planets reside in tight orbits around red dwarfs, which confronts these exoplanets with environments widely different from Earth. Examples are the frequent stellar flares associated with red dwarfs (Vida et al. 2017; Günther et al. 2020) and the tidal interaction which likely leaves these planets tidally locked to their host star. These effects can compromise the habitability of exoplanets around M-type stars.

Rocky exoplanets are also detected in the habitable zones of solar-type¹ stars. An example is Kepler-452b, an exoplanet of R ≈ 1.6 R_⊕ which could be rocky in composition (Jenkins et al. 2015). Kepler-452b orbits around a solar-type star which is ~ 6 Gyr old, which is about 1.5 Gyr older than the Sun. As a consequence, Kepler-452b receives ~ 10% more flux from its host star than Earth, which could result in the orbit of Kepler-452b to lie near the inner edge of the habitable zone inferred from conservative estimates (Kopparapu et al. 2013; Jenkins et al. 2015). This could mean Kepler-452b has entered a so-called moist-greenhouse state where most of the planetary water budget is lost to space (Kasting et al. 1993; Abbot et al. 2012; Kopparapu et al. 2013). However, as Sun-like stars such as Kepler-452 gradually increase in luminosity over time, Kepler-452bmust have been cooler in the past, which could have allowed for the existence of surface liquid water for most of its lifetime (Jenkins et al. 2015). Altogether, the evolution of Sun-like stars can induce significant increases in the surface temperatures of orbiting planets, which may in turn compromise the long-term habitability of these planets.

A negative feedback mechanism which could offset the effects of brightening host stars is the long-term carbon cycle, which can act as a planetary thermostat (Walker et al. 1981; Kasting et al. 1993). This mechanism can control the evolution of the global climate over geological timescales via the recycling of CO₂ between the atmosphere and ocean, crust and mantle. The climate is regulated via temperature-dependent silicate weathering, which removes CO₂ from the atmosphere and provides a negative feedback mechanism for the surface temperature. If surface temperatures become lower, the weathering rate is suppressed, allowing CO₂ to build up in the atmosphere over time via volcanic outgassing, increasing surface temperatures. Conversely, high surface temperatures give rise to an increased weathering rate, removing more atmospheric CO₂ and therefore lowering surface temperatures. The carbon removed from the atmosphere is deposited on the seafloor as carbonate rocks. Subduction of the seafloor allows for the transport of carbon toward Earth’s mantle, or to degas back into the atmosphere via arc volcanism. Altogether, the long-term carbon cycle is thought to have allowed Earth to maintain temperatures on its surface favorable for the evolution of life, despite the significant increase in solar luminosity since Earth’s formation (Sagan & Mullen 1972; Kasting et al. 1993).

Due to the importance of the carbon cycle for Earth’s long-term habitability, the efficiency at which such a feedback mechanism operates is also relevant to consider for the habitability of Earth-like exoplanets. Here, efficiency refers to the possibility of carbon to cycle between the mantle and surface, and the timescale on which this process controls and regulates the atmospheric CO₂ pressure. Observational constraints on the characteristics of these exoplanets remain limited to mass-radius measurements. This makes the planetary characterization required to answer these questions challenging. However, these measurements do reveal a large diversity in the mass and bulk composition of Earth-like exoplanets. For example, Kepler-452b has been shown to likely have an interior that is composed of a larger faction of rock than Earth’s (Jenkins et al. 2015). These differences in interior composition affect the thermal and chemical structure of the interior, and the composition and amount of volatile material outgassed into the atmosphere (Noack et al. 2014; Dorn et al. 2018; Ortenzi et al. 2020; Spaargaren et al. 2020). The thermal budget of a planet is affected by the rate of radiogenic heating, which predominantly occurs via the decay of the long-lived radioactive refractory elements thorium and uranium (Schubert et al. 2001). A first-order estimate for the abundances of the radioactive isotopes in the mantles of exoplanets are the abundances found in the atmosphere of the host star (Unterborn et al. 2015). Up to now, the implications of this diversity in interior structure for the efficiency of long-term carbon cycling and feasibility as a climate-regulation mechanism remain poorly understood.

Crucial for a carbon cycle similar to Earth’s is the presence of plate tectonics, which facilitates the transport of carbon into the mantle via subduction. In addition, plate tectonics replenishes atmospheric CO₂ via volcanic degassing at mid-oceanic ridges and subduction zones (Kasting & Catling 2003). The plate speed and prevalence of plate tectonics on Earth-like exoplanets, however, remain poorly understood (Berovici et al. 2015). While some modeling studies suggest plate tectonics is generally likely on Earth-like exoplanets (Valencia et al. 2007; Foley et al. 2012; Tackley et al. 2013), other studies disagree (O’Neill et al. 2007; Kite et al. 2009; Stamenkovic et al. 2012; Noack & Breuer 2014). The uncertainty in the requirements for plate tectonics to emerge is typically circumvented in planetary evolution models by imposing the presence (e.g., Foley 2015) or absence (e.g., Tosi et al. 2017; Foley & Smye 2018; Foley 2019) of plate tectonics.

Even when assuming the presence of plate tectonics, the diversity in the interiors of Earth-like exoplanets leaves a large parameter space to explore. This originates from the fact that plate tectonics is the surface manifestation of convection in Earth’s mantle (e.g., Schubert et al. 2001), which in turn depends on the thermal and chemical structure of the mantle (e.g., Spaargaren et al. 2020). In addition, the manner in which plate tectonics couples to these interior properties remains poorly understood. Even for Earth, the history of plate tectonics is not fully characterized (Palin et al. 2020). Usually studies considering planets in the plate tectonics regime explore the influence of a single variable on the long-term carbon cycle, such as mantle temperature (Sleep & Zahnle 2001), surface temperature (Foley 2015) or decarbonization fraction in subduction zones (Höning et al. 2019b).

In this work we assessed the effects of the planetary interior on long-term carbon cycling facilitated by plate tectonics. More specifically, we investigated the role of radioactive isotope abundance, core size and mass on the evolution of the atmospheric partial CO₂ pressure. In addition, we assessed the range of these parameters for which plate tectonics provides an efficient mechanism for carbon recycling, while we also derived an indication of the minimum age of a planet required for its atmospheric CO₂ pressure to be completely regulated by the long-term carbon cycle. For these purposes, we developed a two-component model which connects a parametrized mantle convection model to a box model describing the evolution of the long-term carbon cycle. Plate tectonics here acts as the interface between thermal evolution and carbon cycling, where we parametrized the mean plate speed in terms of the Rayleigh number (Schubert et al. 2001; Turcotte & Schubert 2014).

We elaborate on our model setup in Sect. 2, where we derive a parametrization for the plate speed to connect mantle convection with carbon cycling. Subsequently, we present our key results for different radioactive isotope abundance, core size and planet mass in Sect. 3, whose validity is discussed in Sect. 4. In Sect. 5 we discuss the implications of our results for the long-term evolution of atmospheric CO₂ pressure on Earth-like exoplanets. Finally, we summarize our key conclusions in Sect. 6.

2 Model setup

In this section, we present condensed discussions of the thermal evolution and carbon cycling model in Sects. 2.1 and 2.2, while a full technical discussion is provided in Appendices A and B, respectively. Subsequently, we motivate and derive a scaling law for the plate speed as the main link between interior and carbon cycle in Sect. 2.3, while Sect. 2.4 presents the scaling laws required to subsequently extrapolate our model to planets of different core size and mass.

Fig. 1 Evolution of the upper mantle temperature Tm for different surface temperatures Ts and initial upper mantle temperatures Tm0.

2.1 Thermal evolution model

The key aim of our thermal evolution model is to provide a first-order estimate of the evolution of the vigor of mantle convection over time, which both directly and indirectly depends on mantle temperature (e.g. Schubert et al. 2001). Therefore, usage of a parametrized mantle convection model is appropriate. Our model is similar to the models considered in Driscoll & Bercovici (2014) and Höning & Spohn (2016). A full discussion of our model is presented in Appendix A.

We restrict ourselves to a two-component planetary composition; an iron core surrounded by a magnesium-silicate mantle, as these two planetary components are thought to be the main constituents of terrestrial planets (Seager et al. 2007; Zeng & Sasselov 2013). The mantle is assumed to consist of a single, isoviscous, convective layer. The convective mantle is heated externally from below by the core and internally by the decay of radioactive isotopes. Cooling occurs via heat loss though a conductive thermal boundary layer at the surface.

The mantle viscosity plays an important role, as it controls the long-term convective cooling of the mantle. This is demonstrated in Fig. 1, where the evolution of the upper mantle temperature T_m is shown for a set of different initial mantle temperatures (T_m0 = 1800, 1900, 2000 K), while the surface temperature is fixed at a constant value (T_s = 300, 400 K). Additional model parameter values are presented in Table D.1. The initial mantle temperature only plays a role in the first 1–2 Gyr, while later, the evolution of the mantle temperature is governed by the viscosity, independent of initial mantle temperature. On the other hand, the constant surface temperature T_s continues to play a role throughout the entire thermal history due to its effect on the surface heat flux q_u. Higher surface temperatures result in a hotter mantle.

2.2 Carbon cycle model

As a next step, we quantify the effects of the thermal evolution on the atmospheric CO₂ pressure with a box model of the carbon cycle. We here follow a modeling approach similar to, for example, Tajika & Matsui (1992), Sleep & Zahnle (2001), and Foley (2015). A more comprehensive discussion of the model itself is presented in Appendix B.

For the carbon cycle we consider the partition of a total carbon budget over carbon reservoirs for the atmosphere, ocean, oceanic crust, and mantle. The evolution of the carbon content of these reservoirs is found by calculating the carbon added and removed from the reservoirs at each timestep, while the total carbon budget remains conserved. We assume equilibrium of carbon between the atmosphere and ocean reservoir to be instantaneous, as the timescale for equilibrium between these reservoirs is short (~ 10³ yr) compared to the other reservoirs (≳10⁶ yr) (Sleep & Zahnle 2001). Atmospheric and oceanic carbon is removed via both continental and seafloor weathering, and stored on the seafloor in the form of carbonate rocks. Subduction leads to carbon sequestration in the mantle, while a fixed fraction is degassed in the atmosphere via arc volcanism. In addition, mantle carbon is degassed directly in the atmosphere at mid oceanic ridges. Subduction and ridge degassing both depend on the plate speed. Therefore, the rate at which carbon cycles through the interior is governed by the plate speed and the efficiency of the above outlined weathering processes.

2.3 Platespeed parametrization

In this study, we consider the plate speed v_p as the main coupling variable between the carbon cycle and mantle convection models. It affects the rate at which carbon is transported toward the mantle via subduction, and from the mantle via degassing at mid-oceanic ridges. We motivate and derive here a parametrization for the plate speed in terms of the mantle global Rayleigh number and mantle thickness. The behavior of other couplings and their potential effects on the results of this study are discussed in Sect. 4.3.

In the boundary layer model used to develop the thermal evolution model (Appendix A), one can express the mean horizontal flow speed u₀ of mantle fluid in the thermal boundary layer near the planetary surface in terms of the mantle Rayleigh number (Schubert et al. 2001; Turcotte & Schubert 2014): (1)

Here, D denotes the mantle thickness, while β denotes the Nusselt-Rayleigh coupling exponent, and describes the relation between mantle convection and the mantle cooling. Mantle convection is quantified by the Rayleigh number Ra, while mantle cooling is described by the Nusselt number Nu, the ratio between total heat flow and conductive heat flow (Turcotte & Schubert 2014). β connects the Rayleigh and Nusselt number via (2)

Theoretical considerations of Rayleigh-Bénard convection with an isoviscous mantle fluid result in for Earth’s mantle (Turcotte & Schubert 2014). However, the value of β may be lower depending on mantle rheology (Schubert et al. 2001; O’Neill 2020 and references therein). How the precise value of β changes in the mantles of Earth-like exoplanets therefore remains unknown. Throughout this study, a value of is used, while we explore the effects oflower values of β on our model results in Sect. 4.2.

We approximate the Rayleigh number with the global Rayleigh number of the mantle (3)

Here, g denotes the mean mantle gravity, ρ_m the mean mantle density, c_m the mantle heat capacity, α the thermal expansivity, and k the thermal conductivity. In addition, T_c denotes the temperature at the core-mantle boundary (see also Fig. A.1), T_b the temperature above the thermal boundary layer between the mantle and core, T_m the mantle temperature below the upper thermal boundary layer, and T_s the surface temperature. η = η(T_m) is the mantle dynamic viscosity, whose dependence on mantle temperature T_m is described by the Arrhenius law given in Eq. (A.9).

As the upper thermal boundary layer directly corresponds to the lithosphere, we assume that the plate motion is coupled to the motion of the mantle fluid, which implies that the plate speed is linearly coupled to the fluid speed: (4)

Although this relation is generally applied to describe plate tectonics on present-day Earth (Schubert et al. 2001), the dynamics of the plateswith respect to the mantle fluid in reality depends on the material properties of the lithosphere and mantle. In this context, Crowley & O’Connell (2012) found Eq. (4) to be applicable whenever the material strength of the plates is low, such that the plate motion is fully controlled by the convective properties of the mantle. In contrast, if the plate material strength becomes high or the mantle viscosity low, other regimes exist for the plate speed behavior. In these regimes, the plate speed is less sensitive to the properties of the mantle. Instead of mantle properties, lithospheric material properties become an important or even dominant factor in controlling the plate dynamics. Therefore, Eq. (4) is valid only if the mantle and lithosphere properties are similar to the present-day Earth. We discuss the validity and implications of Eq. (4) more elaborately in Sect. 4.3.

Using Eqs. (1) and (4), we write the following scaling relation for the plate speed: (5)

Here D_⊕ and Ra_⊕ denote the present-day Earth mantle thickness and Rayleigh number, respectively, and are given in Table D.1. The value of Ra_⊕ is inferred from a run of the uncoupled thermal evolution model with T_s = T_s⊕ = 285 K, T_m0 = 2000 K, and T_c0 = 3500 K. The resulting behavior of the plate speed as a function of global Rayleigh number is presented in Fig. 2.

It follows from Eqs. (3) and (5) that the plate speed has a strong dependence on mantle temperature T_m through the mantle viscosity η, as given by Eq. (A.9). In addition, the plate speed also depends on core temperature T_c via Eq. (3). Therefore it is to be expected that the initial mantle temperature T_{m 0} and core temperature T_c0 leave their mark on the evolution of the mean plate speed v_p. In order to assess the role of initial interior temperatures in more detail, we perform model test runs where the surface temperature is kept fixed at T_s = 285 K, and show the evolution of the plate speed and mantle temperature for T_m0 = 1600, 1800, 2000 K and T_c0 = 3500, 4000 K in Fig. 3. In the long term, all initial conditions result in a gradually declining plate speed at t ≳ 1 Gyr. Furthermore, the evolution of the plate speed up to t ~ 2 Gyr depends on both the initial mantle and core temperature, while their influence becomes significantly less pronounced at later times. Figure 3 also demonstrates that the evolution of the plate speed is primarily determined by T_m, both due toits effect on the temperature gradient in Eq. (3) and the strong effect of T_m on the mantle viscosity. On the other hand, the evolution of the core temperature T_c affects the evolution of the plate speed both directly via the temperature gradient, and indirectly through the modified evolution of mantle temperature T_m.

In addition to the interior-carbon cycle coupling through the mean plate speed, our model contains another coupling between the CO₂ pressure and the thermal evolution of the interior via surface temperature T_s (Fig. 4). However, initial tests revealed that large (≳100 K) differences in surface temperature must be maintained for prolonged periods of time (~ 1 Gyr) to induce significant differences in the long-term evolution of the mantle temperature and plate speed. Therefore, we do not further discuss the role of this coupling in the rest of this work.

Fig. 2 Plate speed scaling law given in Eq. (5) as a function of the Rayleigh number Ra, along with the thickness of the upper thermal boundary layer (TBL). Values for present-day Earth are indicated for comparison.

2.4 Scaling mass and core size

The planet mass and core size have a large influence on the thermal structure and evolution of the mantle (Valencia et al. 2006). Therefore, we aim to assess how variations in these two planetary properties modify the long-term carbon cycle. For this purpose, we introduce parameter scalings in this section to extrapolate the model discussed so far toward planets with a core size and mass different from Earth.

An important variable in this context is the mean planetary density , which is a function of planet composition (Zapolsky & Salpeter 1969). We note that the density depends on the equation of state of the mantle and core material, and is thus a function of the temperature profile T(r) and pressure profile P(r) of the mantle and core (Valencia et al. 2006; Unterborn et al. 2016). For simplicity, however, we impose a mean mantle density ρ_m and mean core density ρ_c representative for the present-day Earth mantle and core. This simplification limits the parameter space that can be probed reliably. Even more restrictive is the isoviscous approximation, for the mantle, which forces us to consider only planets in the mass range from 0.1 M_⊕ (about the mass of Mars, de Pater & Lissauer 2015) to 2 M_⊕. We elaborate on this upper boundary of the planetary mass range in Sect. 4.1. For the core mass fractions f_c, we consider values between 0.1 and 0.9; an interval which is motivated from the numerical instability of our model at large and small core mass fractions. We note however that pure silicate (f_c = 0) or pure iron (f_c = 1) exoplanets may exist (Elkins-Tanton & Seager 2008; van Hoolst et al. 2019).

The mean mantle density and core density are calculated via (6)

Here, f_c⊕≈ 0.3259 denotes Earth’s core mass fraction (Valencia et al. 2006) and M_⊕≈ 5.97 × 10²⁴ kg Earth’s mass. Furthermore, V_m⊕ and V_c⊕ represent the volume of Earth’s mantle and core, respectively, whose values along with the resulting densities are given in Table D.1.

We approximate the mean mantle gravity g with the planet’s surface gravity, given by . In addition, we scale the total carbon budget R_tot with mantle mass (7)

Here, we choose R_tot,⊕ = 2.5 × 10²² mol as the reference value for Earth (Sleep & Zahnle 2001). Furthermore, we set the total length of mid-ocean ridges or subduction zones L proportional to the planet radius R_m.

Another aspect to consider upon extrapolating is the initial mantle and core temperature throughout M_p -f_c parameter space. Figure 3 shows that the choice of initial mantle temperature plays a large role in the early evolution of the plate speed. The initial temperature throughout the planetary interior is primarily set by an energy balance during accretion between the release of gravitational energy upon the formation and the radiative cooling of the outer surface of the planet (van Hoolst et al. 2019). Subsequently, differentiation provides additional interior heat. In our model, we assume an Earth-like lithosphere and core to be already present at t = 0. Temperature gradients throughout the lithosphere have been found to be only a weak function of planet mass; (Valencia et al. 2007). Using this scaling relation, one finds the upper mantle temperature T_m in a 10 M_⊕ super Earth to differ less than 5% with respect to a 1 M_⊕ (otherwise identical) planet. Altogether, we choose an initial upper mantle temperature of T_{m 0} = 2000 K as a first order approximation to the initial mantle temperature, following for example Noack & Breuer (2014). The latter also impose a temperature increase ΔT_c of roughly 1000 K between the bottom mantle temperature T_b and temperature at the core-mantle boundary T_c. We impose a similar initial temperature gradient through the lower thermal boundary layer as an initial condition. Altogether this allows us to write T_c0 in terms of T_m0, using that T_m increases adiabatically to T_b as given by Eq. (A.5): (8)

With this expression, we find values for T_c0 that are comparable with the values for T_c0 used in Noack & Breuer (2014). However, we note that the value of ΔT_c and also the temperature profile relating T_m to T_b in reality depend on the physical and chemical properties of the mantle. For example, the thermal structure depends on the local behavior of the viscosity, which is in reality not only a function of temperature, but also of mantle material composition and pressure (Stamenkovic et al. 2011, 2012). The high pressures thought to prevail in the lower mantles of massive super Earths could substantially increase the local viscosity. This would in turn impair the efficiency of convection as a heat transport mechanism, and could lead to super-adiabatic temperature profiles in the deeper mantle (Tackley et al. 2013). For simplicity, these effects are not considered in our thermal evolution model. Instead, we limit our study to planets where we expect these effects to be negligible to first order. Therefore, we consider planets up to a maximum mass M_p≤ 2 M_⊕; a limit which we motivate more elaborately in Sect. 4.1. We note that the dependence of T_c0 on planet mass and core size are contained implicitly in Eq. (8) via the mantle thickness D ≡ R_m − R_c and mean mantle gravity g. In addition, D and g directly affect the Rayleigh number and plate speed through Eqs. (3) and (5). Therefore, by varying planet mass and core mass fraction, we effectively vary the mantle thickness and mean mantle gravity in this study.

We explore the behavior of T_c0, D and g as a function of planet mass and core mass fraction in Fig. 5. The initial core temperature remains between 3000 and 4000 K for most of the parameter space considered (0.1 < f_c < 0.9; 0.1 M_⊕≤ M_p ≤ 2 M_⊕). An exception exists for more massive planets with a core fraction below 0.5, where T_c0 can exceed 5000 K with f_c approaching 0.1. This behavior is a result of the dependence of T_c0 on mantle thickness. The variation of the mantle thickness, shown in the center panel of Fig. 5, displays a rapid increase in mantle thickness at small f_c. Though this effect applies to any planet due to the assumed mass-radius relation, the effect is more pronounced for more massive planets due to their larger size. Mantle gravity, which also affects T_c0, is larger for massive planets with large core mass fractions.

Fig. 3 Evolution of plate speed vp (left) and mantle temperature Tm (right) for different initial mantle and core temperatures.

Fig. 4 A visual summary of the model considered in this study.

Fig. 5 Initialcore temperature Tc0 (left), mantle thickness D (center), and surface gravity g (right) as a function of core mass fraction fc and planet mass Mp.

3 Results

This section presents the findings of our model. As a first step, we assess how the distribution of carbon over the various reservoirs evolves over time under the influence of mantle cooling in Sect. 3.1. Subsequently we explore the effects of different interior properties on the long term carbon cycle: mantle radioactive isotope abundance (Sect. 3.2), core size (Sect. 3.3) and planet mass (Sect. 3.4). Thereafter, we investigate the feasibility of carbon cycling mediated by plate tectonics throughout this parameter space in Sect. 3.5, where we introduce the equilibrium timescale as a diagnostic tool.

3.1 Basic behavior

Figure 6 shows the effects of mantle cooling on the long-term carbon cycling, where the evolution of the plate speed, surface temperature and carbon reservoirs are presented for three different initial distributions of carbon over the various reservoirs. In addition, we investigate differences of our results with respect to Foley (2015), where the plate speed is assumed to be a function of surface temperature. Test runs revealed that differences between the evolution obtained from the model of Foley (2015) and a model where the plate speed is fixed at its present-day value of 5 cm yr⁻¹ are negligible in most cases. Noticeable differences were only found when the surface temperature is very high, resulting in lower plate speeds in the model of Foley (2015). The total budget of carbon participating in the carbon cycle is kept constant at R_tot = 2.5 × 10²² mol, while the initial conditions consist of a case where all carbon is initially in the surface reservoir R_s0 = R_a0 + R_o0 = R_tot, in the mantle reservoir R_m0 = R_tot, and an intermediatecase where . Furthermore, we assume T_m0 = 2000 K and T_c0 = 3500 K as initial conditions for the thermal evolution model. These initial conditions serve as our baseline initial conditions throughout the rest of this work, unless stated otherwise. The solar irradiation S_irr (t) is kept fixed at its present-day value S_⊙ in order to isolate the effects ofS_irr = S_irr(t) and by the cooling interior on the CO₂ pressure and T_s. A discussion of the effects of S_irr(t) is presented in Appendix C. All model parameters are presented in Tables D.1 and D.2.

For the two cases with initially CO₂-rich atmospheres, the evolution of the atmospheric CO₂ pressure canbe classified into three phases: Initially, weathering removes carbon from the atmosphere and it gradually accumulateson the seafloor (Figs. 6B, D). Despite the initially high surface temperatures (up to 500 K), liquid water, a requirement for weathering, is still present on the surface due to the higher surface pressure resulting from the thick, CO₂-rich atmosphere (Foley 2015). However, due to the high surface temperature and atmospheric CO₂ pressure, continental weathering is initially supply-limited: the weathering rate is limited by the physical erosion rate, and not on the kinetics of the weathering reaction (e.g. Foley 2015, and references therein). This means continental weathering initially proceeds at a rate , constant with respect to surface temperature and atmospheric CO₂ pressure (see Appendix B). After 10⁶ –10⁷ yr, equilibrium between arc volcanism and weathering is established, while carbon continues to build up in the mantle reservoir (Fig. 6E) until equilibrium between both ridge degassing and arc volcanism and weathering is established after around 10⁸ yr. In this last step, the plate speed plays a pivotal role as it governs the rate at which carbon can be sequestered in the mantle. In the case where the carbon budget is fully in the mantle, the plate speed is also the limiting factor which determines the time required for equilibrium between the various reservoirs to be established as it controls degassing from the mantle reservoir.

We initially find the plate speed to be significantly larger than its present-day Earth value (Fig. 6A) (~ 90 cm yr⁻¹ compared to ~ 5 cm yr⁻¹) due to the initially hot mantle in our model. This causes the time until equilibrium between the reservoirs is established to decrease by an order of magnitude (from ~10⁹ to ~10⁸ yr) with respect to the model of Foley (2015).

3.2 Influence of radiogenic heating

As a next step, we investigate the effects of radiogenic heating on carbon cycling. For this purpose we consider the long-term evolution of partial CO₂ pressure and mantle temperature T_m in Fig. 7 for different abundances of uranium and thorium, which scale linearly with their specific heat productions (parameters q_i in Table D.1). The initial distribution of carbon is set to R_s0 = R_tot, although the effects of this initial condition become negligible within the first 10⁸ yr. Furthermore, all initial conditions and parameter values are kept fixed with respect to Sect. 3.1. For scaling the radiogenic heating provided by the decay of thorium, we use the constraints derived by Unterborn et al. (2015), and scale the thorium abundance f_Th between 0.6 and 2.5 Earth’s thorium abundance. We also apply these limits for the total uranium abundance f_U = f_235U + f_238U as the nucleosynthetic origin of thorium and uranium are correlated (Goriely & Arnould 2001; Frebel et al. 2007; Unterborn et al. 2015).

Figure 7 reveals that increased radioactive isotope abundances result in a higher mantle temperature and hence CO₂ pressure due to the effects of mantle temperature on the plate speed via the Rayleigh number. In addition, a high abundance of uranium results in a more rapid decline of over time. The mantle temperature appears to display a stronger decline for f_U = 2.5, due to the larger contribution of uranium to mantle heating at earlier times. The evolution of follows this more rapid mantle cooling. In addition, differences between models with differing radioactive isotope variations tend to grow over time. Quantitatively, this amounts at t = 10 Gyr to a temperature ranging from 1500 to 1680 K for the coolest and hottest mantle, respectively.This range of mantle temperatures mentioned above translates into a partial CO₂ pressure ranging from 5 Pa to45 Pa (t = 10 Gyr); almost one order of magnitude. It also becomes clear that variations in the abundance of thorium have a larger effect on the thermal evolution and resulting atmospheric CO₂ pressure than uranium, which can be explained by thorium’s long half-life and comparatively high abundance in Earth’s crust. Enrichment or depletion of uranium in the mantle has a similar effect on the mantle temperatures and , albeit less pronounced, in particular at later times. This effect results from the shorter half-life of the uranium isotopes.

As a next step, we investigate the implications of a different isotopic abundance for the planetary climate. In Fig. 8, we consider snapshots att = 4.5 Gyr of , T_s and T_m as a function of f_Th and f_U. Both abundances are varied within the constraints provided by Unterborn et al. (2015). For f_Th = f_U = 2.5 we find a mantle temperature of about 1790 K, compared to 1560 K for the case where f_Th = f_U = 0.6. This results in ranging from ~14 Pa to ~ 134 Pa. The average surface temperature ranges from 283 to 287 K. Furthermore, Fig. 8 reveals that T_m does not increase linearly as a function of radioactive isotope abundance. Instead, the mantle temperature (and hence and T_s) appears to increase more slowly if the amount of radiogenic heating is larger. Overall, a higher isotopic abundance results in a hotter mantle. This in turn enhances the vigor of convection in the mantle due to the decreasing viscosity, resulting in enhanced degassing. Therefore, one can expect to find more long-term accumulation of CO₂ in atmospheres of planets with high abundances of long-lived radioactive isotopes.

Fig. 6 Evolution of plate speed vp (A), atmosphere reservoir content Ra (B), surface temperature Ts (C), crust reservoir content Rk (D) and mantle reservoir content Rm (E) for different initial distributions of the total carbon budget Rtot over the various reservoirs. Different color depict model runs with different initial surface reservoirs Rs (t = 0) = Ra(t = 0) + Ro(t = 0). Solid lines indicate results from this work, while dotted lines represent results from Foley (2015).

Fig. 7 Evolution of partial CO2 pressure (left) and upper mantle temperature Tm (right) for different abundances of uranium (235U and 238U) and thorium (232Th). Abundances are expressed relative to Earth (fU = fTh = 1).

Fig. 8 Present-day partial CO2 pressure (left), surface temperature Ts (center) and mantle temperature Tm (right) as a function of thorium and uranium abundance relative to Earth.

3.3 Influence of core size

In Fig. 9, we explore the effects of the core size on carbon cycling by taking snapshots at t = 0.1, 1, 2, 4.5 and 10 Gyr of the mantle temperature, Rayleigh number and atmospheric CO₂ pressure as a functionof core mass fraction f_c. The core mass fraction is varied between 0.1 and 0.9. This parameter range is motivated from the numerical stability of our model. The former fraction corresponds to a mostly silicate planet which consists largely of a mantle. It has been argued that coreless planets may exist (Elkins-Tanton & Seager 2008). On the other hand exoplanets with an inferred mean density consistent with pure iron have been found as well, which motivates our choice of the upper boundary (van Hoolst et al. 2019). All planets considered are assumed to have an initially CO₂ rich atmosphere (initially, R_s = R_a + R_o = R_tot), and an initial mantle temperature of 2000 K.

The mantle temperatures shown in Fig. 9 indicate that the mantles of planets with massive cores cool faster; a result from their larger surface to volume ratio. The faster cooling results in a rapidly increasing viscosity which reduces the Rayleigh number via Eq. (3). Additionally, a thinner mantle automatically results in a lower Rayleigh number due to the cubic dependence of the Rayleigh number on mantle thickness D. Altogether this results in a rapidly declining Rayleigh number for core fractions f_c ≳ 0.8; a trend followed by due to the dependence of the plate speed on the Rayleigh number. As the Rayleigh number scales with D³, the plate speed scales as for . Therefore, the plate speed also decreases sharply at large core fractions. As a result, the rates of seafloor weathering, ridge volcanism and subduction also become severely limited, and the mantle reservoir decouples from the surface reservoir. On the other hand, continental weathering continues to transport CO₂ from the atmosphere and ocean to the seafloor, resulting in a strongly decreasing atmospheric CO₂ pressure at large core fractions. We note this result is comparable to the conclusions of Noack et al. (2014) for stagnant lid planets. In addition, the lower mantle temperature and more sluggish plate movement at large core fractions could imply that plate tectonics itself is also more difficult to initiate on those planets in the first place. Altogether a large planetary core can be detrimental for the feasibility of a long-term carbon cycle facilitated by plate tectonics, as it severely disrupts the exchange of carbon between the mantle and surface.

Fig. 9 Snapshots over time of upper mantle temperature Tm (left), the global Rayleigh number (center) and partial CO2 pressure (right) as a function of core mass fraction fc.

Fig. 10 Snapshots over time of upper mantle temperature Tm (left), the global Rayleigh number (center) and partial CO2 pressure (right) as a function of planet mass Mp.

3.4 Influence of planet mass

In Fig. 10 we consider snapshots of the mantle temperature, Rayleigh number and CO₂ pressure as a function of planet mass M_p at different timesteps. We consider planets with a mass of 0.1 M_⊕≤ M_p ≤ 2 M_⊕; a parameter range which we discuss in more detail in Sect. 4.1.

We find the atmospheric CO₂ pressure to increase toward more massive planets due to the increasing mantle thickness. This results in the mantle of a 2 M_⊕ planet at t = 4.5 Gyr to be ~ 30 K hotter than that of a 1 M_⊕ planet. The lower viscosity results in a larger Rayleigh number, which promotes degassing via the plate speed. The resulting CO₂ pressure at t = 4.5 Gyr ranges from ~30 Pa to ~110 Pa for a 1 M_⊕ and 2 M_⊕ planet, respectively.

Toward less massive planets, the mantle thickness decreases strongly, as demonstrated in Fig. 5. Therefore, the mantle is able to cool more efficiently. This is in particular true for planets with M_p < 0.5 M_⊕, below which we see an accelerated decrease in mantle temperature. A 0.1 M_⊕ planet has a mantle which is roughly ~90 K cooler than a 1 M_⊕ planet at t = 4.5 Gyr, resulting in a partial CO₂ pressure of only ~0.3 Pa.

For the overall efficiency of the carbon cycle, however, variations in planetary mass appear to have more moderate effects than the core size. This is in particular true for planets which are more massive than Earth. This is a result from the assumed mass-radius relationship , which causes the mantle surface-volume ratio A_m∕V_m to decrease more slowly toward massive planets. However, we note this cooling behavior as a function of M_p is a matter of debate, and may depend on the Nu-Ra coupling parameter β (Eq. (2)) (Stevenson 2003; Seales & Lenardic 2020, 2021). Altogether, we find that the effects of planetary mass on the thermal evolution and hence carbon cycle are more moderate toward more massive planets.

3.5 Equilibrium timescales

Figure 6 shows that after ~10⁸ yr, the vast majority of the total carbon budget is stored in the mantle reservoir for any initial distribution of carbon over the various reservoirs, even if the atmosphere is initially carbon-rich. Such an initially CO₂-rich atmosphere is thought to be a common consequence of the initial solidification of the crust (Elkins-Tanton 2008). This subsequently gives rise to surface temperatures well above the moist greenhouse threshold of 340 K, allowing for large amounts of water loss if those temperatures are maintained long enough (≳ 10⁸ yr) (Kasting 1988; Abbot et al. 2012). If the timescale on which the carbon cycle establishes equilibrium between the various carbon reservoirs is >10⁸ yr, this maybe disadvantageous for the development of a long-term, temperate climate with moderate atmospheric CO₂ pressures. In this section we explore this equilibrium timescale as a function of the interior parameters considered in the previous sections. We also consider this timescale as an indicator of the absolute minimum age an exoplanet must have in order to have a climate which is significantly impacted by a long-term carbon cycle facilitated by plate tectonics.

An estimate for the carbon cycle equilibrium timescale τ_CC can be derived by considering the evolution of the carbon reservoirs R_m, R_k and R_a for different initial carbon distributions. The initial conditions considered are the same as in Sect. 3.1. Subsequently, τ_CC is chosen to be the time where all three different evolutionary tracks of the three reservoirs R_i all deviate no more than 1% from the mean reservoir content inferred from the three different values of R_i(t) for the different initial conditions: (9)

This expression is evaluated for all three reservoirs, R_i ∈ (R_m, R_k, R_a), and the longest time found is chosen as the equilibrium timescale of the given model. With this definition, one can find τ_CC ≈ 116 Myr for the model representing Earth, discussed in Sect. 3.1.

In the context of mantle abundances of uranium and thorium, we find that τ_CC varies moderately (Fig. 11). Values of τ_CC here range from ~100 Myr for f_U = f_Th = 0.6 to ~ 120 Myr for f_U = f_Th = 2.5. The lower value of τ_CC toward mantles rich in radioactive isotopes is a consequence of the higher plate speeds associated with a hotter mantle, which allows for more rapid sequestration of carbon in the mantle reservoir. The variation in τ_CC is moderate due to the fact that equilibrium is already established when the mantle temperature is still close to its initial value, which was chosen to be T_m0 = 2000 K for all models. Therefore the differences in the thermal budgets induced by radiogenic heating are still small, which propagates into the plate speed and τ_CC.

For the core mass fraction f_c, the equilibrium timescale increases by over an order of magnitude over the range of core mass fractions considered (0.1≤ f_c ≤ 0.9) if the planet mass is fixed at M_p = 1 M_⊕ (Fig. 12). This increase is the consequence of two effects, the first being the comparatively low thickness of the mantle D at large core fractions. As the plate speed scales linearly with mantle thickness, planets with a thinner mantle have lower initial plate speeds, as the mantle temperature is initially nearly identical for all models. The second effect is the high surface-volume ratio of the mantle. This allows for efficient mantle cooling at larger core fractions (see also Fig. 9), reducing the Rayleigh number and plate speed even further over time. It becomes clear that on planets with large f_c, carbon exchange with the mantle is limited due to the low rates of plate subduction and seafloor spreading.

As a function of planet mass at constant f_c = f_c⊕ = 0.3259, τ_CC decreases from τ_CC ≈ 290 Myr for a 0.1 M_⊕ planet to τ_CC ≈ 100 Myr for a 2 M_⊕ planet, a consequence of the same two effects. More massive planets have mantles which are thicker and hotter, resulting in more vigorous convection, and hence higher plate speeds. This leads to a higher rate of carbon sequestration into the mantle, and hence results in a shorter equilibrium timescale for more massive planets. It becomes clear that carbon exchange with the mantle facilitated by plate tectonics is efficient for a wide range of planetary interiors.

In addition to the physical properties of the planet considered, we note that the equilibrium timescale also depends on the parameter values used in our model. An important parameter in this context is the efficiency of carbon sequestration, which is contained inf, the fraction of carbon which degasses via arc volcanism upon subduction (see also Appendix B). A high value of f means a large fraction of the subducted carbon is degassed during subduction. This limits the rate at which carbon is sequestered in the mantle, increasing τ_CC. In our study, we use the value f = 0.5 (Table D.2), although estimates for f vary from below 0.1 up to values as large as 0.7 (Sleep & Zahnle 2001; Dasgupta & Hirschmann 2010; Ague & Nicolescu 2014; Kelemen & Manning 2015; Foley 2015). For our model, we find τ_CC to vary from ~60 to ~200 Myr within this range, as is shown in Fig. 13. This variation is significant when compared with the variation in τ_CC found for different radioactive isotope abundances, planet mass and core mass fraction. The uncertainty in the arc volcanism degassing fraction f thus induces a considerable uncertainty in estimates of the equilibrium timescale.

Another parameter which could affect the values found for τ_CC is the exponent b, which controls the direct dependence of continental weathering on the atmospheric CO₂ pressure (Eq. (B.6)). The value of b = 0.55 used in this study (Table D.2) follows from Foley (2015) and Driscoll & Bercovici (2013). On one hand, this is considerably larger than values used by Walker et al. (1981) and Zahnle & Sleep (2002) (b = 0.3), or by Berner (1991) (b = 0.22). On the other hand, these values do not imply biological enhancement of weathering (e.g. Schwartzman & Volk 1989). Since biological productivity increases with temperature and atmospheric CO₂, a potential biosphere amplifies the CO₂-dependence of weathering and thereby climate stability (Caldeira & Kasting 1992; Höning 2020), which would effectively increase the value of b. In any case, a lower value of b weakens the feedback provided by continental weathering, and thereby increases the equilibrium timescale. Figure 13 shows that lower values of the weathering exponent b result in a longer equilibrium timescale, increasing from τ_CC ≈ 115 Myr for b = 0.55 to τ_CC ≈ 127 Myr for b = 0.2. This is smaller than the variation found in Fig. 11 for different mantle radioactive isotope abundances, and well bellow the variation in τ_CC found as a function of M_p and f_c. Altogether the effect of b on τ_CC remains limited.

Fig. 11 Equilibrium timescale τCC as a functionof mantle abundance of uranium (235U and 238U) and thorium.

Fig. 12 Equilibrium timescale τCC as a functionof core size fc and planet mass Mp.

Fig. 13 Equilibrium timescale τCC as a functionof weathering exponent b and arc degassing fraction f.

4 Discussion

4.1 Extrapolation to massive planets

Exploring the effect of planet mass on the carbon cycle, we have assumed the mass-radius relation , and neglected any effects of the pressure or temperature on the density. However, these assumptions will limit the range of planetary masses that can be meaningfully explored with this model. Valencia et al. (2006) find a mass-radius relation for planets with a Mercury-sized core and for an Earth-sized core, with the precise value of the exponent also depending on the equation of state of mantle and core. They also note that for planets less massive than Earth, the effects of pressure and temperature on the density of the mantle and core becomeless important, resulting in a convergence of the exponent in the mass-radius relationship toward , the value associated with constant density. This suggests that for planets less massive than Earth, the constant-density assumption is acceptable.

Toward more massive planets, we impose that the relative difference between the radius estimates must be smaller than typical uncertainties associated with observational radius measurements. In this context, Otegi et al. (2020) use a relative uncertainty of σ_R ∕R ≤ 8% as a selection criterion for robust radius measurements of exoplanets. The deviation of with respect to the relation exceeds this difference at M_p ≈ 3.6 M_⊕.

Besides the effects of pressure on the mean density, the increasing pressure in the interior of massive planets has also effects on the mantle viscosity. Through our work, we assumed the mantle to be isoviscous. However, from literature, conclusions for the viscosity profile throughout the mantle seem to differ for higher-mass planets. Some studies suggest that the isoviscous approximation remains valid for more massive planets (Karato 2011), while othersfind a strongly increasing viscosity along the same temperature profile throughout the mantle (Stamenkovic et al. 2011). This discrepancy may result from the lack of knowledge of the behavior of mantle rocks under high pressure (Stamenkovic et al. 2012). In general, these studies model the pressure-dependence of the viscosity viaextensions of the Arrhenius-type law for purely temperature-dependent viscosity used in our model (10)

Here V_eff(P) denotes the effective activation volume, which determines the coupling strength of pressure to viscosity, and is determined by the material properties of the mantle rocks. Altogether this may lead to an increasing viscosity as a function of mantle depth, which can have the formation of a stagnant lid at the core-mantle boundary as a result (Stamenkovic et al. 2012). The formation of such a lower stagnant lid is thought to ensue once the viscosity contrast between the top and bottom of the mantle exceeds values of 10⁴ (Solomatov 1995; Solomatov & Moresi 1997). The viscosity contrast in Earth’s mantle is estimated between 2 and 10, well below this range (Yamazaki & Karato 2001; Paulson et al. 2005; Driscoll & Bercovici 2014). On the other hand, planets with masses of M_p = 5 M_⊕ and M_p = 10 M_⊕ have inferred viscosity contrasts in the range of ~10⁵−10⁷ and ~10¹⁰−10¹⁴, respectively (Stamenkovic et al. 2011). A range of contrasts is stated as Stamenkovic et al. (2011) also consider variations in the temperature at the core-mantle boundary.

Altogether there appears to be no straightforward way in estimating the viscosity contrast in a mantle without imposing more detailed assumptionson the mineralogy of the mantles considered. However, the results from Stamenkovic et al. (2011) indicate that an upper boundary of M_p = 5 M_⊕ is without doubt too high for our isoviscous mantle model to produce reliable results. To be conservative, we limited ourselves to an upper mass limit of M_p = 2 M_⊕ for our parameter space. As temperature and pressure effects are thought to become less important for planets less massive than Earth, we limit the minimum mass considered in our parameter space to M_p = 0.1 M_⊕, about the mass of Mars (de Pater & Lissauer 2015).

4.2 The role of coupling exponent β

The Nu-Ra coupling exponent β has been kept fixed at in this study. This value is based on the assumption that the mantle can be treated as an isoviscous fluid undergoing Rayleigh-Bénard convection, which is also a key assumption underlying our thermal evolution model. β is an extremely important parameter as it describes the coupling strength between heat flow and convective vigor and therefore the feedback strength between mantle convection and mantle heat loss. However, the precise value of β in the mantles of Earth-like exoplanets is unknown. Therefore, we discuss which values for β are possible in the mantles of Earth-like exoplanets. In addition, we show the implications for the long-term evolution of the atmospheric CO₂ pressure.

While follows from theoretical considerations, laboratory experiments suggest β ≈ 0.3 to be appropriate for Earth’s mantle (Turcotte & Schubert 2014). However, β has been found to vary between 0 and depending on mantle rheology, thermal state and tectonic regime (Unterborn et al. 2015; O’Neill 2020 and references therein). Therefore, β could indirectly even be a function of time. For simplicity, however, we restrict our discussion to constant β, and consider the effects of values within the range of 0 to on mantle temperature T_m and in Fig. 14. Lowering the value of β weakens the effect of mantle convection on mantle cooling and plate speed. This means that the mantle temperature increases significantly over time for small β, up to above 3500 K for β = 0 at t = 10 Gyr. We note that at these temperatures, widespread mantle melt becomes an important cooling mechanism for the mantle. However, the cooling effects of melting have not been included in our model, but will in reality result in lower mantle temperatures than shown in Fig. 14 (Richter 1985; Driscoll & Bercovici 2014).

The evolution of the CO₂ pressure is also sensitive to the value of β due to its dependence on the mean plate speed v_p. The increasing mantle temperature increases the Rayleigh number via the viscosity, resulting in a higher plate speed and hence atmospheric CO₂ pressure for the intermediate values β = 0.1, 0.2, 0.3. However, a lower value of β also weakens the sensitivity of the plate speed on the Rayleigh number and hence mantle temperature. This explains why we find the highest atmospheric CO₂ pressure in Fig. 14 at intermediate values of β. β = 0 eliminates the dependence of the plate speed on the Rayleigh number entirely, resulting in the plate speed remaining fixed at its present-day Earth value throughout the entire evolution. It becomes clear that the applicability of our model coupling is limited for smaller values of β. However, given the similarities to Earth imposed on the planets considered in this study and mantle boundary layer model considered, appears to be a reasonable choice for this parameter.

Fig. 14 Long-term evolution of the upper mantle temperature Tm and partial atmospheric CO2 pressure for different values of the Nu-Ra coupling exponent β.

4.3 Platespeed scalings and additional couplings

In this study, we assumed a power-law relation between the mean plate speed v_p and global Rayleigh number Ra to connect the plate speed to mantle cooling. Furthermore, the plate speed was assumed to scale linearly with the flow speed of the mantle fluid. However, we note that other parametrizations for the plate speed exist, as the evolution of the plate speed does not necessarily have to be dominated by the thermal evolution of the mantle. In addition, the plate speed is not the only variable which couples the carbon cycle to the thermal evolution of the interior. In this section we discuss other parametrizations of the plate speed and their potential effects on our results. In addition, we highlight a few couplings which are not considered explicitly in this study.

In Sect. 3.1 we compared our results to the study of Foley (2015), where the plate speed was imposed to depend on surface temperatureT_s via the linear relation (11)

Here a₁ and a₂ are fitting parameters from earlier work (Foley & Bercovici 2014). However, we find that the plate speed does not vary significantly over the history of the planet, as depicted in Fig. 6. This different behavior of plate speed evolution with the evolution found with our parametrization is a consequence of the fact that the long-term behavior of the mean plate speed remains poorly understood. One class of models suggests faster plate tectonics in the past as a result of a hotter, more strongly convecting mantle (Höning & Spohn 2016). This could explain the observed D/H fractionation on Earth (Kurokawa et al. 2018). Other work suggests that the plate speed remained rather constant as a function of mantle temperature due to the effects of volatile depletion on mantle rheology (Korenaga 2003; Korenaga et al. 2017). However, the importance of these effects on the long-term evolution of plate tectonics on Earth is still a topic of debate, making an adequate description for these processes on Earth-like exoplanets challenging.

The long-term behavior of the plate speed in this work is assumed to be controlled by the Rayleigh number via Eq. (5). This expression relies crucially on Eq. (4), the assumption that the plate speed depends linearly on the mantle fluid speed. However, in case a planet has a high lithosphere strengths or a low mantle viscosity, other solutions for the plate speed exist in the so-called sluggish-lid regime. In this case, the plate dynamics is not fully coupled to the dynamics of the mantle fluid, such that v_p < u₀, and the material properties of the lithosphere control the movement of tectonic plates (Crowley & O’Connell 2012). The linear relation can be recovered in the presence of weak lithospheres and low material strength contrasts with respect to the mantle underneath. Altogether a linear relation between plate speed and fluid speed may not hold for exoplanets with mantle compositions and properties that differ significantly from Earth. In addition, a significantly hotter mantle than Earth’s may also violate Eq. (4), since the corresponding lower viscosity may decouple the plate motion from the mantle fluid flow (Valencia et al. 2007; Berovici et al. 2015). Contrary, it has also been argued that a hotter mantle may ultimately have a higher viscosity when the effects of dehydration on mantle rheology are taken into consideration (Korenaga 2003). Altogether, the precise behavior of plate movement as a function of temperaturewould require more detailed modeling of the mantle and lithosphere composition and their effects on the mantle viscosity and lithospheric strength with respect to Earth. However, it is important to emphasize that our results would not be affected significantly: the rate of mantle degassing directly depends on the rate at which mantle material enters the source region of partial melt, which is represented by the convection rate u₀ rather than by the plate speed v_p (c.f. Eq. (4)). The same applies to the rate of seafloor weathering, which is a function of the crustal production rate, not necessarily of the plate speed (see e.g. Höning et al. 2019b). Therefore, our results would still hold, even in a case where the plate speed speed does not linearly depend on the convection rate.

Other couplings between the thermal state of the mantle and the climate exist. For example, the decarbonation fraction in subduction zones has been studied as a function of the mantle temperature (Höning et al. 2019b). Similarly, these authors find declining atmospheric CO₂ pressures over time as a consequence of mantle cooling. Furthermore, the mid-ocean ridge melting depth d_melt has been argued to become larger in response to a hotter mantle (Sleep & Zahnle 2001). This would enhance the amount of ridge degassing at earlier times (Eq. (B.11)), while gradually decreasing over time as the mantle cools. Altogether the inclusion of the dependence of decarbonization fraction or melting depth on mantle temperature would amplify the declining trend in CO₂ pressure over time found in this study.

The removal of atmospheric CO₂ predominantly depends on the weathering rate F_weather. This rate is also indirectly connected to the plate speed. The efficiency of weathering is dependent on the amount of fresh, weatherable rock exposed to the atmosphere via the erosion rate, which is in the long term determined by the amount of uplift caused by plate tectonics (e.g. Kasting & Catling 2003; Foley 2015; Höning et al. 2019a). Therefore, one can expect weathering fluxes to decline for lower plate speeds. The coupling between continental weathering rate and thermal evolution, however, is left for future studies.

5 Implications for exoplanets

In this work we aimed to assess the effects of the mantle abundance of radioactive isotopes, core mass fraction and planet mass on the long-term carbon cycles of Earth-like exoplanets. In this section we therefore address the implications of our results for Earth-like exoplanets.

The coupling between mantle cooling and the plate speed suggests gradually decreasing plate speeds over time. As the planet ages and the mantle cools down, we find declining atmospheric CO₂ pressure. We note that this decline in CO₂ pressure likely amplifies the decline in CO₂ pressure in response to the increasing surface temperatures induced by the increasing luminosity of the host star (see Appendix C). This means that if a planet has a carbon cycle facilitated by plate tectonics, one would expect the atmospheres of rocky planets around older G-type stars to be depleted in CO₂ with respect to their younger counterparts. An example of such a star which has aged considerably is HD 16417. The age of this G1V-type star has been estimated at 7.0 ± 0.4 Gyr (Baumann et al. 2010), and is known to have at least one Neptune-sized (M_p = 22.1 ± 2.0 M_⊕) planet (O’Toole et al. 2009). If such a system were to have a terrestrial planet undergoing carbon cycling as known on Earth, its atmosphere could show an atmospheric CO₂ depletion, which could be an indicator of the presence of an Earth-like carbon cycle on such a planet. This could present an independent indicator for the existence of carbon cycling which could be distinguished with for example statistical comparative planetology (Bean et al. 2017; Graham & Pierrehumbert 2020).

Initially high plate speeds resulting from the hot interior allow for rapid exchange of carbon between the planetary surface and mantle, which allows for a rapid transition from a disequilibrium atmosphere to an atmosphere with its CO₂ pressure controlled by long-term carbon cycling. This transition occurs within 120 Myr for a wide range of radioactive isotope abundances, core sizes and planet masses. It indicates that if plate tectonics starts on a young, Earth-like exoplanet, its atmospheric CO₂ pressure could in principle become fully controlled by the long-term carbon cycle within 120 Myr. Therefore an exoplanet does not have to be old to have its atmospheric CO₂ pressure regulated by the long-term carbon cycle. This implies that around younger stars such as HD 141937 (1.3 ± 0.9 Gyr old, Baumann et al. 2010), one may find planets whose atmospheric CO₂ pressure is already controlled by carbon cycling. However, this is only possible if plate tectonics starts early and weathering is efficient at removing atmospheric CO₂.

Unterborn et al. (2015) found a spread in the abundance of thorium from 0.6 to 2.5 solar in the atmospheres of a sample of fourteen solar twins and analogs. These stellar atmospheric abundances were derived by means of spectral line fitting. As both thorium and uranium are refractory elements and have a correlated nucleosynthetic origin, we assumed that the abundances of these species in the mantles of terrestrial planets orbiting these stars scale linearly with the abundance of thorium in the atmosphere of the parent star. We found that after 4.5 Gyr, a mantle depleted in radioactive isotopes (f_Th = f_U = 0.6 Earth’s mantle abundance) can have a partial CO₂ pressure almost an order of magnitude below the partial pressure found on a planet with a mantle rich in radioactive isotopes (f_Th = f_U = 2.5). The abundance of thorium was here found to have a larger effect than the abundance of uranium, in particular at later times. This implies that planets around thorium-rich stars such as HD 102117, HD 141937 and HD 160691 (Unterborn et al. 2015) may favor the development of a warmer climate. The higher plate speeds associated with the hotter mantle also facilitate more efficient transport of carbon into and from the mantle. A mantle initially rich in radioactive isotopes also allows the mantle to remain hotter over longer time spans (Fig. 7), and thus maintain these higher plate speeds. This means that efficient carbon cycling is possible over a larger timespan. However, Fig. 11 shows that at earlier times, the equilibrium timescale does not depend as strongly on mantle radioactive isotope abundance, but rather on the initial mantle temperature. Altogether, rocky exoplanets in the plate tectonics regime around thorium- or uranium-rich stars may generally be more efficient in carbon cycling than Earth, in particular at later times.

The core size and mass of a planet also have important consequences for the carbon cycle facilitated by plate tectonics. The rapidly increasing equilibrium timescale in Fig. 12 as a function of core size suggests that this mode of carbon recycling is becoming less efficient on planets with a large core. The mantle reservoir starts to decouple from the surface due to the rapid mantle cooling (Fig. 9), with the equilibrium timescale approaching 1 Gyr for f_c ≳ 0.8. Altogether, plate tectonics operates less efficiently on planets with large core sizes; a conclusion consistent with Noack et al. (2014), who find that a planet is more likely to be in the stagnant lid regime when the iron core is large, making an Earth-like carbon cycle impossible. For planets that are less massive, we found similar behavior; small planets cool quickly and hence have more sluggish convection in their mantles, while a hotter interior for planets more massive than Earth would allow for faster carbon cycling on these planets.

6 Conclusions and outlook

In this work we studied the role of the planetary interior in carbon cycling facilitated by plate tectonics and the resulting evolution of the atmospheric CO₂ pressure. Interior properties considered are mantle radioactive isotope abundance, planet mass and core mass fraction, parameters which vary significantly for Earth-like exoplanets. In order to assess the effects of these interior parameters on the long-term carbon cycle, we developed a parametrized planetary evolution model which links the thermal evolution of the interior to carbon cycling. The mean plate speed is used as the key coupling variable between the two models. We extrapolated Earth’s long-term carbon cycle for a range of these interior parameters. Our findings can be summarized as follows:

• Including the effects of mantle cooling on the plate speed results in gradually declining atmospheric CO₂ pressures over time, due to the gradually decreasing degassing. This decline amounts up to an order of magnitude in over 10 Gyr; a conclusion qualitatively similar to Höning et al. (2019b). However, the erosion rate E_max also declines for lower plate speeds, and thus may cause the continental weathering rate to decline in response to lower plate speeds. Feedback between the plate speed and these variables could be important to consider in future studies, and assess whether this feedback stabilizes the atmospheric partial CO₂ pressure against decreased degassing.

• A long-term carbon cycle driven by plate tectonics could operate efficiently on planets with amounts of radiogenic heating in their mantles different from Earth. However, planets with their mantles enriched in radioactive isotopes with respect to Earth, may favor the development of warmer climates resulting from a more CO₂ rich atmosphere. This is in particular the case for a planet with a higher thorium abundance. In addition, the carbon cycle operates more efficiently on planets rich in radioactive isotopes, motivating the characterization of planetary systems around stars whose atmospheres are rich in thorium or uranium.

• A long-term carbon cycle regulated by plate tectonics may not be feasible for planets with core mass fractions f_c ≳ 0.8. This can mainly be attributed to the decrease of the Rayleigh number over time in response to fast mantle cooling. This effect reduces the effectiveness of plate tectonics in sequestering carbon in the mantle through subduction. Altogether this result emphasizes the importance of high-precision mass-radius measurements of Earth-sized exoplanets, which can help in constraining the core sizes of these exoplanets.

• Planets with plate tectonics may favor higher atmospheric CO₂ pressure for higher total mass or smaller core fractions than Earth’s. Both result in a hotter mantle with a higher Rayleighnumber, which promotes degassing through higher plate speeds. Future studies should focus on extrapolating these results toward planet masses above 2 M_⊕, which requires the incorporation of the pressure-dependence of the mantle viscosity and density in the thermal evolution model.

• Carbon cycling facilitated by plate tectonics can regulate the atmospheric CO₂ pressure from disequilibrium values comparatively fast; an equilibrium timescale between 100 and 200 Myr was found throughout most of the parameter space explored. This means that provided plate tectonic initiates, an exoplanet does not necessarily have to be old to have an atmosphere significantly altered by a long-term carbon cycle.

Altogether we find that the thermal evolution of the planetary interior has a significant effect on the evolution of the atmospheric CO₂ pressure. However, for a large fraction of the parameter space explored in this study, this does not significantly impede the efficiency at which the long-term carbon cycle is able to regulate the atmospheric CO₂ pressure.

Acknowledgements

We would like to thank Bradford Foley for his help in benchmarking the carbon cycle model, and an anonymous reviewer for constructive comments. D.H. has been supported through the NWA StartImpuls.

Appendix A Thermal evolution model

The evolution for the mantle temperature T_m and core temperature T_c are calculated by integrating the following energy conservation equations for the mantle and core over time (e.g. Schubert et al. 2001):

Here, V_m and A_m denote the volume and outer surface area of the mantle, while ρ_m and c_m its density and specific heat capacity. Furthermore, quantities with subscript c denote the same quantities for the core.

H(t) denotes the volumetric specific radiogenic heating per unit mantle volume provided by the radioactive isotopes ²³⁵U, ²³⁸U, ²³² Th and ⁴⁰K, (A.3)

Here Q₀ denotes thetotal present-day Earth heat production per unit volume. t_⊕ and t_h,i denote the age of Earth and half-life of radioactive isotope i, respectively.q_i is the percentage of Q₀ provided by radioactive isotope i.

In addition, q_u and q_l denote the heat fluxes through the upper and lower thermal boundary layer, respectively: (A.4)

Here k denotes the mantle thermal conductivity, T_s the surface temperature and T_b the temperature at the top of the lower thermal boundary layer (see Fig. A.1). T_b is related to T_m via the linearized adiabatic temperature profile (A.5)

Here, α denotes the thermal expansivity, g the mean gravitational acceleration throughout the mantle and D = R_m − R_c the mantle thickness. Furthermore, we assume the thickness of the lower thermal boundary layer δ_l ≂ δ_u. Furthermore, δ_u is dictated by a local instability criterion (Stevenson et al. 1983) where we impose that the local Rayleigh number Ra_u of the thermal boundary layer is equal to the critical Rayleigh number: (A.6)

Here, η(T_m) denotes the mantle viscosity. One thus finds for δ_u: (A.7)

For general Nu-Ra coupling exponent β (see Sect. 2.3), the above relation can be written as (Driscoll & Bercovici 2014) (A.8)

Fig. A.1 Overview of the temperature profile (blue) and thermal structure considered in the thermal evolution model

Furthermore, we assume the mantle to be isoviscous, with its viscosity being given by the Arrhenius law (A.9)

where E^⋆ denotes the activation energy and R_g is the ideal gas constant. η_⊕ denotes a reference viscosity at reference temperature T_m,⊕, for which we choose values representative for the present-day Earth upper mantle. An overview of all values used for the parameters introduced in this section is presented in Table D.1.

Appendix B Carbon cycle model

The carbon cycle model we use closely follows the model considered by Foley (2015), who developed this model to investigate the response of the carbon cycle to different amounts of land coverage and carbon participating in the cycle. We here present a technical overview of this model, while we refer to Foley (2015) for a full motivation and explanation of the model parameters used.

We consider the partition of a constant total carbon budget R_tot participating in the long-term carbon cycle over various reservoirs for the atmosphere R_a, ocean R_o, oceanic crust R_k and mantle R_m. R_i here denotes the amount of carbon present in reservoir i in moles.

Here, F_ridge denotes the carbon flux (in moles per unit time) due to degassing of carbon from the mantle at mid-oceanic ridges. F_arc = fF_sub is the degassing of CO₂ from subduction zones via arc volcanism, and F_sub denotes the carbon flux into subduction zones, with f being the fraction of the carbon which is degassed during subduction. We note that f is not well constrained (Sleep & Zahnle 2001; Dasgupta & Hirschmann 2010; Ague & Nicolescu 2014; Foley 2015) and depends on the thermal structure of subduction zones, which evolves over time (Johnston et al. 2011; Höning et al. 2019b). However, for simplicity, we keep f fixed at value of 0.5 throughout this work. The factor of in Eqs. (B.1) and (B.2) results from the partial release of carbon back into the atmosphere upon the deposition of carbon on the seafloor in the form of carbonate rocks (Kasting & Catling 2003).

The partition of carbon between the atmosphere and ocean is assumed to be instantaneous, and dictated by Henry’s law, which can be written as (B.4)

Here, denotes the partial atmospheric CO₂ pressure, and k_c the solubility of CO₂ in water. can be written in terms of R_a via (B.5)

denotes the molar mass of CO₂, g the surface gravity, and A_m denotes the surface area of the planet. Combination of Eqs. (B.4) and (B.5) allows one to infer R_a and R_o at each timestep, with R_a + R_o being calculated at each timestep via Eq. (B.1).

For continental weathering we distinguish between two weathering regimes: a kinetically-limited regime and a supply-limited regime (West 2012; Foley 2015). However, we note the existence of more recent treatments for the continental weathering rate which for example distinguish an additional “runoff-limited” regime (Graham & Pierrehumbert 2020; Hakim et al. 2020). In the kinetically-limited regime, the weathering rate is limited by the weathering reaction rate and hence depends on surface temperature T_s and partial atmospheric CO₂ pressure . However, if the weathering rate becomes very high, the weathering rate is limited by the rate at which weatherable rock is exposed to the atmosphere. We use a parametrization which captures the behavior of continental weathering in both regimes, following (West 2012; Foley 2015), (B.6)

Here, F_w⊕ denotes the present-day Earth weathering rate, f_land the fraction of the surface area covered by land, E_a the activation energy for the weathering reaction, R_g the ideal gas constant, and T_s the average surface temperature. Values denoted with ⊕-subscript indicate the respective parameter for present-day Earth. Furthermore, a and b are constants, and P_sat is the saturation pressure for water vapor, given by (B.7)

where P_sat0 denotes the reference saturation pressure for water vapor at temperature T_sat0. m_w denotes the molar mass of water and L_w the latent heat of water. The continental weathering rate is scaled with a saturation vapor pressure term in order to account for variations in runoff for different T_s (Driscoll & Bercovici 2013; Foley 2015).

When weathering becomes supply-limited, the dependence on T_s and disappears. Different parametrizations describing the continental weathering rate in the supply-limited regime exist (e.g. Riebe et al. 2004; Foley 2015). In this work, we use the expression derived by Foley (2015), (B.8)

Here E_max denotes the globally averaged erosion rate, for which we follow Foley (2015) and use an upper bound value which is kept constant over time. Furthermore ρ_r denotes the average regolith density, χ_cc is the fraction of reactable cations in the surface rock, and the average molar mass of molecules in the rock which participate in silicate weathering.

For seafloor weathering, we use a parametrization which depends on and plate speed v_p, although wenote that different parametrizations exist, for example in terms of the ocean floor temperature and ocean pH (Coogan & Dosso 2015; Krissansen-Totton & Catling 2017). In case of seafloor weathering dependent on and v_p, one can write the scaling law (Sleep & Zahnle 2001; Mills et al. 2014; Foley 2015) (B.9)

F_sfw⊕, v_{p ⊕} and here denote the present-day Earth seafloor weathering rate, plate speed and partial CO₂ pressure, respectively.

The subduction flux F_sub, can be written in terms of the surface area of seafloor entering subduction zones per unit time multiplied by the density of carbon on the seafloor, (B.10)

where L denotes the total length of subduction zones.

Lastly, the degassing flux at mid-oceanic ridges, F_ridge is given by (B.11)

Here, f_d denotes the fraction of the carbon in the upwelling mantle material that is released into the atmosphere, and d_melt the scaling depth at which mantle rocks start to melt due to depressurization (Schubert et al. 2001).

Lastly we relate to T_s. For this purpose we use the parametrization derived by Walker et al. (1981). (B.12)

Here, T_e the effective temperature and T_e⊕ the present-day Earth effective temperature. The effective temperature of a planet is given by (de Pater & Lissauer 2015) (B.13)

where S_irr denotes the solar flux, A the planetary albedo, and σ the Stefan-Boltzmann constant. We note that S_irr is a functionof time, depending on the evolution of the host star and migration of the planet’s orbit (Gough 1981; Lubow & Ida 2010). However, initial model runs revealed that the Sun has a very strong influence on the evolution of atmospheric CO₂ as its evolution provides a strong external forcing on the surface temperature and therefore weathering. For simplicity, we therefore fixed the solar flux at its present-day Earth value S_⊕ = 1360 W m⁻² throughout this work, while the effects of variable S_irr are considered in Appendix C. An overview of all model parameters associated with the carbon cycle model is presented inTable D.2.

Appendix C Effects ofstellar evolution on atmospheric CO₂ content

The cooling of the interior is shown in Sect. 3.1 to have important effects on the long-term behavior of the carbon reservoirs before equilibrium is established. However, another process which becomes relevant on timescales longer than t ~ 10⁹ yr is the gradual increase of the stellar luminosity S_irr (Gough 1981; Foley 2015). Therefore, we consider the effects both the thermal evolution and stellar evolution on the evolution of the atmospheric CO₂ pressure, and assess their relative importance for driving the long-term evolution of the atmospheric CO₂ pressure.

In this analysis, we restrict ourselves to Sun-like stars, such that the evolution of the stellar irradiation is given by the following parametrization derived by Gough (1981) for the Sun: (C.1)

We note however that the precise evolution of stellar irradiation over time is a function of both stellar spectral type and planetary migration. Incorporating Eqs. (C.1) into (B.13), we show the evolution of the partial atmospheric CO₂ pressure and surface temperature T_s in Fig. C.1. We here compare a case where S_irr is fixed at its present-day value with a case where S_irr is allowed to evolve over time. We set R_s0 = R_tot and R_m0 = R_k0 = 0 as initial distribution of carbon for all these models. Furthermore, we use the same parameter values and initial conditions as for the results in Fig. 6. Though we only focus on the long-term behavior, it should be noted that the converging behavior toward equilibrium for different initial distributions of carbon remains similar when the evolution of the Sun is considered.

Fig. C.1 Evolution of the atmospheric partial CO2 pressure andmean surface temperature for a case where Sirr = S⊙ and a case where Sirr evolves according to Eq. (C.1).

Figure C.1 reveals two effects of the evolution of the Sun and interior on the planetary evolution. First of all, the evolution of the Sun makes a crucial difference for the evolution of , with atmospheric CO₂ pressure dropping from values ~10⁴ Pa at t ≈ 1 Gyr to Pa at t ~ 10 Gyr. This decrease in atmospheric CO₂ results from the fact that the brightening of the Sun over time increases the surface temperature. A higher surface temperature in turn

increases the amount of continental weathering, given by Eq. (B.6), which subsequently results in depletion of atmospheric CO₂, dampening the effect of the Sun on the surface temperature. However, this dampening effect by the carbon cycle appears to be insufficient to sustain surface temperatures above the freezing point of water (T_s < 273 K) at earlier times (t ≲3 Gyr), and below the moist greenhouse threshold (T_s > 340 K) at later times (t ≳ 8.5 Gyr). The second message of Fig. C.1 is with respect to the coupled model; the gradually declining degassing due to the cooler interior may actually help in depleting the atmosphere of CO₂ to offset the gradually brightening Sun, while helping in sustaining a more CO₂-rich atmosphere at earlier times. However, it appears that this stabilizing effect is only very minor, and other mechanisms are required to mitigate the effects of the faint Sun at earlier times, and the bright Sun at later times, as is also noted by for example Sleep & Zahnle (2001).

Appendix D Parameter tables

D.1 Thermal evolution model

D.2 Carbon cycle model

References

All Tables

All Figures

	Fig. 1 Evolution of the upper mantle temperature Tm for different surface temperatures Ts and initial upper mantle temperatures Tm0.
In the text

	Fig. 2 Plate speed scaling law given in Eq. (5) as a function of the Rayleigh number Ra, along with the thickness of the upper thermal boundary layer (TBL). Values for present-day Earth are indicated for comparison.
In the text

	Fig. 3 Evolution of plate speed vp (left) and mantle temperature Tm (right) for different initial mantle and core temperatures.
In the text

	Fig. 4 A visual summary of the model considered in this study.
In the text

	Fig. 5 Initialcore temperature Tc0 (left), mantle thickness D (center), and surface gravity g (right) as a function of core mass fraction fc and planet mass Mp.
In the text

Fig. 6 Evolution of plate speed vp (A), atmosphere reservoir content Ra (B), surface temperature Ts (C), crust reservoir content Rk (D) and mantle reservoir content Rm (E) for different initial distributions of the total carbon budget Rtot over the various reservoirs. Different color depict model runs with different initial surface reservoirs Rs (t = 0) = Ra(t = 0) + Ro(t = 0). Solid lines indicate results from this work, while dotted lines represent results from Foley (2015).

In the text

	Fig. 7 Evolution of partial CO2 pressure (left) and upper mantle temperature Tm (right) for different abundances of uranium (235U and 238U) and thorium (232Th). Abundances are expressed relative to Earth (fU = fTh = 1).
In the text

	Fig. 8 Present-day partial CO2 pressure (left), surface temperature Ts (center) and mantle temperature Tm (right) as a function of thorium and uranium abundance relative to Earth.
In the text

	Fig. 9 Snapshots over time of upper mantle temperature Tm (left), the global Rayleigh number (center) and partial CO2 pressure (right) as a function of core mass fraction fc.
In the text

	Fig. 10 Snapshots over time of upper mantle temperature Tm (left), the global Rayleigh number (center) and partial CO2 pressure (right) as a function of planet mass Mp.
In the text

	Fig. 11 Equilibrium timescale τCC as a functionof mantle abundance of uranium (235U and 238U) and thorium.
In the text

	Fig. 12 Equilibrium timescale τCC as a functionof core size fc and planet mass Mp.
In the text

	Fig. 13 Equilibrium timescale τCC as a functionof weathering exponent b and arc degassing fraction f.
In the text

	Fig. 14 Long-term evolution of the upper mantle temperature Tm and partial atmospheric CO2 pressure for different values of the Nu-Ra coupling exponent β.
In the text

	Fig. A.1 Overview of the temperature profile (blue) and thermal structure considered in the thermal evolution model
In the text

	Fig. C.1 Evolution of the atmospheric partial CO2 pressure andmean surface temperature for a case where Sirr = S⊙ and a case where Sirr evolves according to Eq. (C.1).
In the text