THE INFLUENCE OF PRESSURE-DEPENDENT VISCOSITY ON THE THERMAL EVOLUTION OF SUPER-EARTHS

For mantle viscosity contrasts larger than ∼3–4 orders of magnitude a stagnant layer forms on top of the convecting mantle (e.g., Solomatov 1995). Typical viscosity contrasts in the silicate mantles of terrestrial planets can be on the order of many magnitudes due to the interior temperature distribution, placing these planets into the so-called SL regime. Just below the SL is the upper thermal boundary layer δ_u (Figure 2). The viscosity contrast across this upper thermal boundary layer has been shown to be small, i.e., a factor of ∼10 (Davaille & Jaupart 1993; Grasset & Parmentier 1998) and we follow the approach of Dumoulin et al. (2005) assuming that the critical Rayleigh number of the upper boundary layer Ra^u_crit is approximately the critical Rayleigh number of an isoviscous system Ra^u_crit ≈ Ra_crit^iso ≈ 500 (Dumoulin et al. 2005). In that case the corresponding Rayleigh number for the upper boundary layer is defined for the lowest viscosity of the layer η(R_m), which is at the base of δ_u at the radiusR = R_m(P_m, T_m).

This defines the upper thermal boundary layer as a function of the critical Rayleigh number for isoviscous convection, i.e., Ra(δ_u) = Ra^u_crit ≈ Ra^iso_crit for a temperature- and pressure-dependent viscosity η(P, T), compare Equation (9)

$$\begin{equation} \delta _u = \left({\frac{{\kappa _m \eta ({R_m } ){\rm Ra}_{{\rm crit}}^u }}{{\alpha _m \rho _m g\left| {T_m - T_l } \right|}}} \right)^{1/3}. \end{equation} \tag{ 15 }$$

The heat flux out of the mantle q_m,out is then defined as the conductive heat flux through the upper thermal boundary layer:

$$\begin{equation} q_{m,{\rm out}} = k_m \frac{{T_m - T_l }}{{\delta _u }}.\end{equation} \tag{ 16 }$$

In the case of PT, we assume that T_l is equal to the surface temperature T_s (e.g., Breuer & Spohn 2003). In the case of SL convection, we need to specify the temperature T_l at the base of the SL. As stated above, the viscosity contrast through the upper thermal boundary layer is approximately one order of magnitude (note that this limitation in the viscosity contrast for the upper thermal boundary layer is not required when modeling PT). Hence, the temperature T_l defining the base of the SL can be calculated as a function of the upper mantle temperature T_m and the viscosity law (Equation (8)):

$$\begin{eqnarray} T_l \approx \frac{{T_m \cdot ({E^* + P_l V_{{\rm eff}}^* } )}}{{E^* + P_m V_{{\rm eff}}^* + R_g T_m \ln 10}}&\quad {\rm with}\quad P_m = g \cdot \rho _m \cdot ({L + \delta _u } )\nonumber\\ &\quad\ \ {\rm and}\quad P_l = P_m - g\rho _m \delta _u.\nonumber\\ \end{eqnarray} \tag{ 17 }$$

For a solely temperature-dependent viscosity and R_gT_mln 10 ≪ E*, Equation (17) can be transformed into

$$\begin{equation} T_l \approx T_m - \frac{{\ln 10 \cdot R_g T_m^2 }}{{E^* }}, \end{equation} \tag{ 18 }$$

which excellently approximates the experimental results for T_l by, e.g., Davaille & Jaupart (1993). We then calculate the evolution of the lithospheric thickness L(t) via the energy balance at the base of the SL neglecting volcanic heat transfer (Schubert et al. 1979; Spohn 1991):

$$\begin{eqnarray} \rho _m c_m (T_m - T_l)\frac{{\it dL}}{{\it dt}} &=& {-} q_{m,{\rm out}} - k_m \left({\frac{{\partial T}}{{\partial R}}} \right)_{R = R_l } \nonumber\\ &&\quad {\rm with}\quad T_l = T_l ({T_m } ). \end{eqnarray} \tag{ 19 }$$

To solve this SL growth equation (Equation (19)) we need to determine the upper mantle temperature T_m(t) and the thermal gradient at the SL base $({\partial T}/ {\partial R})_{R = R_l }$. The thermal gradient at the base of the SL can be calculated by integrating the time-independent heat conduction equation:

$$\begin{equation} \frac{1}{{R^2 }}\frac{\partial }{{\partial R}}\left({R^2 k_m \frac{{\partial T}}{{\partial R}}} \right) = - Q_m \end{equation} \tag{ 20 }$$

with the SL base temperature T_l and the surface temperature T_s imposing a determined boundary problem. Note that for simplicity we use the steady-state heat solution because the typical relaxation time in the lithosphere is of the order of τ = L²/κ_m ≈ 100 Myr. This is much smaller than the relaxation time of the mantle, which is of the order of a few Gyr.

To calculate the upper mantle temperature T_m(t) we need to determine the heat flux into the convecting mantle q_m,in. For this we have to specify the lower mantle energy budget—here pressure effects are crucial and we suggest a parameterization that considers a potential stagnant zone in the lower mantle.

The heat fluxes out of the core q_c and into the convective mantle q_{m, in} are described via the thermal gradients at the top and the base of the total lower mantle conductive zone $\tilde \delta _c$, respectively (see Figure 2),

$$\begin{equation} q_{m,{\rm in}} = - k_m \left({\frac{{\partial T}}{{\partial R}}} \right)_{R = R_b } \end{equation} \tag{ 21 }$$

$$\begin{equation} q_c = - k_m \left({\frac{{\partial T}}{{\partial R}}} \right)_{R = R_c }.\end{equation} \tag{ 22 }$$

These heat fluxes can be determined from the conductive thermal profile in $\tilde \delta _c$. As the stagnant CMB-lid and $\tilde \delta _c$ can comprise a substantial part of the mantle, we have to consider the radiogenic heat sources Q_m and the time-dependent heat equation:

$$\begin{equation} \frac{1}{{R^2 }}\frac{\partial }{{\partial R}}\left({R^2 k_m \frac{{\partial T}}{{\partial R}}} \right) = - Q_m + \rho _m C_m \frac{{\partial T}}{{\partial t}}. \end{equation} \tag{ 23 }$$

To solve for Equations (21)–(23), we need to determine the thickness of $\tilde \delta _c$. The conductive lower mantle layer $\tilde \delta _c$ is assumed to be at the verge of convective stability. Thus, we can compute the thickness of $\tilde \delta _c$ similarly to the upper thermal boundary layer, by specifying the critical Rayleigh number Ra^c_crit of $\tilde \delta _c$ based on a characteristic viscosity η_c and the temperature contrast |T_c − T_b| of this layer:

