Habitability from Tidally Induced Tectonics

A planet with a continuous and substantial source of tidal heating may have a partially molten interior.⁵ Just like the volcanism on Io, this melt may reach the surface through plumes advecting heat out to the surface. It is thus important to calculate the amount of heating a planet can experience. Based on the theory of tides (Murray & Dermott 1999) the tidal heating available to a planet depends on the eccentricity of the system e, the semimajor axis a, the mass of the star M_⋆, the mean motion n, the radius of the planet R_p, the specific dissipation parameter Q, the ratio of elastic to gravitational forces $\bar{\mu }$ ($\approx {({10}^{4}{\rm{k}}{\rm{m}}/{{R}}_{{p}})}^{2}$) and the gravitational constant G:

$$\begin{eqnarray}&&\displaystyle \frac{{dE}}{{dt}}=\displaystyle \frac{63}{4}\displaystyle \frac{{e}^{2}n}{\bar{\mu }Q}{\left(\displaystyle \frac{{R}_{p}}{a}\right)}^{5}\displaystyle \frac{{{GM}}_{\star }^{2}}{a}.\end{eqnarray} \tag{ 1 }$$

The least well constrained parameter is Q although typical values of ∼100 are commonly used for icy/rocky planets. In reality the dissipation parameter is a function of the planet's structure and should vary radially instead of being a single value. However, as a first step, we take the same approach as in most studies and use a constant value. Equation (1) can also be written as

$$\begin{eqnarray}\displaystyle \frac{{d}{E}}{{d}{t}} & \approx & 10\times {\left(\displaystyle \frac{{e}}{0.01}\right)}^{2}\left(\displaystyle \frac{100}{{Q}}\right){\left(\displaystyle \frac{{{R}}_{{p}}}{{{R}}_{{\rm{E}}}}\right)}^{7}\\ & & \times {\left(\displaystyle \frac{{{M}}_{\star }}{{{M}}_{\odot }}\right)}^{5/2}{\left(\displaystyle \frac{0.1{\rm{a}}{\rm{u}}}{{a}}\right)}^{15/2}\,[{\rm{T}}{\rm{W}}].\end{eqnarray} \tag{ 2 }$$

Therefore, an Earth-sized planet orbiting a Sun-like star at 0.1 au and with a low eccentricity of 0.01 would dissipate 10 TW or ∼0.02 W m⁻², a tidal heat flux 100 times less than Io's flux (Veeder et al. 1994; Spencer et al. 2000). For comparison, Earth's flux is estimated at 0.09 W m⁻² (Davies & Davies 2010). To reproduce the measured heat flux for Io at its location with respect to Jupiter with Equation (1), the dissipation factor would have to be chosen to be Q_Io = 150.

Given that we are interested in planets in the habitable zone that can experience tidal heating, we need to consider planets around M stars or even brown dwarfs. Planets around G-type stars that are in the habitable zone are too far away to experience any tidal heating from the star, whereas M stars have their habitable zone close enough that tidal dissipation can matter. In fact, tidal dissipation in planets around M dwarfs has previously been studied in the context of orbital and thermal evolution (Driscoll & Barnes 2015). In fact, Barnes et al. (2009) claimed that there is a limit to how much tidal heating a habitable planet may experience based on the assumption that too much volcanism and high resurfacing rates would be inhospitable. In contrast, we propose that in these conditions a new mechanism for climate regulation may kick in, rendering the planet habitable.

For a star with M_⋆ = 0.08 M_⊙, such as Trappist-1, an Earth-sized planet orbiting at 0.02 au with an eccentricity of 0.01 would experience 3000 TW or 6 W m⁻² of tidal energy. For an assumed eccentricity of 0.01 planets d and e would have tidal fluxes of 0.96 and 0.30 W m⁻², respectively. These numbers quickly grow with eccentricity, so that planets around M stars can have large tidal heating fluxes. The important aspect for the model we propose is to determine whether there is enough energy to yield a partially molten mantle. One could use Io for a simple comparison to infer whether or not planets such as Trappist-1 d, e, and f could have a partially molten mantle by comparing the amount of tidal energy available to the tidal energy carried to the surface.