$$\begin{equation} \tilde \delta _c = \left({\frac{{\kappa _m \eta _c Ra_{{\rm crit}}^c }}{{\alpha _m \rho _m g| {T_c - T_b } |}}} \right)^{1/3}. \end{equation} \tag{ 24 }$$

The critical Rayleigh number of such a layer depends on the viscosity contrast Δη through it. This has originally been derived for pressure-dependent viscosities, which increase with depth, by Schubert et al. (1969), and for solely temperature-dependent viscosities, which decrease with depth, by Stengel et al. (1982) (see as well, e.g., Dumoulin et al. 2005 and Solomatov 1995). Based on these findings, we approximate the critical Rayleigh number of a layer with a viscosity contrast Δη by

$$\begin{equation} {\rm Ra}_{{\rm crit}} ({\Delta \eta } ) \approx \max \big\{ {{\rm Ra}_{{\rm crit}}^{{\rm iso}},\phi \cdot [ {\ln ({\Delta \eta } )} ]^4 } \big\},\end{equation} \tag{ 25 }$$

where Ra^iso_crit ≈ 500 is the critical Rayleigh number for an isoviscous system and ϕ is a constant.

This equation holds when the Rayleigh number is defined for the lowest viscosity of the system and when the viscosity monotonically increases or decreases (Figure 3). The value of ϕ varies approximately between ∼10 and ∼60 depending on the boundary conditions, the viscosity law, and the amount of internal heating (e.g., Dumoulin et al. 2005; Schubert et al. 1969; Solomatov 1995; Turcotte & Schubert 2002). We choose for our calculations ϕ equal to 11.74 according to Dumoulin et al. (2005)—additional tests have shown that a variation in ϕ between 10 and 60 does not change our findings. The value ϕ · [ln (Δη)]⁴ corresponds to the asymptotic limit of the critical Rayleigh number for large viscosity contrasts, but is generally used for viscosity contrasts above ∼3–4 orders of magnitude (e.g., Solomatov 1995; Fu et al. 2010). We note that Fu et al. (2010) use for smaller viscosity contrasts of (Δη ⩽ 10³)a critical Rayleigh number of ∼1000–2000 but define the reference Rayleigh number at the average viscosity of the system. Due to the larger reference viscosity for the Rayleigh number we find this approach to be approximately consistent with our method.

Zoom In Zoom Out Reset image size

We can define the lower mantle conductive zone $\tilde \delta _c$ using Equations (24) and (25), setting the critical Rayleigh number equal to Ra^c_crit = Ra_crit(Δη) with Δη = max (η(R_c)/η(R_b), η(R_b)/η(R_c)) and the characteristic viscosity η_c = min {η(R_c), η(R_b)}:

$$\begin{equation} \tilde \delta _c = \left({{\rm Ra}_{{\rm crit}} ({\Delta \eta } ) \cdot \frac{{\kappa _m \cdot \min {\{ {\eta ({R_c } ),\eta ({R_b } )} } \}}}{{\alpha _m \rho _m g| {T_c - T_b } |}}} \right)^{1/3}. \end{equation} \tag{ 26 }$$

For V*_eff = 0 the minimal and hence characteristic viscosity is η(R_c) at the CMB (for illustration see Figures 2 and 3). In this case we find that, although the initial temperature contrast and thus the related viscosity contrast across $\tilde \delta _c$ can be large, Δη quickly (<500 Myr) decreases below one order of magnitude due to the planet's efficient heat transport, resulting in Ra_crit ≈ Ra^iso_crit and $\tilde \delta _c = \delta _c$ just as expected for a standard boundary layer model without pressure dependence. This simplifies for V*_eff = 0 the lower mantle heat flux to $q_{m,{\rm in}} = q_c = k_m ({{T_c - T_b }})/{{\tilde \delta _c }} = k_m ({{T_c - T_b }})/{{\delta _c }}$.

When pressure effects lead to a viscosity increase through $\tilde \delta _c$ with depth, as shown in Figure 2, our approach in determining $\tilde \delta _c$ is analogous to the determination of the upper conductive zone, consisting of the SL and the upper thermal boundary layer δ_u as suggested by Solomatov (1995). In this strongly pressure-dependent viscosity case, the minimal and therefore characteristic viscosity is η(R_b) at the top of $\tilde \delta _c$.

5.1.1. Reference Model

The reference model is based on the initial reference temperature profiles specified in Table 3 and explained in Section 4. The thermal evolution is shown in Figures 4–7: non-pressure-dependent SL planets (model A with V*_eff = 0) show a period of 2–5 Gyr in which the mantle and the core heat up, to strongly cool afterward. The initial phase of core heating is a consequence of the SL mode and V*_eff = 0 (for V*_eff = 0 PT planets will continuously cool; compare with Figure 12). The upper lid reduces the heat transport out of the mantle; instead heat is transported efficiently into the core. As a consequence of decreasing mantle heat sources with time the core finally starts to cool down after a few Gyr. We find that the initial phase of core heating increases with increasing planetary mass for V*_eff = 0. This is due to the relatively low initial CMB temperatures as well as the increasing amount of heat sources in the mantle with increasing planetary mass. For this pressure-independent case we find—as expected—thin lower thermal boundary layers in the order of 1–10 km, which generally decrease with increasing planetary mass. The decrease in the thickness of the lower thermal boundary layer with increasing mass is a consequence of the higher lower mantle temperatures (Figure 4) and the associated lower viscosities.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For the pressure-dependent models B and C a different thermal evolution can be observed depending on the planet's mass: for the Earth model (Figure 4, M = 1) pressure effects are relatively small, and we still observe a typical lower thermal boundary layer, which grows with the planet's cooling from 30 to 100 km in 13 Gyr. However, the effects of pressure on the viscosity are strong enough to reduce the convective vigor in the mantle. This is expressed by generally hotter interior temperatures in the mantle and core in comparison to model A (compare Figures 4 and 7 for M = 1). Nevertheless, independent of the effective activation volume V*_eff, the M = 1 SL model predominantly cools with time (Figures 4 and 7).