According to Equation (1) tidal production increases as R_p⁵, while the heat flow carried from the mantle to the atmosphere increases as R_p², so that to zeroth order larger planets should be more likely to have an interior melt region. Even at very low eccentricity values of 0.01, Trappist-1 d would experience similar tidal heat fluxes to Io. Thus, it is possible that the Trappist-1 planets in the habitable zone have persistent melt in their interior.

This melt would be advected in the form of volcanic conduits cooling at the surface, carrying with it volatiles and degassing them into the ocean/atmosphere (see Figure 1). This continuous basalt extraction would form a layer of a certain thickness, limited by a partially molten mantle at the bottom. To estimate both the thickness of the basalt layer and the resurfacing velocity that are used in the climate model we use the model of O'Reilly & Davies (1981).

The resurfacing velocity, taken to be the same as the subsidence velocity, is

$$\begin{eqnarray}&&v=\displaystyle \frac{{q}_{t}}{\rho (L+{C}_{p}{\rm{\Delta }}T)},\end{eqnarray} \tag{ 3 }$$

where q_t is the tidal heat flux (tidal energy per unit surface area), ρ the density of the mantle, L the latent heat, C_p the heat capacity, and ΔT = T_m – T_s the difference between the melting and surface temperatures. The thickness of the basalt layer is

$$\begin{eqnarray}&&{d}_{\mathrm{bas}}=\displaystyle \frac{\kappa }{v}\mathrm{log}\left(1+\displaystyle \frac{{kv}}{\kappa }\displaystyle \frac{{\rm{\Delta }}T}{{q}_{t}\left(\tfrac{1}{a}-1\right)}\right),\end{eqnarray} \tag{ 4 }$$

where κ is the thermal conductivity, k the thermal diffusivity, and a the ratio of tidal heat flux to total heat flux. While the resurfacing rate depends only on the tidal heat flux, the thickness of the basalt layer is also determined by the contribution of tidal energy flux to total flux, a, which is unknown. To illustrate this point, we have calculated the resurfacing velocity and basalt layer thickness for different values of a. Figure 2 shows what this thickness would be as a function of contribution of tidal heat to the total heat flux. As tidal heating becomes the dominant source (i.e., a approaches 1) the thickness becomes infinitely large (as shown by O'Reilly & Davies 1981). We have also set a minimum contribution by considering a steady state where the total heat flux is from tidal heating and radioactive decay. By setting a maximum to the radioactive decay of 2 W m⁻² calculated from radioactive sources on Earth at 3.8 Gyr ago, we obtain a minimum value of a. Trappist-1 d, with a calculated tidal flux of ∼1 W m⁻² according to this simple model, would have resurfacing velocities of 6 cm yr⁻¹ and a crustal layer of at least 2 km.

Zoom In Zoom Out Reset image size

However, without a detailed calculation for the thermal history of these planets, which is beyond the scope of this paper, it is not possible to know how much of the total flux is coming from tidal heating, so as a simple first step we take v and d_bas to be set independent of each other.

Melt—Our model also requires melt to be continuously advected to the surface, which depends on how deep the melt is produced and the chemistry of the melt. On Earth, due to plate tectonics, decompression melting takes place near the surface as upper mantle material adiabatically ascends into mid-ocean ridges, and melt gets carried out as volcanism. However, if melt is produced at depth, it can stay within the mantle due to the crossover density of melt and solid at pressures of ∼13 GPa or ∼350 km on Earth (Ohtani et al. 1995). In fact, Noack et al. (2017) have argued on the basis of this crossover density that stagnant-lid planets with large masses or large core-mass fractions would not exhibit volcanism: owing to larger surface gravities, the pressure for the crossover density would be reached at shallower depths. This could pose a problem for our proposed mechanism if tidal dissipation creates melt deeper than the crossover density for the specific planet. The pressure for this crossover density depends on the composition of the solid (e.g., basalt versus peridotite, or iron content) and the water content of magma (Sakamaki et al. 2006). The higher the iron content of the solid or the drier the magma, the shallower the crossover density. Sakamaki et al. (2006) show that, for Earth, hydrated magma occurring from melting basalt rock in a mid-ocean ridge has a crossover density near 410 km (∼24 GPa) if water is present at the 2% level, and deeper at the 8% level. By extrapolating their data and using the preliminary reference Earth model (Dziewonski & Anderson 1981), we obtain a depth for crossover density for 8% water in magma of 1400 km (∼60 GPa) for Earth. Given as well that tidal strain is largest near the surface, it is reasonable to think that tidal heating may produce melt near the surface, and thus our proposed mechanism may work for a subset of planets.

In broad strokes, the carbon–silicate cycle is the process by which CO₂ is drawn from the atmosphere by removing C when it interacts with silicate rocks via water and carbonic acid, which releases bicarbonate, Ca⁺⁺ (and Mg⁺⁺) ions, either on the continents or directly with basalt on the newly formed ocean crust. If it is on the continents, these ions eventually find their way to the ocean floor, carried down by rivers and ground water. On the ocean floor they form carbonates that trap carbon (that was once in the atmosphere) into rocks, which eventually get subducted into the mantle at convergent margins. The net reaction is captured in the Urey reaction (Urey 1952):

$$\begin{eqnarray}&&{{\rm{CO}}}_{2}({\rm{g}})+{{\rm{CaSiO}}}_{3}({\rm{s}})\Longleftrightarrow {{\rm{SiO}}}_{2}({\rm{s}})+{{\rm{CaCO}}}_{3}({\rm{s}}),\end{eqnarray} \tag{ 5 }$$

where the right arrow describes how carbonates are formed. Once at depth, these carbonate rocks can react back via metamorphism, releasing carbon dioxide that eventually makes it back into the atmosphere via volcanism (the reverse reaction, left arrow).

On Earth, this cycle is enabled by plate tectonics that continuously exposes fresh rock on the ocean floor, as well as occurring on the continents driven by erosion via orogeny and the build-up of topography. Rock weathering has been well studied (Walker et al. 1981; Berner & Caldeira 1997; Berner 2004), while sea weathering has gained more interest in recent years (Coogan & Gillis 2013; Coogan & Dosso 2015). The carbonate reactions of sea weathering may be controlled in a different manner to continental weathering, with secondary carbonates and alkalinity playing an important role (Coogan & Dosso 2015). The contribution of each weathering type to the total weathering process on Earth, or how it might have changed throughout Earth's history, is still unknown (Mills et al. 2014), although recent estimates by Krissansen-Totton & Catling (2017) suggest that continental weathering has been dominant in the last 100 Myr.

For the planets proposed here, we envision that tidally driven tectonism can continuously expose fresh rock under an ocean, but that any continental shelves, if present, may not exhibit rock weathering owing to a sustained lack of erosion. In this case, the only mechanism for drawing down CO₂ in these planets would be sea weathering.

While some authors have considered sea weathering to be dependent on CO₂ concentrations alone (Sleep & Zahnle 2001; Foley 2015), laboratory experiments (Brady & Gíslason 1997) and isotopic constraints from oceanic carbonates (Coogan & Dosso 2015) suggest a temperature dependence in the form of an Arrhenius law. The temperature of the deep ocean, where Ca leaching and carbonate production are taking place, depends in turn on the atmospheric temperature, and hence on the atmospheric CO₂ concentrations, because atmospheric temperature determines how much cold surface water enters the thermohaline circulation system (Brady & Gíslason 1997). This dependence of sea weathering on atmospheric temperature and CO₂ is crucial in allowing for a climate-stabilizing feedback. We adopt the same equation as Mills et al. (2014) to describe sea weathering as a function of atmospheric CO₂ levels and atmospheric temperature T,

$$\begin{eqnarray}&&{W}_{\mathrm{sea}}=\omega {W}_{\mathrm{sea}}^{{\rm{E}}}{\left(\displaystyle \frac{{p{\mathrm{CO}}_{2}}^{\mathrm{atm}}}{{p}_{{\rm{E}}}}\right)}^{\alpha }\exp \left(\displaystyle \frac{E}{R}\left(\displaystyle \frac{1}{{T}_{{\rm{E}}}}-\displaystyle \frac{1}{T}\right)\right),\end{eqnarray} \tag{ 6 }$$

where R is the gas constant, E is the activation energy, W_sea^E is the baseline estimate of the sea weathering rate (in bar Myr⁻¹ or mol Myr⁻¹) for a reference state that we set at first to be the equilibrium atmospheric partial pressure of carbon dioxide p_E and atmospheric equilibrium temperature T_E, and ω is a function that lumps together all other quantities that sea weathering might depend on. The feedback mechanism comes from the dependence on atmospheric temperature. Any deviations from the equilibrium temperature would drive much higher or lower levels of sea weathering, depending on whether the planet is hot with high levels of atmospheric CO₂ or cold, respectively.

In addition, sea weathering may include a dependence (through ω) on the velocity at which fresh rock is exposed, which in the case of the Earth is the ridge spreading velocity, which has has changed over time. For Earth this term is often written as $\omega ={\left(\tfrac{v}{{v}_{\oplus }}\right)}^{\beta }$, where v_⊕ is the present-day velocity. While many authors (Sleep & Zahnle 2001; Mills et al. 2014; Foley 2015) use β = 1, Krissansen-Totton & Catling (2017) propose β ∼ 0.5. In either case, this contribution is of order unity. Other factors may also be affecting the weathering rate (e.g., alkalinity, grain size, etc.), and thus for simplicity we take ω ∼ 1; or conversely our results for sea weathering can be taken as scaling with ω.

Two terms in Equation (6) are poorly known: the direct dependence on atmospheric carbon dioxide pressure α and the activation energy E. Brady & Gíslason (1997) proposed α = 0.23 and E = 41 kJ mol⁻¹, while Coogan & Dosso (2015) proposed E = 92 kJ mol⁻¹, and Krissansen-Totton & Catling (2017) considered a range of E = 40–110 kJ mol⁻¹. We vary both α and E, as well as perform a stability analysis to see the effect on the climate-stabilizing feedback we propose here.

Alkalinity—Coogan & Dosso (2015) have argued that carbonate mineral precipitation is largely controlled by alkalinity, and because rock dissolution involved in sea weathering increases alkalinity the effect is to efficiently drive carbon sequestration. We note that Equation (6) does not have an explicit dependence on alkalinity, but an implicit one, because it depends on global temperatures that influence deep-water ocean temperature, which in turn controls the dissolution rates of mafic and ultramafic rocks (Krissansen-Totton & Catling 2017).

To build our model of the carbon cycle we divide our planet into three distinct reservoirs: the mantle, the basalt layer, and the ocean and atmosphere together (similar to Sleep & Zahnle 2001, Foley 2015). See Figure 1 for a cartoon representation of the tectonic process we propose in this study.

Magmatic volcanism originates in the molten mantle layer and carries melt within pipes that outgas CO₂ into the reservoirs of the ocean and atmosphere at a rate proportional to the resurfacing rate v. The flux of CO₂ degassed into the atmosphere–ocean reservoir is

$$\begin{eqnarray}&&D={p{\mathrm{CO}}_{2}}^{\mathrm{man}}\left(\displaystyle \frac{{{fvA}}_{p}}{{V}_{\mathrm{man}}}\right),\end{eqnarray} \tag{ 7 }$$

where ${p{\mathrm{CO}}_{2}}^{\mathrm{man}}$ is the partial pressure of carbon dioxide in the mantle, V_man is the volume of the mantle, A_p is the area of the planet, and f is the fraction of CO₂ within the melt that gets outgassed. On the other hand, the sink for the atmosphere–ocean reservoir is the flux of CO₂ that is weathered at the bottom of the ocean, namely the sea weathering rate described in Equation (6). By drawing out CO₂ from the ocean, sea weathering effectively draws out CO₂ from the atmosphere, given that the partitioning between the two reservoirs is set by how soluble CO₂ is. Following Foley (2015), we use Henry's law for solubility to determine how much CO₂ is in the atmosphere (${p{\mathrm{CO}}_{2}}^{\mathrm{atm}}$) versus the ocean (${p{\mathrm{CO}}_{2}}^{\mathrm{oc}}$),

$$\begin{eqnarray}&&{p{\mathrm{CO}}_{2}}^{\mathrm{atm}}={p{\mathrm{CO}}_{2}}^{\mathrm{atm}+\mathrm{oc}}-{p{\mathrm{CO}}_{2}}^{\mathrm{oc}}\\ &&\,\,=\,\displaystyle \frac{{K}_{H}}{{\mu }_{{\mathrm{CO}}_{2}}}\displaystyle \frac{{p{\mathrm{CO}}_{2}}^{\mathrm{oc}}}{\left(\tfrac{p{{\rm{H}}}_{2}{\rm{O}}}{{\mu }_{{{\rm{H}}}_{2}{\rm{O}}}}+\tfrac{{p{\mathrm{CO}}_{2}}^{\mathrm{oc}}}{{\mu }_{{\mathrm{CO}}_{2}}}\right)},\end{eqnarray} \tag{ 8 }$$

where K_H is the solubility constant, $p{{\rm{H}}}_{2}{\rm{O}}$ is the content of water in the ocean, and μ is the molar mass. The content of water in the ocean is calculated by assuming the mass of the Earth's ocean, M_oc.

While sea weathering is a sink for the atmosphere–ocean reservoir, it behaves like a source for the basaltic crust. Once C is sequestered into the basalt layer, it starts getting buried by subsequent melt deposited at the top. The C subsides until it reaches the bottom of the basalt layer, at which point it is delaminated into the molten mantle layer. To keep all the quantities in the same metrics, we use C and CO₂ interchangeably, knowing that C resides in rocks while CO₂ is the gaseous form (degassing from the mantle).

The rate at which CO₂ in the basalt is foundered or delaminated into the mantle is

$$\begin{eqnarray}&&{F}_{\mathrm{found}}={p{\mathrm{CO}}_{2}}^{\mathrm{bas}}\displaystyle \frac{v}{{d}_{\mathrm{bas}}}.\end{eqnarray} \tag{ 9 }$$

Therefore the source of carbon dioxide for the mantle is the carbon dioxide foundered from the basalt layer, and the sink is the flux being outgassed through volcanism. Thus, the equations describing this system are

$$\begin{eqnarray}\displaystyle \frac{d}{{dt}}{p{\mathrm{CO}}_{2}}^{\mathrm{atm}+\mathrm{oc}} & = & D({p{\mathrm{CO}}_{2}}^{\mathrm{man}})-{W}_{\mathrm{sea}}({p{\mathrm{CO}}_{2}}^{\mathrm{atm}},T)\\ \displaystyle \frac{d}{{dt}}{p{\mathrm{CO}}_{2}}^{\mathrm{man}} & = & {F}_{\mathrm{found}}({p{\mathrm{CO}}_{2}}^{\mathrm{bas}})-D({p{\mathrm{CO}}_{2}}^{\mathrm{man}})\\ \displaystyle \frac{d}{{dt}}{p{\mathrm{CO}}_{2}}^{\mathrm{bas}} & = & {W}_{\mathrm{sea}}({p{\mathrm{CO}}_{2}}^{\mathrm{atm}},T)-{F}_{\mathrm{found}}({p{\mathrm{CO}}_{2}}^{\mathrm{bas}}).\end{eqnarray} \tag{ 10 }$$

Given that the total content of CO₂ for the planet is fixed, ${p{\mathrm{CO}}_{2}}^{\mathrm{atm}+\mathrm{oc}}+{p{\mathrm{CO}}_{2}}^{\mathrm{man}}+{p{\mathrm{CO}}_{2}}^{\mathrm{bas}}={p{\mathrm{CO}}_{2}}^{\mathrm{tot}}$, one of the three equations is redundant.

The equilibrium carbon dioxide values for the mantle, atmosphere, and ocean are

$$\begin{eqnarray}&&{p}_{{\rm{E}}}^{\mathrm{man}}={W}_{\mathrm{sea}}^{{\rm{E}}}\displaystyle \frac{{V}_{\mathrm{man}}}{{{fvA}}_{p}}\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{p}_{{\rm{E}}}^{\mathrm{ao}}={p{\mathrm{CO}}_{2}}^{\mathrm{tot}}-{W}_{\mathrm{sea}}^{{\rm{E}}}\left(\displaystyle \frac{{d}_{\mathrm{bas}}}{v}+\displaystyle \frac{{V}_{\mathrm{man}}}{{{fvA}}_{p}}\right)\end{eqnarray} \tag{ 12 }$$

and ${p}_{{\rm{E}}}^{\mathrm{oc}}={p}_{{\rm{E}}}^{\mathrm{ao}}-{p}_{{\rm{E}}}$ in dimensional form.

For typical values refer to Table 1. To be able to solve the system of Equations (10) or (27), we need to specify the equilibrium conditions for the atmosphere, namely the equilibrium partial pressure of CO₂, p_E, that sets the equilibrium atmospheric temperature. While Menou (2015) took p_E = 330 μbar, the pre-industrialization partial pressure of carbon dioxide, as the equilibrium value for the Earth, Haqq-Misra et al. (2016) argued that a more appropriate value for an abiotic Earth would be p_E = 0.01 bar, by including all the carbon dioxide currently stored in the soil. We also take the equilibrium partial pressure to be p_E = 0.01 bar.

Table 1.  Model Parameters

Parameter	Symbol	Value	Reference
Sea weathering dependence on CO2	α	0.23	(1)
Sea weathering at equilibrium	WseaE	0.12 bar Myr−1	calculated
Equilibrium atmospheric temperature	TE	288 K	assumed
Equilibrium atmospheric CO2 partial pressure	pE	0.01 bar	(2)
Activation energy	Esea	41 kJ	nominal value, (1)
Outgassing CO2 fraction	f	0.01	nominal value, assumed
CO2 molecular weight	${\mu }_{{\mathrm{CO}}_{2}}$	44.1 mol g−1	calculated
H2O molecular weight	${\mu }_{{{\rm{H}}}_{2}{\rm{O}}}$	18.015 mol g−1	calculated
Mass of ocean	Moc	1.4 × 1024 kg
Solubility constant	Kh	235.48 bar	calculated
Radius of planet	Rp	6371 km	(3)
Radius of core	Rc	3400 km	(3)
Mass of planet	Mp	5.972 × 1024 kg	(3)
Zenith angle	μ	0.232	calculated
Sun's luminosity today		1361 W m−2
Thickness of basalt crust	dbas	10 km	nominal value, assumed
Resurfacing velocity	v	1 cm yr−1	nominal value, assumed
Total carbon content	${p{\mathrm{CO}}_{2}}^{\mathrm{tot}}$	215 bar	(4)
Characteristic timescale	τ	3 Myr	calculated

Note. References: (1) Brady & Gíslason (1997), (2) Haqq-Misra et al. (2016), (3) Stacey & Davis (2008), (4) Sleep & Zahnle (2001).

Download table as:  ASCIITypeset image

This value for the atmosphere, the solubility constant for the ocean, and the amount of water in Earth's ocean, set the equilibrium CO₂ value for the atmosphere–ocean system via Equation (8) to be ${p}_{{\rm{E}}}^{\mathrm{ao}}=0.038$ bar.

We also note that in line with astrophysical studies we use the units of bars instead of moles, and the conversion factor we use is

$$\begin{eqnarray}&&{p}{{\rm{C}}{\rm{O}}}_{2}\,[{\rm{b}}{\rm{a}}{\rm{r}}]=1.0197\times {10}^{\mathrm{-5}}\times {{\mu }}_{{{\rm{C}}{\rm{O}}}_{2}}\displaystyle \frac{{g}}{4{\pi }{{R}}_{{p}}^{2}}\times {\rm{C}}\,[{\rm{m}}{\rm{o}}{\rm{l}}],\end{eqnarray} \tag{ 13 }$$

where g is the planet's gravity and R_p is the planet's radius.

To obtain the atmospheric temperature of a planet provided the partial pressure of carbon dioxide, we use the energy balance model (EBM) of Haqq-Misra et al. (2016). They provide a fit to the outgoing longwave radiation (OLR), and the Bond albedo at the top of the atmosphere for four different type of stars including G and M stars. These quantities depend on the partial pressure of atmospheric carbon dioxide, temperature, and zenith angle μ. We simplify the model to capture global parameters by fixing the zenith angle to a value μ = 0.232 rad that would provide a global surface equilibrium temperature of T_E = 288 K for p_E = 0.01 bar for an Earth-like scenario. The zeroth-order EBM we use equates the planet's thermal radiation to the net insolation

$$\begin{eqnarray}&&{\rm{OLR}}({p{{\rm{CO}}}_{2}}^{{\rm{atm}}},T)=\displaystyle \frac{{S}_{\star }}{4}(1-A({p{{\rm{CO}}}_{2}}^{{\rm{atm}}},T)),\end{eqnarray} \tag{ 14 }$$

where S_⋆ is the solar constant at the planet's location, and A is the albedo at the top of the atmosphere. This quantity also depends on the ground albedo, which we take from Williams & Kasting (1997) and modify for M stars. For the Earth, the ground albedo is taken to be a discontinuous function of T, where for temperatures less than 273 K, the ground is considered to be frozen and exhibit a high albedo. At temperatures below 263 K the water in the atmosphere is assumed to condense out entirely, owing to the abrupt transition of the Earth entering a snowball state. However, a recent study by Checlair et al. (2017) on planets around M stars suggests that, unlike on Earth, ice coverage would proceed gradually, avoiding sudden snowball transitions, owing to the special spatial insolation pattern they receive. In our model, we use both the same albedo function for Earth and a modified version that precludes the snowball state following Checlair et al. (2017). To achieve this, we allowed for a gradual ice coverage as temperatures decrease below 273 K, retaining 30% of the land exposure at 150 K, and allowed for cloud albedo at all temperatures.

The timescale governing Equation (14) is much shorter than the long geophysical timescale involved in the carbon–silicate cycle of Equations (10) (Menou 2015). This means that in our model, surface temperature is calculated instantaneously from the amount of $p{{\mathrm{CO}}_{2}}^{\mathrm{atm}}$ at each timestep in the integration of Equations (10).

Figure 3 shows the results from the EBM model for the equilibrium case ${p{\mathrm{CO}}_{2}}^{\mathrm{atm}}=0.01$ bar. We show the case of the Earth (orange line) for comparison purposes. Earth's climate exhibits two stable points, one near 225 K where the planet would be in a snowball state, and the temperate 288 K, as well as one unstable point, near 273 K. Furthermore, the snowball state at 225 K is thought to be transient, given that rock weathering is considered to be suppressed at these temperatures so that outgassing from volcanoes eventually deglaciates the planet. It is worth mentioning that if sea weathering can operate at these low temperatures in a large enough way to balance volcanism, this snowball state could be stable in the long term. Sea weathering would have to happen on the ocean floor below a thick crust of ice, which still allows for carbon dioxide to diffuse from the atmosphere to the ocean, and volcanism would most likely have to be sluggish.

Zoom In Zoom Out Reset image size

Considering planets around M dwarfs, we also modified the EBM to restrict snowball states. This allows only one state with a stable temperature (purple line).

We investigate whether an adequate surface temperature can be kept during the evolution of an M star over billions of years after the activity of the first billion. Compared to G-type stars such as our Sun, which brighten over time after hydrogen ignition takes place, M stars get dimmer. This means that the habitable zone moves in with time, and hence planets that are found in the habitable zone were once too hot. This raises issues as to the likelihood of these planets being desiccated or wet after the intense EUV/XUV star fluxes. While earlier studies (Barnes et al. 2013) estimated planets around M dwarfs to be completely dry, recent works (Schaefer et al. 2016; Bolmont et al. 2017) that use a better prescription for the atmospheric loss (Tian 2015) suggest that some water is retained, making them prospective planets for habitability.

For the evolution of an M dwarf, we use a spline fit to the model by Baraffe et al. (2015) for a star of mass M_star =0.08 M_Sun such as Trappist-1. Figure 4 shows the evolution of Trappist-1 compared to that of our Sun.

Zoom In Zoom Out Reset image size