With increasing planetary mass the cooling behavior differs significantly between the pressure-dependent and the pressure-independent models (Figures 4–7 for M = 5 and M = 10): we find a continuous, but slow heating of the core for both pressure-dependent models. The change in temperature is less than ∼50–150 K in 13 Gyr, leading to a negligible decrease in the CMB viscosity for models B and C with time (Figure 6). The heating of the core is the consequence of the conductive heating of the lower layer by the radiogenic heat sources. The core represents a large thermal reservoir without internal heat sources; hence the temperature in the lower mantle increases faster than in the core. Interestingly, for our pressure-dependent viscosity models the core does not heat up as fast as it does for V*_eff = 0 in the first 2–5 Gyr. This can be explained by the observation that for the pressure-dependent models a large fraction of the lower mantle is conductively heated, distributing the energy produced by radiogenic heat sources over the whole lower mantle and the core. For V*_eff = 0 the heat transport is more efficient due to vigorous convection and the energy is transported directly into two distinct thermal sinks, the upper mantle and the core, keeping the average lower mantle cooler than is the case for the pressure-dependent viscosity models (compare in Figure 6 the thermal profiles for V*_eff = 0 with V*_eff > 0).

The conductive zones in the lower mantle $\tilde \delta _c$ occupy for the studied super-Earths and pressure-dependent models a large fraction of the planetary mantle. These conductive zones decrease with the initial phase of heating of the upper mantle in the first ∼3 Gyr (see Figure 7), and later grow with time as the upper mantle cools. The variation is, however, less than 15%, resulting in approximately constant $\tilde \delta _c$ thicknesses. $\tilde \delta _c$ increases in size with planetary mass and effective activation volume for pressure-dependent models (Figure 5 shows the thicknesses of $\tilde \delta _c$ averaged over 13 Gyr). The actual values are ∼4000 km and ∼5000 km for the M = 10 planet, and ∼700 km and ∼3000 km for M = 5 for models C and B, respectively.

The thermal evolution of the core for the pressure-dependent models B and C of the M = 10 planet is mainly independent of the convective mantle. A reduction of the effective activation volume, and hence an increase of convective vigor in the mantle from model B to C does not reduce the CMB and core heating, although the thickness of the lower mantle conductive zone is reduced by about 20%–40% (Figure 6). This effect is also attributed to the large thermal diffusion times in the thick conductive zones (see also discussion). Hence, the strong pressure-dependent viscosity thermally decouples the lowermost mantle and core from the convective mantle. This is different for the M = 5 planet, where pressure effects are weaker and a change from model B to C leads to a less efficient core heating by ∼30% (Figure 4).

Figure 6 shows the depth profiles of temperature and viscosity for an M = 10 planet with the three different viscosity models at three different times. Note that the adiabatic gradient depends on the upper mantle temperature T_m (e.g., see Stamenković et al. 2011), unlike in the classic linear approximation for smaller planets (e.g., Grott & Breuer 2008), and varies with time due to a change of T_m. Please further note that the innermost slope of the conductive thermal profile for models B and C (Figure 6) changes slowly with time due to the small thermal diffusion times. For larger thermal diffusivities (i.e., due to an increase of the thermal conductivity; see Section 6.4) this innermost slope diverges faster from the initial thermal slope.

For V*_eff = 0, the initial phase of heating of the upper mantle is less efficient with increasing planetary mass, but at later times (>5 Gyr) more massive planets show higher upper mantle temperatures (Figure 7). This behavior can be explained by the initially stronger heating of the core for more massive planets as previously discussed (Figure 4)—less energy is available to heat the upper mantle at times before ∼3 Gyr. At later times (>3 Gyr) the super-heated core starts to cool and releases the stored energy. Due to the earlier stronger heating of the core with increasing mass more energy is available for larger planets to be released into the upper mantle.

The thermal evolution of the upper mantle is more complex when pressure effects are taken into account (Figure 7): for the pressure-dependent models B and C, the peak upper mantle temperatures still decrease with planetary mass. However, the general thermal evolution of the upper mantle is modulated by various factors. This complexity can be demonstrated for a specific planet with a fixed mass M and fixed initial temperature distribution considering an increase in the effective activation volume (i.e., comparing the viscosity model C with the viscosity model B for M = 10): the lower conductive zone in case of model B is thicker than in model C, and hence the effectively convecting mantle with the approximate thickness D_eff ≈ R_l − R_b is smaller in comparison to model C. The thermal evolution of the upper mantle then crucially depends on (see Equation (12) and Figure 2):

• 1.  
  the effective volume Ω_eff of the convective mantle, which approximately scales as ∼D_eff · A, with A being the average surface in the convective mantle,
• 2.  
  the energy flux due to radiogenic heat sources in the fully convective region D_eff, which scales as ∼Ω_eff · Q_m · A⁻¹∼ D_eff · Q_m ∼ D_eff,
• 3.  
  the heat flux q_m,in into this fully convecting zone out of $\tilde \delta _c$, and
• 4.  
  the heat flux out of the mantle q_m,out,

which all strongly depend on the effective activation volume and the time-dependent ratio of the effective mantle thickness to the thickness of the conductive lower mantle zone $r_{{\rm eff}} (t) = {D_{{\rm eff}} (t )}/ {\tilde \delta _c (t )}$. The ratio r_eff decreases with increasing effective activation volume for a fixed planetary mass with the same initial thermal profile: in this case the heat flux into the convecting part of mantle q_m,in generally increases due to a stronger heating of the conductive zone in the lower mantle $\tilde \delta _c$, and the heat flux out of the mantle is reduced due to an increase in the thickness of the upper thermal boundary layer (indicating the reduced convective vigor in the upper mantle). This implies that for a specific planetary mass the upper mantle temperature T_m(t) should grow when the effective activation volume is increased. This is observed for intermediate pressure effects and smaller planets like our model Earth (Figure 7) when r_eff is still large.

For more massive planets additional effects complicate the thermal evolution of the upper mantle when the effective activation volume is further increased (comparing the viscosity model C with the viscosity model B for M = 5, 10 in Figure 7): r_eff strongly decreases, reducing the heat contribution due to radiogenic heat sources in Ω_eff. This leads generally to cooler upper mantle temperatures for model B in comparison to model C. However, for M = 10 we observe after ∼7 Gyr (Figure 7) that the upper mantle temperature is again larger for the viscosity model B in comparison to the viscosity model C. This can be explained with the growth of the lower mantle conductive zone as well as the strong increase of the upper mantle SL thickness (L increases by a factor of three in 13 Gyr) with time, leading for the viscosity model B and M = 10 to a thin convective mantle (∼500 km). The effectively convecting mantle volume becomes so small that all incoming energy is directly converted into solely increasing T_m. This complex behavior also indicates that a simple scaling from planetary mass to the total heat flux is difficult when pressure effects of the viscosity are considered.

5.1.2. Variation of Initial Mantle and Core Temperatures

We first keep the reference initial CMB temperatures but increase the initial upper mantle temperature to the maximal peridotite solidus temperature in this zone, i.e., T_m(0) = 2300 K (see Section 4, Table 3): for all viscosity models the initially hotter upper mantle cools quickly and converges to the same upper mantle temperatures as for the initially cooler reference model—after about 2 Gyr the difference in upper mantle temperatures is less than 5% between the two models. This leads in the first ∼2 Gyr to a smaller conductive zone in the lower mantle by less than 25% for the initially hotter upper mantle model.

The impact of a varied initial upper mantle temperature on the thermal evolution of the core depends on the viscosity model and the planetary mass: the larger the lower mantle conductive zone (hence the more massive the planet or the larger the effective activation volume) the smaller is the impact from an initially hotter upper mantle on the lower mantle and the core. For both pressure-dependent models the initially hotter upper mantle heats the core additionally by ∼5%–30% during the entire evolution. For V*_eff = 0 the thermal exchange between the upper mantle and the core is stronger, and we observe a stronger heating of the core during the first 2–5 Gyr. The maximal core heating is enlarged by about 40%–100%. However, due to the large CMB temperatures above 5000–6000 K this difference in core heating for all three viscosity models corresponds only to a small variation of the core temperature (<5%), which illustrates that the planet's dynamics is not strongly affected by the increase in initial upper mantle temperature.

The situation is different when varying the initial CMB temperature: with increasing initial CMB temperature the thickness of the lower mantle conductive zone $\tilde \delta _c$ strongly decreases (Figure 8). For both super-Earth models (M = 5, 10) conductive zones in the lower mantle thicker than 500 km exist for initial CMB temperatures below ∼10,000–17,000 K for model B, and below ∼6000–12,000 K for model C, respectively. The decrease of the thickness of the lower mantle conductive zone with increasing initial CMB temperature reduces the efficiency of core heating and allows the core to start cooling above a certain initial CMB temperature (Figure 9). We find that the M = 10 planet starts to cool at initial CMB temperatures above ∼12,000 K for the viscosity model B and at initial CMB temperatures above 9000 K for model C. For the M = 5 planet and model B cooling starts at initial CMB temperatures of ∼7500 K; and for model C at initial CMB temperatures above ∼5600 K.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 9 shows the time evolution of the CMB temperature for the M = 5 planet and the viscosity models A and C with varying initial CMB temperatures up to the melting temperature of MgSiO₃ perovskite at ∼570 GPa (∼13,500 K; Stamenković et al. 2011). For a pressure-dependent viscosity the cooling is much less efficient in comparison to V*_eff = 0 even for high initial CMB temperatures. It is important to note that for the pressure-dependent models the CMB temperatures can differ by several thousand degrees kelvin after 13 Gyr (about 3000 K for the M = 5 planet in Figure 9) depending on the initial CMB temperature—the difference is much smaller without pressure dependence of the viscosity.

We examine the thermal evolution of planets with varying mass assuming the planets to be in the PT instead of the SL regime. PT increases the cooling rate of the upper mantle in contrast to SL convection, which results for pressure-independent rheologies also into an effective cooling of the core (Figures 10 and 11 for V*_eff = 0; compare with, e.g., Breuer & Moore 2007). This is in contrast to SL convection where the deep interior remains comparatively warm as already suggested for smaller terrestrial planets (e.g., Breuer & Moore 2007).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

With increasing V*_eff and increasing planetary mass, however, we observe that the cooling behavior of the core and lowermost mantle becomes increasingly similar between the SL and the PT regime. This might be counterintuitive at first, but it reflects the ineffective convection of the lower mantle with increasing V*_eff for both convection regimes. For super-Earths where CMB-lids are present, the core is even insulated by a conductive zone and hence not directly influenced by the convection in the upper mantle—and thus not affected by a change in the convective regime. This is shown in Figure 10 for the M = 10 planet with V*_eff = 1.7 cm³ mol⁻¹ and T_c(0) = 6100 K: in contrast to SL convection, PT leads to a stronger cooling of the upper mantle by about 600 K in 13 Gyr, and thus to an increase of the viscosity above the lower conductive zone $\tilde \delta _c$. As a consequence $\tilde \delta _c$ increases by about 10%, leading to an even better thermal insulation of the core. We observe this similarity in lower mantle and core cooling between the SL and the PT regimes for all studied super-Earths with viscosity models B and C, independent of initial CMB temperature (up to the melting temperature of MgSiO₃ pv at the CMB). For smaller planets (below M = 1) or super-Earths with smaller effective activation volumes than suggested by model C the effects of pressure on the viscosity are smaller and we observe that PT can again more efficiently cool the core in comparison to SL convection.

Interestingly, with increasing effective activation volume the cooling of the core can become similar between an M = 1 planet in an SL and an M = 1 planet in a PT regime: for V*_eff = 0 we observe a typical behavior where PT can much more effectively cool the core in comparison to SL convection (Figure 11). Independent of initial CMB conditions we obtain for V*_eff = 0 and a PT Earth model a present-day CMB temperature of ∼2450 K, about 1300 K below the current estimate (see Kawai & Tsuchiya 2009; Stacey & Davis 2004; Tateno et al. 2009). The additional core cooling due to PT in comparison to SL convection is for V*_eff = 0 about 760 K in 13 Gyr, for V*_eff = 1.6 cm³ mol⁻¹ it is 380 K, and for V*_eff = 2.5 cm³ mol⁻¹ the core cooling due to PT approaches the one caused by SL convection.

Note that for V*_eff = 2.5 cm³ mol⁻¹, i.e., the assumed effective activation volume at Earth's CMB for MgSiO₃ perovskite (Stamenković et al. 2011), and a PT mode, the present-day CMB temperature is, for initial CMB temperatures up to 10,000 K, between 3550 and 3800 K. This value is similar to the current CMB temperature estimate of the Earth (see Kawai & Tsuchiya 2009; Stacey & Davis 2004; Tateno et al. 2009) in contrast to V*_eff = 0 and active PT (∼2450 K). As a comparison: for V*_eff = 0 our Earth model could only reach higher CMB temperatures if we model the Earth throughout 4.5 Gyr in the SL mode, but even then we obtain maximal CMB temperatures between 3300 and 3400 K today. This of course neglects the possibility of mantle layering as discussed later.

The 1D model is a first approach to modify classic parameterized 1D convection models by including strong pressure effects on the viscosity. In the following, we compare the 1D results with those of a 2D spherical convection model to examine whether our approach can describe satisfactorily the thermal evolution with temperature- and pressure-dependent viscosity. It is not aimed to derive or adjust the parameterization from the 2D results as this requires a substantial study far beyond the scope of this paper.

We have performed a variety of 2D spherical convection simulations assuming the Boussinesq approximation and compare them with our 1D parameterized calculations. We use the same parameters as in the reference model for the 1D runs (Tables 1–4) apart from the reference viscosity and the initial temperature profile to sustain numerical stability. Details of the 2D convection model including the chosen parameters are described in Appendix B and in Table 5.

To distinguish the convective and stagnant regions in the 2D model, we use the Peclet number, Pe, which is the ratio between the energy transported by convection and by diffusion, Pe = uD/κ_m with u being the velocity, D the mantle thickness, and κ_m the thermal diffusivity. For Peclet numbers smaller than one, convection velocities are almost zero and we define the region to be stagnant.

We further distinguish between transitional and convective regions. For Peclet numbers between 1 and 100 the convection is highly sluggish (depending on planetary mass Pe = 100 corresponds to mantle velocities below 0.1–1 mm yr⁻¹). The effectively convecting region is defined for Pe numbers larger than 100. The laterally averaged values for the CMB-lid and transitional lid thicknesses have to be considered with care as information about local partial stagnant zones might be lost, or as very low convective velocities might be associated with a zone where, based on our definition, no stagnant (CMB) lid exists.

All 2D convection models confirm the formation of a stagnant lower mantle zone, a CMB-lid, during the thermal evolution of super-Earths with 5 or 10 Earth masses for the reference temperatures. The CMB-lid thickness is found to increase strongly with planetary mass, effective activation volume, and reference viscosity, and to decrease for larger interior temperatures as suggested already by the 1D models. For our Earth model (M = 1) stagnant zones at the CMB are not present but small transitional zones are observed assuming a relatively high reference viscosity of 10²⁵ Pa s (Table 5; note that the transitional zones would decrease if we assume a more appropriate reference viscosity of 10²¹ Pa s).

In the following, we show the most significant 2D results for M = 5 and M = 10 planets at times when the convective vigor peaks as well as after 13 Gyr of thermal evolution to illustrate the general convective behavior assuming pressure-dependent viscosity in the mantles of super-Earths (Figures 13 and 14). Note that for convection to start it takes up to 1–3 Gyr depending on initial thermal profile and reference viscosity.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 13 shows an M = 10 planet with standard initial temperatures and a reference viscosity of 10²¹ Pa s. The convective velocities reach a maximum of ∼3.5 cm yr⁻¹ at ∼7.5 Gyr and decrease to maximal ∼0.5 cm yr⁻¹ after 13 Gyr of evolution. The CMB-lid comprises more than ∼60% of the entire mantle and the CMB heats slowly by ∼50K during 13 Gyr, indicating the conductive heating of the lower mantle and core as observed and discussed for the 1D model. Although the uppermost mantle convects, we find for all studied M = 10 cases, independent of the chosen initial upper mantle or CMB temperature (see Appendix B), that a transitional zone comprises almost the whole convecting mantle, indicating a generally sluggish convection for the M = 10 planet. This is indicated by the rather low convective velocities thoughout the planet's mantle (Figures 13(b) and (e)).

The M = 5 planet shows more vigorous convection in comparison to the M = 10 planet. Figure 14 shows results for a 5 Earth mass planet and V*_eff = 2.5 cm³ mol⁻¹, T_c(0) = 5100 K, and η_ref = 10²² Pa s at t = 7.5 Gyr, when the convection is strongest. For this case we find thinner CMB-lids, which comprise about 40% of the entire mantle, and a transitional zone that comprises about 70% of the non-stagnant mantle.

Our 2D results show that for the suggested temperature- and pressure-dependent rheologies and the standard temperature profiles convection strength in super-Earths decreases with planetary mass, leading to a stagnant lower mantle and to a highly sluggish convection above large CMB-lids (Table 5).

As described in Appendix B, we have performed additional runs with CMB temperatures increased by 1000–2000 K, as well as runs with an initially hotter upper and/or hotter average mantle. But even here we generally find a strong reduction of convective vigor for more massive planets as well as increasing CMB-lids with planetary mass. For those initially hotter planet models we specifically find the following results, which are all in excellent agreement with our 1D model:

• 1.  
  An increase of CMB temperature (2000 K for M = 10, 1000 K for M = 5) does not significantly affect the thickness of the CMB-lid or the transitional lid.
• 2.  
  An initially hotter upper mantle (by reducing the initial upper thermal boundary layer thickness, see Appendix B) does not significantly affect the CMB-lids, although it helps to earlier ignite upper mantle convection.
• 3.  
  An initially hotter average mantle (by increasing the initial intermediate mantle temperature with planetary mass, see Appendix B) contributes in reducing the CMB-lid thicknesses, but still results in large transitional lids, which comprise a significant fraction of the convective mantle.

The results of the 2D and 1D models are in good agreement. In particular, the thermal profiles in the lower mantle and the CMB temperatures are very similar between the two models for M = 5 and M = 10: for instance, the difference in core heating assuming the reference temperatures is less than 10% and the CMB temperatures differ by no more than 2% at all times during 13 Gyr.

The thermal profiles, however, differ in the upper part of the lower mantle conductive zone and the uppermost convective mantle (Figures 15(a) and (c)). This difference can be explained by the assumption made for the 1D model that the zone above $\tilde \delta _c$ is adiabatic (isothermal) between the thermal boundary layers and by not separating the lower thermal boundary layer and the CMB-lid. The first assumption, i.e., the convective temperature profile in our 1D model being adiabatic, implies a vigorously convecting mantle layer. Note that we use the Boussinesq approximation in the 2D model, which neglects adiabatic effects and viscous heating, and results in an almost isothermal convecting zone when the convection is vigorous. A comparison between the 2D and the 1D model, which intrinsically accounts for an adiabatic gradient, indicates that the Boussinesq approximation in the 2D model underestimates lower mantle temperatures. This temperature underestimation is not significant, indicating that the Boussinesq approximation can be used as a first step (e.g., O'Neill & Lenardic 2007; Stein et al. 2011). A more detailed analysis is required to account for the adiabatic and viscous effects in the 2D model, but is beyond the scope of this work. For a direct comparison with the Boussinesq approximated 2D model, we adjust our corresponding 1D model calculations by assuming that the adiabatic temperature increase in the convecting zone is negligible, and thus is isothermal. The 2D thermal profile shows an isothermal upper mantle at the time of strongest convection (Figure 15(a)) but later in time the temperature increases with depth across the upper mantle due to sluggish convection (Figure 15(c)). Thus, the 1D model overestimates the convective vigor in the convecting zone (above $\tilde \delta _c$) as we always assume strong, non-sluggish convection in this layer.

Zoom In Zoom Out Reset image size

The conductive behavior in the lower mantle is on the other hand overestimated in the 1D model because we cannot separate the lower thermal boundary layer and the CMB-lid. This leads to a smaller cooling of the lower mantle in the 1D model compared to the 2D model but the difference is comparatively small. In addition, a direct comparison between the CMB-lid based on the 2D model and the 1D lower thermal conductive zone $\tilde \delta _c$ is difficult. Nevertheless, we find that $\tilde \delta _c$ is always in between the 2D CMB-lid and transitional lid, increases with planetary mass, effective activation volume, and reference viscosities, and decreases with interior temperatures in a similar manner as observed with the 2D model.

Although the differences in the results between the 1D and 2D models can be relevant for a detailed study of the thermal evolution of super-Earths, we find that the general thermal trends, in particular of the lower mantle and the core, are represented satisfactorily with the 1D model. Being aware of the 1D model's limitations (see also below), we can use it to discuss general consequences of a pressure-dependent viscosity on the thermal evolution of super-Earths.

In the following, we want to address some limitations of the 1D model and their possible impact on our findings. The critical Rayleigh number based on the viscosity contrast (Equation (25)) was derived for monotonically increasing or decreasing viscosities (Figure 3). For some cases, however, we find that due to radiogenic heating in the conductive zone of the lower mantle $\tilde \delta _c$ the viscosity profile deviates with time from a monotonic behavior (Figure 6). This deviation shifts the actual minimum of the viscosity profile to greater depths below R_b and leads to an effectively larger viscosity contrast of this subsystem. We performed all runs adjusting for this viscosity shift with time to greater depths, by using for the calculation of $\tilde \delta _c$ (Equation (26)) the minimal viscosity of the layer as the characteristic viscosity of the system instead of the viscosity at R_b as well as the maximal viscosity contrast of the layer to calculate Ra_crit(Δη). We find that the increase of the viscosity contrast due to this shift compensates the effects of the reduced characteristic viscosity of the system. Thus, the thermal evolution does not significantly differ (deviation in T_m and T_c is less than 1%), although $\tilde \delta _c$ is reduced by ∼10%–50% and shows a generally stronger decrease with time. We hence conclude that this deviation is not large enough to change the thermal trends observed by our model. In addition, we have to be aware that classic Rayleigh numbers are defined for a linear temperature increase, which is due to conductive heating of the lower mantle not any longer given (see Figure 6). However, we find this deviation to be small, and observe no essential change of the thermal evolution if we use the maximal temperature contrast of the system instead of ΔT = T_c − T_b. Note that the deviation from a linear temperature increase is smaller considering also the depth dependence of the thermal conductivity. In this case, the lower mantle cools more effectively and the thermal profile in the lower mantle conductive zone becomes increasingly linear, leading to a monotonically increasing or decreasing viscosity (see Section 6.4).

The results of the 2D model with a thick CMB-lid (e.g., Figures 13–15) indicate that convection in the upper mantle is sluggish in contrast to what we assume in the 1D model. An overestimation of the convective strength is also possible when the conductive layer at the CMB is thin due to high CMB temperatures. This scenario is illustrated in Figure 16 and occurs for the pressure-dependent models B and C at initial CMB temperatures above ∼8000 K and 9000–11,000 K for the M = 5 and the M = 10 planets, respectively. For these cases, the viscosity contrasts in the convective adiabatic layer above $\tilde \delta _c$ are between 10⁵ and 10¹⁰ and the absolute viscosities can reach values between 10²⁵ Pa s and 10³⁵ Pa s—conditions not favorable for convection in the zone above $\tilde \delta _c$ (the red zone in Figure 16 in the online journal).

Zoom In Zoom Out Reset image size

This, on the other hand, suggests that the cooling of the interior might be less efficient than predicted by our model. In addition, in this super-heated core scenario (Figure 16) the stability of the conductive layer $\tilde \delta _c$ at the CMB is difficult to predict: the viscosity decreases with pressure across the conductive layer and the viscosity contrasts can reach ∼10 orders of magnitude, but the least stable part of $\tilde \delta _c$ is at its bottom where the viscosity is minimal. This suggests that the uppermost part of $\tilde \delta _c$ is stagnant and that the lowermost part of $\tilde \delta _c$ may be unstable.

In the present 1D model, we do not distinguish between the lower thermal boundary layer and the CMB-lid. This introduces an error into the energy conservation equation of the mantle (Equation (12)) as we underestimate the volume and hence the amount of radiogenic heat sources in the convective zone. As shown by Stengel et al. (1982) the viscosity-contrast-dependent critical Rayleigh number used to determine $\tilde \delta _c$ (Equation (25)) can be partially explained by assuming that the lower thermal boundary layer's thickness is described by

$$\begin{equation} \delta _c \approx \frac{c}{{\ln ({\Delta \eta })}} \cdot \tilde \delta _c.\end{equation} \tag{ 28 }$$

The constant c depends on the viscosity law (compare with Solomatov 1995). For a viscosity exponentially decreasing with temperature Stengel et al. (1982) find that c = 8, and hence Equation (28) suggests that the lower thermal boundary layer comprises only a fraction of the lower mantle conductive zone if the viscosity contrast through $\tilde \delta _c$ is above ∼3 × 10³. Adjusting the energy equations (Equation (12)) to include δ_c as part of the convective mantle using Equation (28), we find qualitatively the same trends as observed without a direct determination of the lower thermal boundary layer. We have modeled several cases accounting for this correction and find that in comparison to the standard 1D model the upper mantle is maximally ∼150 K (∼5%) hotter due to the larger amount of heat sources in the convective zone. This leads to stronger core heating or smaller core cooling rates, respectively. For the reference model the core is up to ∼50 K hotter after 13 Gyr of evolution for the pressure-dependent models B and C. As an additional consequence of the hotter convective mantle, the total lower mantle conductive zones are about 5%–15% smaller.

However, we find that the CMB-lids predicted by Equation (28) for the 1D model are a factor of two smaller than those suggested by 2D models for the same cases. Moreover, although it has been shown that the viscosity contrast through the upper thermal boundary layer is smaller than one order of magnitude (e.g., Davaille & Jaupart 1993) we do not obtain the same viscosity contrast through δ_c estimated with Equation (28). Instead we find that the viscosity contrast through the suggested lower thermal boundary layer starts initially at ∼10³–10⁴. For our reference model and massive planets with large CMB-lids such as the M = 10 super-Earth with V*_eff = 2.5 cm³ mol⁻¹ this viscosity contrast falls after 5 Gyr below one order of magnitude. For less massive planets or weaker pressure effects the viscosity contrast is only after ∼10 Gyr below one order of magnitude. This might suggest that the viscosity contrast through the lower thermal boundary layer reaches one order of magnitude in the steady-state. Our results, however, indicate that this limit cannot be reached in geologically relevant times.

Interestingly, the comparison between 2D and 1D CMB-lid thicknesses is far better when the constant of c ∼ 8 in Equation (28) is reduced to c ∼ 1. In this case the viscosity contrast through δ_c is closer to one order of magnitude, indicating that the lower thermal boundary layer, which is effectively taking part in the convection, is for our viscosity law indeed smaller than suggested by Equation (28). Future studies are required to better constrain the thickness of the lower thermal boundary layer associated with a temperature- and pressure-dependent Arrhenius-type viscosity.

The 1D model assumes that the characteristic Rayleigh number defining the lower mantle conductive zone equals the viscosity-contrast-dependent critical Rayleigh number ${\rm Ra}({\tilde \delta _c } ) = {\rm Ra}_{{\rm crit}} ({\Delta \eta } )$. To test the sensitivity of our results on the critical Rayleigh number, we additionally perform all runs with a constant value for ${\rm Ra}({\tilde \delta _c } )$, setting it to our isoviscous value of ${\rm Ra}({\tilde \delta _c } ) = {\rm Ra}_{{\rm crit}}^{{\rm iso}} = 500$—corresponding to a lower limit of the possible values of ${\rm Ra}({\tilde \delta _c } )$.

The main differences between a constant characteristic Rayleigh number model and our standard viscosity-contrast-dependent model are illustrated in Figure 17: for the standard reference temperatures we find almost no differences in the thermal evolution (the CMB temperature T_c(t) and the upper mantle temperature T_m(t) vary by less than 1% between the two models). The differences increase for larger initial CMB temperatures and smaller planets. Here we find smaller lower mantle conductive zones (by 5%–30%) and an increase of cooling rates (see Figure 17(a) for M = 5 in an SL mode, and Figure 17(b) for M = 1 in a PT mode). For super-Earths the core starts cooling at 100–1000 K lower initial CMB temperatures than for ${\rm Ra}({\tilde \delta _c } ) = {\rm Ra}_{{\rm crit}} ({\Delta \eta } )$, but nonetheless cooling is far less effective than for V*_eff = 0.

Zoom In Zoom Out Reset image size

In conclusion, the difference between 1D thermal evolution models with ${\rm Ra}({\tilde \delta _c } ) = {\rm Ra}_{{\rm crit}} ({\Delta \eta } )$ and ${\rm Ra}({\tilde \delta _c } ) = 500$ is not significant and all observed thermal trends and effects remain unaltered. These results indicate that the viscosity controls substantially the dynamics and the thermal evolution.

In the models presented we investigate the effects of a temperature- and pressure-dependent viscosity on the thermal evolution of super-Earths, but we neglect any temperature and pressure dependence of the thermal mantle conductivity k_m, the thermal mantle expansivity α_m and the mantle density ρ_m.. Their change with depth, due to pressure and temperature, is however smaller compared to the variation of the mantle viscosity with depth.

To test the impact of an additional variation of k_m, α_m, and ρ_m with temperature and pressure, we proceed to compute the upper thermal boundary layer, the upper mantle SL, and the heat flux out of the mantle as before with standard values for k_m, α_m (Table 1), and the average mantle density ρ_m. However, for the lower mantle, i.e., for the computation of the lower mantle conductive zone (Equation (26)), the lower mantle conductive zone heat equation (Equation (23)), and the heat flux into the mantle and out of the core (Equations (21) and (22)) we distinguish two simplified test cases:

• 1.  
  We use higher k_m values for the lower mantle (varied from the standard value of k_m = 4 W m⁻¹ K⁻¹ to 20 W m⁻¹ K⁻¹ and 40 W m⁻¹ K⁻¹) but keep our standard α_m and the average mantle density ρ_m.
• 2.  
  We use the results from Stamenković et al. (2011), who computed k_m, α_m, and ρ_m for an Earth-like 80 vol.% MgSiO₃ perovskite and 20 vol.% MgO periclase mantle composite (in fact we assume periclase to only significantly alter the phonon thermal conductivity), and use for every time step of the thermal evolution average values of k_m, α_m, and ρ_m evaluated at the average temperature and pressure in $\tilde \delta _c$.

We generally observe that with increasing thermal conductivities the lower mantle cools more efficiently, leading to larger $\tilde \delta _c$ due to the larger viscosity. This results in an increasingly linear temperature profile in this zone (Figure 18 for M = 10 in an SL mode). However, the more efficient cooling of the lower mantle does not imply a stronger cooling of the core. In fact, with increasing thermal conductivity heat is transported also more efficiently into the core, which results in an even stronger heating of the core (Figure 18). There is, however, a threshold, i.e., increasing the thermal conductivity of the lower mantle above a critical value would cause the core to only cool. For the chosen example in Figure 18 this critical conductivity is ∼100 W m⁻¹ K⁻¹, which is about four times larger than the thermal conductivity suggested by Stamenković et al. (2011) for the corresponding pressure and temperature regime.

Zoom In Zoom Out Reset image size

The results of Stamenković et al. (2011) suggest for the specific M = 10 planet (Figure 18) that the average total thermal conductivity in $\tilde \delta _c$ (phonon + radiation + electrons) is about 20–25 W m⁻¹ K⁻¹ (5–6 times larger than our standard value of 4 W m⁻¹ K⁻¹), that the average thermal expansivity is ∼4 × 10⁻⁶ K⁻¹ (∼5 times smaller than our standard value) and that the average density in $\tilde \delta _c$ is similar to the average mantle density ρ_m. We find the effects of the five-times-smaller thermal expansivity to be minor compared to the effects of the five-times-larger thermal conductivity in comparison to our standard model. The model with varied k_m, α_m, and ρ_m is similar to the k_m = 20 W m⁻¹ K⁻¹ model in Figure 18. This is due to the fact that the thermal expansivity only affects the thickness of the lower mantle conductive zone $\tilde \delta _c$ (Equations (24) and (26)), but the thermal conductivity governs as well the heat transport through $\tilde \delta _c$ (Equation (23)), the heat flux out of the core (Equation (22)), and the heat flux into the convecting mantle (Equation (21)).

In general, we conclude that the temperature and pressure dependence of k_m, α_m, and ρ_m only strengthens our conclusion that a stagnant layer at the CMB can be present in super-Earths. In particular, higher values of the thermal conductivity in the lower mantle result in more efficient cooling of the lower mantle and hence lead to increasing $\tilde \delta _c$. This even further supports the idea that heat in the deep interior of super-Earths is transported via conduction rather than by convection for lower initial interior temperatures.

The propensity of super-Earths to have PT has been controversially discussed in the recent past (Korenaga 2010; O'Neill & Lenardic 2007; Stein et al. 2011; Valencia et al. 2007; Valencia & O'Connell 2009; Van Heck & Tackley 2011). In these models, the pressure dependence of the viscosity is generally neglected and it is assumed that the entire mantle convects vigorously. Our present models show that if the pressure dependence of the viscosity is considered the effectively convecting mantle can be significantly reduced in thickness depending on initial conditions. Hence, one can also expect the propensity of PT to be affected, which is influenced, among other factors, by the yield stress τ_y and the convective shear stress τ_c beneath the lithospheric plate of thickness L. A necessity for PT is plate yielding, which can only be achieved when the convective stress overcomes the yield stress of the lithospheric plate. Due to the fact that for V*_eff > 0 we generally observe a reduction of convective vigor with increasing planetary mass, as well as a reduction of convective mantle thickness, the mantles of super-Earths with large enough CMB-lids resemble in their convective structure actually less massive planets. Assuming that for V*_eff = 0 less massive planets have smaller chances of driving PT, as has been suggested by Valencia et al. (2007), would lead to the conclusion that for massive super-Earths the propensity of PT could start to decrease with increasing mass. Interestingly, O'Neill & Lenardic (2007) suggest a decrease of surface mobility for super-Earths already for viscosities with nearly V*_eff = 0 and Stein et al. (2011) even for a decrease in lower mantle viscosity due to a change of diffusion mechanism from vacancy to interstitial control as suggested by Karato (2011). However, as already discussed, for V*_eff > 0 the convective structures in super-Earths strongly depend on the planet's initial conditions and hence it might be difficult to obtain a relation between the propensity of PT and the planetary mass.

Moreover, independent of viscosity law and the mass of the planet, the propensity of PT might be even stronger influenced by the amount of volatiles, i.e., water in the lithosphere, the age of the planet, lithospheric buoyancy as suggested by Kite et al. (2009), and most of all by the amount of radiogenic heat sources, which might crucially vary from planet to planet and could as well depend on the time of planetary formation in relation to the age of the universe.

Classic dynamo generation requires convection in a partially liquid metallic core by thermal or compositional buoyancy, which on the other hand requires that the core is cooling (for a description of thermal and compositionally driven dynamos, see, e.g., Stevenson 2009). The present results indicate that core cooling in super-Earths is strongly reduced for pressure-dependent rheologies. In fact for our reference models, assuming temperatures that are currently suggested for super-Earths (Papuc & Davies 2008; Sotin et al. 2007; Valencia et al. 2006), the core is continuously heating—preventing dynamo generation for these cases. For higher initial temperatures, the core may cool with time but possibly with a heat flux below the critical heat flux necessary for driving a thermal dynamo. Such behavior is already discussed by Tachinami et al. (2011) for super-Earth conditions with an olivine rheology. However, this does not generally exclude chemical dynamo generation, for which the growth of a solid iron core could induce convection and hence a magnetic field—in that case the core also needs to cool, but the rate can be smaller than for a thermal dynamo and the core temperature needs to be partially (but not completely, see below) under the liquidus of the iron-rich core. Moreover, the steep melting curves for MgSiO₃ perovskite and MgO by Stamenković et al. (2011) associated with large effective activation volumes for mantle rock also indicate large melting gradients for iron. Hence, additionally, the melting temperatures of iron could strongly increase with pressure and make it more difficult for more massive planets to have molten iron cores, a necessity for dynamos. Thus, although we cannot exclude magnetic field generation in super-Earths, assuming a temperature- and pressure-dependent viscosity strongly limits the conditions preferable for dynamo action to a small parameter range. Note that this conclusion does not change if the planet is in the PT regime (compare with Section 5.2).

Present observations of super-Earths are mainly limited to mass and radius and do not provide direct information on interior dynamics. The atmospheric chemistry, however, which could be investigated by combined future observations (as with Spitzer, Hubble, the possibly upcoming EChO mission, or with the James Webb Space Telescope), is influenced by the planet's outgassing history and volcanism—with both being affected by the interior composition, structure, and rheology.

If the viscosity is strongly pressure-dependent then super-Earths might be even undifferentiated, i.e., no separation of an iron core and a silicate mantle, as indicated by Stamenković et al. (2011)—note that for this study we assume that differentiation took place. An undifferentiated interior would lead to a potentially different outgassing history, and hence atmospheric composition, in comparison to a differentiated body, which is more likely to form when the effective activation volume is small.

Assuming, however, a differentiated interior, we find that for SL planets melt production decreases with both increasing planetary mass and effective activation volume, indicating a reduced volcanic activity for SL super-Earths in comparison to an SL M = 1 planet. Using the peridotite solidus data from Hirschmann (2000), the duration of upper mantle melt production for SL planets decreases with planetary mass for V*_eff = 0: melt production is possible for an M = 1 planet until 9 Gyr after the planet's formation, but an M = 10 super-Earth can produce melt only during the first 6 Gyr. This reduced duration of melt production for V*_eff = 0 with increasing planetary mass can be mainly explained by the increase of gravity that enlarges the radial gradient of the melting curve, and additionally by the decreasing peak upper mantle temperatures for larger M (see Figure 7).

An increase of the effective activation volume does not severely alter the upper mantle temperature for a specific planet (see Figure 7) but leads to a strong increase of the SL thickness for super-Earths (when V*_eff is increased from V*_eff = 0 to V*_eff = 2.5 cm³ mol⁻¹ the lid thickness L for M = 1 is almost unaltered, but for M = 5, 10 L increases by a factor of two or four, respectively. This is due to the higher viscosity, which is a consequence of the increased pressure effect). This leads to an additional reduction of the thermal gradient in the upper mantle for M = 5, 10 planets in comparison to V*_eff = 0, hence melt production in the upper mantle is even more suppressed for increasing V*_eff. For M = 5 and M = 10 melt in the upper mantle can only be produced during the first 4 Gyr or the first 1 Gyr, respectively. A basal magma ocean seems unlikely due to the large MgSiO₃ perovskite and MgO melting temperatures found by Stamenković et al. (2011).

Hence if super-Earths are more likely to be in an SL regime, as suggested by our findings, then volcanism ceases earlier on those planets. Note that in case of PT planets volcanism is an ongoing process due to the continuous formation of oceanic and continental crust. Therefore, if detailed spectroscopic data from super-Earths become feasible, then the atmospheric composition might contribute to distinguish different rheological models.