(2021) 73:167
Kameyama ﻿Earth, Planets and Space
https://doi.org/10.1186/s40623-021-01499-w

Open Access

FULL PAPER

Linear analysis on the onset of thermal
convection of highly compressible fluidswith variable viscosity and thermal conductivity
in spherical geometry: implications
for the mantle convection of super‑Earths
Masanori Kameyama1,2*

Abstract
In this paper, we carried out a series of linear analyses on the onset of thermal convection of highly compressible
fluids whose physical properties strongly vary in space in convecting vessels either of a three-dimensional spherical shell or a two-dimensional spherical annulus geometry. The variations in thermodynamic properties (thermal
expansivity and reference density) with depth are taken to be relevant for the super-Earths with ten times the Earth’s
mass, while the thermal conductivity and viscosity are assumed to exponentially depend on depth and temperature,
respectively. Our analysis showed that, for the cases with strong temperature dependence in viscosity and strong
depth dependence in thermal conductivity, the critical Rayleigh number is on the order of 1­ 08–109, implying that the
mantle convection of massive super-Earths is most likely to fall in the stagnant-lid regime very close to the criticalcondition, if the properties of their mantle materials are quite similar to the Earth’s. Our analysis also demonstrated
that the structures of incipient flows of stagnant-lid convection in the presence of strong adiabatic compression are
significantly affected by the depth dependence in thermal conductivity and the geometries of convecting vessels,
through the changes in the static stability of thermal stratification of the reference state. When the increase in thermal
conductivity with depth is sufficiently large, the thermal stratification can be greatly stabilized at depth, further inducing regions of insignificant fluid motions above the bottom hot boundaries in addition to the stagnant lids along
the top cold surfaces. We can therefore speculate that the stagnant-lid convection in the mantles of massive superEarths is accompanied by another motionless regions at the base of the mantles if the thermal conductivity strongly
increases with depth (or pressure), even though their occurrence is hindered by the effects the spherical geometriesof convecting vessels.
Keywords: Mantle convection, Super-Earths, Adiabatic compression, Temperature-dependent viscosity, Depthdependent thermal conductivity, Critical Rayleigh number, Stable thermal stratification*Correspondence: kameyama@sci.ehime-u.ac.jp
2
Japan Agency for Marine-Earth Science and Technology (JAMSTEC),
Yokohama, Kanagawa 236‑0001, Japan
Full list of author information is available at the end of the article

Introduction
Super-Earths are extra-solar planets which have small
masses (compared to that of Jupiter) and high mean density larger than 5000 kg/m3 (Howard et al. 2010). Since
the first discovery of such massive terrestrial planets(Rivera et al. 2005), a growing effort has been devoted tothe studies in the internal (radial) structure (e.g., Seager© The Author(s) 2021. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which
permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to theoriginal author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images orother third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line
to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory
regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of thislicence, visit http://​creat​iveco​mmons.​org/​licen​ses/​by/4.​0/.

Kameyama ﻿Earth, Planets and Space(2021) 73:167

et al. 2007; Valencia et al. 2007b; Sotin et al. 2007; Wagner et al. 2011), the thermal evolution (e.g., Papuc andDavies 2008; Kite et al. 2009; Korenaga 2010; Gaidos et al.
2010; Tachinami et al. 2011; Čížková et al. 2017) and thehabitability (e.g., Noack et al. 2017; Tosi et al. 2017; Foley2019) of super-Earths. During the last decades, the studies on mantle dynamics have been playing an importantrole in addressing these issues (e.g., Shahnas et al. 2018;van den Berg et al. 2018).
It can be intuitively understood that, because of theirlarge sizes, the convection in the mantles of massive
super-Earths is most likely to occur in a quite different
manner from those in the Earth and smaller terrestrial
planets. For example, the variations of physical propertiesof mantle materials are expected to be larger in the muchthicker mantles of super-Earths than in the Earth’s (e.g.,Karato 2011). Indeed, the effects of strong spatial variations in physical properties (such as viscosity and thermal conductivity) on the convecting flow patterns havebeen widely studied so far, under the assumption that theeffects of other physical mechanisms such as compressibility can be minor (e.g., Kameyama and Ogawa 2000;
Miyauchi and Kameyama 2013). In particular, the effectsof strong temperature dependence in viscosity have been
intensively studied by numerical models of convectionin three-dimensional spherical geometry (e.g., Yao et al.
2014; Yanagisawa et al. 2016; Guerrero et al. 2018).
On the other hand, it has been well acknowledged thatthe effect of adiabatic compression affects the nature ofthermal convection (Jarvis and McKenzie 1980; Kinget al. 2010; Liu and Zhong 2013). This effect is expected
to be of great importance in the convection in the thickmantles of massive super-Earths, considering that thepressure can reach on the order of terapascals in theirinteriors (Tsuchiya and Tsuchiya 2011). Indeed, our earlier studies (Kameyama and Kinoshita 2013; Kameyama
et al. 2015; Kameyama 2016) showed that the nature ofthermal convection in the thick mantles of massive superEarths is strongly affected by the interplay between theeffects of adiabatic compression and those of spatial variations in physical properties of mantle materials. However, in these studies, the effects of spherical geometriesof planetary mantles have been thoroughly ignored, since
the analyses had been carried out using the convecting
vessels with the Cartesian (planar) domains. It is therefore very important to estimate the effects of the strong
adiabatic compression and those of strong variations inphysical properties of mantle materials on the convecting
flow structures in spherical geometries, in order to gain
the fundamental insights into the mantle convection ofmassive super-Earths.
In this study, we carry out a series of linear analyses on
the thermal convection in spherical domains of a highly

Page 2 of 17

compressible fluid with the spatial variation in physical properties, to investigate how the spherical geometries of planetary mantles can affect the critical stateand condition of thermal convection in the presence of
the interplay between the adiabatic compression andspatial variations in physical properties such as viscosity and thermal conductivity. In particular, we will focus
on the influences of adiabatic changes in temperatureon the incipient convection, by using a simple 1D modelof a basally heated fluid. By further comparing the present results with our earlier one in the Cartesian (planar) model (Kameyama 2016), we aim at deepening theinsights into the thermal convection in the mantles ofmassive terrestrial super-Earths both from the thermodynamic and fluid-dynamic aspects.

Model descriptions
The onset of thermal convection is considered for ahighly compressible fluid of infinite Prandtl number, asa model of mantle convection of a super-Earth with ten
times the Earth’s mass. The numerical model is the sameas those of our earlier one (Kameyama 2016) except thatthe convecting vessels are taken to be either a threedimensional (3D) spherical shell or a two-dimensional(2D) spherical annulus whose inner and outer radii arermin and rmax , respectively, and rmin /rmax = 0.5. The temperature T of the fluid layer is fixed to be T = Ttop and
Tbot (> Ttop) at the top and bottom boundaries, respectively. In order to keep the present analysis as simple
as possible, we ignored the effects of internal heating
and spatial variations in specific heat Cp of the fluid andthe gravity g throughout the layer. In addition, we fixed
Ttop /�T = 0.1 (where T ≡ Tbot − Ttop is the temperature difference across the layer), indicating that the temperature at the top surface of the planet is similar to the
Earth’s.
Physical properties of modeled compressible fluid

As a model of pressure dependence in thermodynamic
properties of mantle materials of massive super-Earths,
we employ the same variations in thermal expansivity
α and reference density ρ with depths as in our earlier
works (Tachinami et al. 2014; Miyagoshi et al. 2014, 2015;
Kameyama 2016; Kameyama and Yamamoto 2018). The
thermal expansivity α decreases with depth, and its radial
profile is taken to be similar to that estimated for MgOfrom an ab initio calculation (Tsuchiya and Tsuchiya2011) for the mantles of super-Earths with ten times
the Earth’s mass. That is, the value of α(r = rmin ) ≡ αbot
at the bottom boundary is as small as about 1/14 of thevalue α(r = rmax ) ≡ αtop (= 4 × 10−5 K−1) at the top surface, and the value of its average αave is taken to beKameyama ﻿Earth, Planets and Space

αave
1
=
αtop
rmax − rmin

(2021) 73:167

 r=rmax
r=rminα(r)dr ≃ 0.1574.
αtop

Page 3 of 17(1)

On the other hand, the reference density ρ is assumed
to linearly increase with depth, yielding the valueof ρ(r = rmin ) at the bottom boundary larger thanρ(r = rmax ) ≡ ρtop at the top surface by about 3.17
(Valencia et al. 2006). We note that the linear increasein ρ assumed here is not likely to spoil the significanceof the present analysis despite its geophysical unsoundness, since the vertical profile of (reference) density has avery minor effect on the nature of critical states of thermal convection of highly compressible fluids (Kameyama
et al. 2015).
In addition to the variations in thermodynamic properties (α and ρ ), we also assume the spatial variations intransport properties. In this study, the thermal conductivity k is taken to be exponentially dependent on depthas,rmax − r,k(z) = ktop exp ln(rk ) ×(2)
rmax − rmin
where ktop is the thermal conductivity at r = rmax ,and rk is a dimensionless constant equal to the ratio ofk(r = rmin ) ≡ kbot to ktop. For a geophysical significance,
we restrict ourselves to the cases where rk ≥ 1; the thermal conductivity k increases with depth for rk > 1.
According to the recent estimates both by theoretical calculations (Dekura et al. 2013; Dekura and Tsuchiya 2017,2019) and by laboratory experiments (Haigis et al. 2012;
Ohta et al. 2012, 2017; Lobanov et al. 2017), the thermal
conductivity of mantle materials can increase by up to
about one–two orders of magnitude across the mantles
of super-Earths. In this study, we varied rk up to 100.
On the other hand, the viscosity η is assumed to be
exponentially dependent on temperature T as,
Tbot − T,
η(T ) = ηbot exp ln(rη ) ×(3)
Tbot − Ttop
where ηbot is the viscosity for T = Tbot , and rη is a dimensionless constant given by rη = η(T = Ttop )/η(T = Tbot ).
In (3), for simplicity, we ignored the dependence on pressure (or depth) in viscosity, since it is still highly controversial if the viscosity of mantle materials increases or not
as the pressure increases up to the levels relevant to theinterior of super-Earths (Stamenković et al. 2011; Karato
2011; Wagner et al. 2012; Tackley et al. 2013; Noack andBreuer 2014). In this study, we restrict ourselves to thecases where rη ≥ 1; the viscosity η decreases with temperature for rη > 1. Using this assumption, we will focus on
the effects of the reduction of viscosity in the deep mantle of the super-Earths proposed in several earlier studies(Karato 2011; Stein et al. 2011; Wagner et al. 2012; Tackley et al. 2013).
Here, some comments are necessary on the depth
dependence in thermal conductivity k in this study. First,
the exponential increase in k with depth given by (2) ischosen simply from analytical and numerical convenience
as in our earlier studies (Kameyama and Kinoshita 2013;
Kameyama 2016; Kameyama and Yamamoto 2018) and,in other words, may not be necessarily supported by studies on mineral physics conducted so far (e.g., Dalton et al.
2013; Hsieh et al. 2017, 2018; Deschamps and Hsieh 2019).
Second, it is well acknowledged that the thermal conductivity k depends not only on depth (or pressure) but also onother parameters such as composition and temperature.
For example, the thermal conductivity k of ferro-periclase
is much higher and more sensitive to pressure (about a factor of 10 increase in the pressure range of the Earth’s mantle) than that of bridgmanite. Similarly, the values of k ofboth minerals are sensitive to the iron content. In addition,
the temperature dependence in k is expected to compensate for the increase in k with depth, since high temperatures tend to reduce the contributions of heat transport bylattices (or phonons). In short, the present study may overestimate the values of rk ≡ kbot /ktop in silicate mantles ofsuper-Earths by ignoring the dependence of k on composition and temperature.
Basic equations

Under the truncated anelastic liquid approximation(TALA), the basic equations can be written in a dimensionless form as:

0 = ∇ · (ρv),

(4)

 

2
0 = −∇p + η ∇ ⊗ v + v ⊗ ∇ − (∇ · v)I· ∇ + Raρ αT er ,
3

(5)


Di
∂T
ρ
= ∇ · k∇T − ρvT − Diρ α T vr +η(ε̇II )2 ,
∂t
Ra
(6)

where the quantities with overbars denote those of reference state which vary only in the radial direction,and the nondimensionalization is carried out with thelength scale of d ≡ rmax − rmin , time scale of d 2 /κtop(where κtop = ktop /ρtop Cp is the value of thermal diffusivity at the top surface), and temperature scale ofT ≡ Tbot − Ttop. Here, ⊗ stands for the tensor productof vectors, v ⊗ ∇ is the transpose of ∇ ⊗ v , I is the identity tensor, er is the unit vector in radial direction positive outward, v is the velocity vector, vr ≡ (v · er ) is thevelocity in the radial direction, p is pressure, and ε̇II is thesecond invariant of strain rate tensor. In the right-handside of (6), the first, second, and third terms represent theKameyama ﻿Earth, Planets and Space(2021) 73:167

Page 4 of 17
effects of apparent heat transport (advection plus conk = exp [ln(rk ) × (rmax − r)], η = exp ln(rη ) × Tbot − T ,duction), adiabatic temperature change, and viscous dis(11)sipation (frictional heating), respectively.
while those for infinitesimal perturbations are given by,
The present formulation includes two nondimensional parameters. The first one is the dissipation number Di of
0 = ∇ · ρv ′ ,(12)
the modeled fluid defined by,
 

2
(13)
0 = −∇p′ + η ∇ ⊗ v ′ + v ′ ⊗ ∇ − (∇ · v ′ )I· ∇ + Raρ αT ′ er ,3

αtop gd
.
Di ≡
Cp

(7)

The second one is the Rayleigh number Ra of thermal
convection defined by,Ra ≡ρtop αtop �Tgd 3.
ηbot κtop(8)

Note here that, as indicated by (3), the viscosity is scaledby the value ηbot at the bottom boundary while for otherproperties the reference values are those at the top surface. The values of Di and Ra for the mantles of terrestrial
planets are expected to increase with the mass M of theplanets, coming from the increase in the mass, gravity,
and the mantle thickness of the terrestrial planets whose
composition is similar to the Earth’s (e.g., Seager et al.
2007; Valencia et al. 2007a; Wagner et al. 2011). According to the scaling relationship [ g ∝ M 0.5 and d ∝ M 0.28;
Valencia et al. (2007a)] together with Di ∼ 0.7 and
Ra = O(107 ) for the Earth’s mantle (Schubert et al. 2001),
one can estimate Di ∼ 5 and Ra = O ­(108–109) for the
mantles of super-Earths whose mass is ten times larger
than the Earth’s (Tachinami et al. 2014; Miyagoshi et al.
2014, 2015, 2018). In the present study, we fixed Di = 5,
while the values of Ra estimated above can be comparedwith those of the critical Rayleigh numbers derived fromthe linear stability analysis.
In the linear stability analysis presented in this paper, we
solve for the temporal evolution of an infinitesimal perturbation superimposed to a reference state described by astationary (motionless) state with steady one-dimensional
heat conduction in the radial direction. In the following,
the quantities with overbars and primes denote those ofreference state and perturbation, respectively. The dimensionless equations for the reference state are written as,(9)v = 0,∇· k∇T = 0,

T (r = rmax ) = Ttop ,T (r = rmin ) = Tbot = 1 + Ttop ,

(10)

 dT

∂T ′
1
= ∇ · k∇T ′ −v′ − Diα T vr′ .
∂t
ρ
dr r

(14)

In deriving these equations, the second-order terms ofinfinitesimal (primed) quantities are ignored. This is thereason why the effect of viscous dissipation is eliminatedin (14).
At the top and bottom boundaries, we kept the perturbations in temperature T ′ and the radial velocity vr′ tobe zero. For the horizontal velocity at the boundaries,
on the other hand, we assumed a free-slip condition.
Numerical method

The equations for infinitesimal perturbations (12) to(14) are recast into the one-dimensional differential
equations in the radial direction, by applying the spectral expansion in the horizontal planes using either the
spherical harmonics or sinusoidal function dependingon three- and two-dimensional spherical geometries.
The derived one-dimensional differential equations
are discretized with an eighth-order accuracy in space
using equally spaced 257 grid points. The eigenvalue
problem for a given condition (such as the sphericalharmonic degree ℓ or the horizontal wavenumber Kof perturbations) is solved for a critical Rayleigh number Rac as well as the radial profiles of eigenfunctions
for infinitesimal perturbations by a numerical method
developed in Kameyama et al. (2015). By increasing
the values of ℓ or K, we search the one which yields
the absolute minimum of the critical Rayleigh number Rac0 and the values of ℓc0 or Kc0 corresponding toRac0 . See Appendix for the detailed descriptions for thederivation and solution of one-dimensional differential
equations.

Results and discussion
Static stability of thermal stratification for reference(steady) states

We first see how the structure and static stability ofthe thermal stratification of the reference state varydepending on the model geometries and on the depthKameyama ﻿Earth, Planets and Space(2021) 73:167

Page 5 of 17

dependence in thermal conductivity k . By the solid lines
in Fig. 1a, we show the radial profiles of reference temperature T (r) derived from Eq. (10) for different model
geometries and various values of rk . [Note that the temperature dependence in viscosity η affects neither thethermal state nor its static stability of the reference state
given by Eq. (10)]. Figure 1b shows the magnitude of its
radial gradient |dT /dr|. To further highlight the effects ofmodel geometries, we also show in Fig. 1 the similar plotsof T and its vertical gradient dT /dz for the Cartesiangeometry derived in our earlier work (Kameyama 2016).
From the plots in Fig. 1a, b, we can see that the radial
profiles of reference temperature T are significantly
affected by the values of rk as well as by the model geometries. The comparison of the plots for different rk andgiven geometries shows that the values of T in the fluidlayer are higher for larger rk . The increase in T for large rk
comes from the changes in the radial profiles of |dT /dr|within the fluid layer depending on rk [see also Eq. (10)]:a larger increase in k with depth induces a larger decreasein |dT /dr| with depth. This in turn makes |dT /dr| larger
in the shallow cold part and smaller in the deep hotpart of the fluid layer (see Fig. 1b), which consequently
raises T in the entire layer (see Fig. 1a). The comparisonof the plots of T for 3D and 2D spherical geometries, on
the other hand, shows that for a given rk the values of T
in the fluid layer are lower in 3D geometry than in 2D
geometry. The difference in T due to the model geometries comes from the difference in the changes in thesurface area with depth. The decrease in the surface areawith depth is more significant in 3D cases (∝ r −2) than in
2D cases (∝ r −1). Owing to this effect, the magnitude of|dT /dr| in the shallow part of fluid layer is smaller in 3D
cases than in 2D cases and, at the same time, the magnitude in the deep part is larger (see Fig. 1b). This results inlower T in the entire layer for 3D cases than in 2D cases.
The effect of the changes in the surface area with depth
can be further confirmed from the comparison with the
Cartesian cases (Kameyama 2016). Among the threemodel geometries presented in Fig. 1, indeed, the valuesof T for given rk are the highest for Cartesian cases where
the effect of the changes in the surface area with depth isabsent.
To see the static stability of the thermal stratification
more clearly, we show in Fig. 1c the radial profiles ofpotential temperature T pot of the reference states. Here,T pot at the dimensionless position r is calculated by ′r =rmax

T pot (r) = T (r) exp −Di

r ′ =rα(r ′ )dr ′ .

(15)

We also show by thick gray lines the ranges of depths
where T pot increases with height. Note that in theseregions the thermal stratification given by T is stable, as
is inferred from the Brunt–Väisälä frequency N: Here thedimensionless value of N is given in the present numerical configuration by






 dT 
 −  dT  ,N 2 (r) = α(r) 
(16)


dr s
dr 
where a nondimensionalization
is carried out by a fre
quency scale of g/d , and (dT /dr)s is the adiabatic temperature gradient whose dimensionless value is defined
by


dT
≡ −Di α T .
(17)
dr s

In the regions where N 2 > 0 (i.e., N is real), the thermal
stratification is stable against an imposed infinitesimal
displacement of fluid parcels in the direction of gravity.
Indeed, by differentiating (15) with respect to r we get






 dT 
dT pot −  dT  T pot ,
= 
(18)


dr
dr s
dr  T

which indicates that dT pot /dr is positive when|dT /dr|s > |dT /dr| and, hence, N 2 > 0. In addition, thestatic stability of thermal stratification is stronger forlarger N 2 and dT pot /dr.
The comparison between the plots in Fig. 1c clearlyshows that the static stability of the thermal stratification is affected not only by the values of rk but also by thegeometries of fluid layer. For both 3D and 2D sphericalgeometries, regions of stable thermal stratification withdT pot /dr > 0 emerge in the deeper part of the fluid layer
when rk ≥ 30. This is because a large rk reduces the magnitude of |dT /dr| in the deeper part and, hence, enhances
the static stability of thermal stratification [see also Eq.
(18)]. These regions become thicker and the thermal
stratification there becomes more stable for larger rk . We
note, on the other hand, from the comparison between
3D and 2D cases that for given rk ≥ 30 the regions ofstable stratification are thinner for 3D cases than thosefor 2D cases. This is again due to the difference in thechanges in the surface area with depth. For 3D cases, themagnitude of |dT /dr| in the deeper part of the fluid layer
becomes larger than in 2D cases, which in turn reduces
the static stability of thermal stratification at depths. The
effect of model geometries on the static stability can befurther confirmed from the comparison with the Cartesian cases (Kameyama 2016) where the emergence ofregions of stable stratification is greatly enhanced withrespect to that with spherical cases.
We will conclude this subsection by estimating thecondition under which the regions of stable stratification

Kameyama ﻿Earth, Planets and Space(2021) 73:167

Page 6 of 17

Fig. 1 Thermal state in the fluid layer of different model geometries estimated for various values of rk in the range of 1 ≤ rk ≤ 100. a Plots against
height from the bottom surface of the reference temperature T . Also shown by thin dashed lines for comparison are the plots of several adiabatsagainst height. b Plots against height of the magnitude of the gradient of temperature |dT /dr| (or |dT /dz|) in the radial (vertical) directions. c Plotsagainst height of the potential temperature T pot. Also shown by thick gray lines are the ranges of depths where the thermal stratification given by T
is stable owing to the increase in T pot with heightKameyama ﻿Earth, Planets and Space(2021) 73:167

Page 7 of 17

emerge at the bottom of thefluid layer. This condition is
given by dT /dr s > dT /dr  at r = rmin . Using (17), thisinequality can be further rewritten as dT 1

Di > Dithres ≡.
(19)
αbot Tbot  dr r=rmin

In Table 1, we summarize the values of dT /dr r=r calmin
culated for several rk and 3D spherical shell geometry
with rmin /rmax = 0.5 (see Fig. 1a, b), together with theresulting values of Dithres. In calculating Dithres we used
αbot ∼ 1/14 and Tbot = 1.1 in dimensionless units. In
Table 1 we also show the values of the threshold mass
Mthres of planets above which Di > Dithres in their mantles, based on the scaling relationship of Di ∝ gd ∝ M 0.78
described in the previous section. Note, however, that inthese estimates we employed the same decrease in thermal expansivity α with depth as that for M = 10ME
(where ME is the Earth’s mass). That is, these estimates
may not be very accurate for the cases with M = 10ME.
In particular, the values of Dithres and Mthres are most
likely to be smaller than those estimated here for thecases with M < 10ME, where both αave /αtop andαbot /αtop are expected to be larger than those for
M = 10ME and, in other words, the effects of adiabatic
compression are more significant.
Linear stability analyses for the onset of thermal
convection

Next we carry out linear stability analyses for the onsetof thermal convection of highly compressible fluids in thepresence of spatial variations in physical properties (temperature dependence of viscosity η and depth dependence in thermal conductivity k) in the convecting vessels
with different geometries. In particular, we compared thechanges in the critical states and structures of incipient
Table 1 The values of the temperature gradient |dT /dr|r=rminat the bottom of the fluid layer of 3D spherical shell ofrmin /rmax = 0.5 calculated with several values of rk = kbot /ktop,
the threshold values of dissipation number Dithres and planetary
mass Mthres divided by the Earth’s mass ME above which regionsof stable thermal stratification occur at the bottom of the fluid
layer
rk





dT /dr Dithres

Mthres /ME

r=rmin

1

2.000

26.14

83.34

3

1.246

16.28

45.43

10

0.6593

8.616

20.09

24.24

0.38265

1030

0.3327

4.348

8.360

100

0.1429

1.867

2.829

flows with those in the static stability of thermal stratification for the reference state, to investigate how thenature of incipient convection of compressible fluids withspatial variations in physical properties is affected by the
adiabatic temperature change and, in other words, theoccurrence of regions of stable stratification.
In Fig. 2, we show the plots against the viscosity contrast rη across the fluid layer of (a) the absolute minimum
value of critical Rayleigh number Rac0 and (b) the spherical harmonic degree ℓc0 or the horizontal wavenumber
Kc0 corresponding to Rac0, obtained for different model
geometries and various values of rk . Figure 2a shows thatin general the values of Rac0 become larger for larger rη
and rk , which can be intuitively understood from an overall increase in viscosity and thermal conductivity with rη
and/or rk . Figure 2b shows, on the other hand, that thevalues of ℓc0 and Kc0 rapidly increase with rη for rη  O­(103–104), regardless of the model geometries and thevalues of rk . This increase in ℓc0 or Kc0 with rη has beenrecognized by the transition of “stagnant-lid” regime ofthermal convection (Stengel et al. 1982; Kameyama et al.
2013). Taken together with the plots in Fig. 2a, the values of Rac0 for the stagnant-lid regime and rk ≥ 30 are aslarge as O ­(108–109), which is very close to the estimateof the Rayleigh number of mantle convection of superEarths whose mass is about ten times the Earth’s. Thissuggests that in massive super-Earths the stagnant-lid
type of mantle convection may occur under the conditions very close to the critical ones. In addition, among
the three geometries presented here, the 3D sphericalgeometry needs the largest values of rη for the transition into “stagnant-lid” regime, which is consistent withan earlier finding by linear (Kameyama et al. 2013) andnonlinear (Guerrero et al. 2018) studies on the effects ofcurvature of fluid layer.
We next focus on the nature of critical conditions forstrong temperature dependence in viscosity (rη  O­(103–104) where the thermal convection falls in the stagnant-lid regime. For both the 2D and 3D spherical geometries, the changes in Rac0 and those in ℓc0 and Kc0 withrη are different depending on whether rk ≤ 10 or rk ≥ 30
and, hence, whether the regions of stable thermal stratification occur or not in the reference state (see also Fig. 1).
For the cases with rk ≤ 10 lacking in stable regions, thevalues of Rac0 become almost independent of rk for largerη (Fig. 2a), and the values of ℓc0 and Kc0 are smaller forlarger rk (Fig. 2b). For the cases with rk ≥ 30 in the presence of stable regions, by contrast, the values of Rac0 arealways larger for larger rk , and the values of ℓc0 and Kc0
are smaller for larger rk . We also note from the plots for
the Cartesian geometry that the similar changes in Rac0and Kc0 with rη occur depending on whether rk ≤ 10 orKameyama ﻿Earth, Planets and Space(2021) 73:167

Page 8 of 17

Fig. 2 The plots against rη of a the absolute minimum value of critical Rayleigh number Rac0 and b the spherical harmonic degree ℓc0 or
horizontal wavenumber Kc0 of perturbation corresponding to Rac0, obtained for different model geometries and various values of rk in the range of1 ≤ rk ≤ 100

Kameyama ﻿Earth, Planets and Space(2021) 73:167

Page 9 of 17

rk ≥ 30 and, in other words, whether the regions of stable
thermal stratification exist or not.
To see how the depth dependence in thermal conductivity k affects the critical states of stagnant-lid convection, we show in Fig. 3 the schematic illustration on
two-dimensional planes of the incipient flows obtained
for several values of rk and various model geometries
together with rη = 108. (For the cases in 3D sphericalgeometry, we show the flow structures in meridian planes
assuming the incipient flows occurring with only zonal
components whose spherical harmonic order is m = 0.)
Shown in color are the distributions in the perturbationsin temperature T ′ whose color scales are indicated at thebottom of the figure. Shown by solid lines are the contours of a “potential” ψ of the mass flux ρv ′ which satisfies
ρv ′ = ∇ × (ψen ),(20)

where en is the unit vector perpendicular to the plane.
(Note that the function ψ is a generalization of thestream function for the Boussinesq or extended Boussinesq approximation where ρ = 1.) As can be seen from
Fig. 3a, the incipient flows for rk = 1 occur as a wellknown stagnant-lid convection. Indeed, thick lids of coldand highly viscous fluid develop along the top surface
where the fluid motion is insignificant. The convection
occurs only beneath the thick and stagnant lids, which
results in a small vertical extent of convection cells. The
horizontal extent of convection cells also becomes muchsmaller than the thickness of entire layer, in accordancewith the reduction in their vertical extent. This resultsin large ℓc0 and Kc0 of the infinitesimal perturbation (seealso Fig. 2b).
The comparison of Fig. 3b and c with 3a clearly showsthat the structures of incipient flows of stagnant-lid convection are significantly affected by the increase in rk ,particularly through the changes in the flow strength
not only in uppermost part but also in the lowermost
part of fluid layer. On the one hand, the changes in theflow strength in the uppermost part is characterized by
the thinning of stagnant lids along the top surface by theincrease in rk . This is because for larger rk the temperature T of the reference state is higher (see also Fig. 1a)and, hence, the regions of high viscosity η becomethinner owing to the strong temperature dependence inviscosity. On the other hand, the change in the lowermost
part can be seen as an emergence of the regions of insignificant fluid motion just above the bottom boundary
when rk is as large as 100. The emergence of such distinct
regions for very large rk results from the changes in thestatic stability of thermal stratification for the reference
state. As can be seen from Fig. 1c, thick regions of stable
thermal stratification occur above the bottom boundaries
for rk = 100, while no such regions occur for rk ≤ 10.
The insignificant fluid motion at depth is consistent with
our earlier studies based on the linear (Kameyama 2016)
and nonlinear (Kameyama and Yamamoto 2018) analyses
in Cartesian geometry.
The comparison of Fig. 3b and c with 3a also suggeststhat the changes in the values of ℓc0 and Kc0 with increasing rk for large rη in Fig. 2b can be understood from thechanges in the vertical extent of incipient convection
cells. For example, the increase in rk from 1 to 10 makes
the stagnant lids along the top surface thinner and, hence,the convection cells beneath the lids taller. This in turnincreases the horizontal extent of convection cells and, inother words, decreases ℓc0 and Kc0 of infinitesimal perturbations. Increasing rk from 10 to 100, by contrast, notonly leads to thinner stagnant lids along the top surface
but also thickens the regions of insignificant fluid motion
above the bottom boundary. Owing to these combined
effects, the vertical extent of convection cells becomes
smaller as rk increases to 100, leading to an increase in ℓc0
and Kc0 of infinitesimal perturbations.
In Fig. 4, we show the variations in the radial profilesof eigenfunctions for radial velocity V for the criticalstates as functions of rη obtained for the various model
geometries and several values of 1 ≤ rk ≤ 100. Shownin colors and by thin contour lines in each plot of Fig. 4
are the values of V normalized by their maxima Vmax for
given rη, while indicated by the green lines are the heights
at which V = Vmax for given rη. In the cases where theregions of stable thermal stratification occur above thebottom boundary, we show by the thick gray lines theheights below which the thermal stratification is stable.
Fig. 4 shows that a cold and stagnant lid develops along
the top surface for sufficiently large rη, regardless of the(See figure on next page.)
Fig. 3 Schematic illustrations on two-dimensional planes of the structures of incipient flows obtained for the cases with strong
temperature-dependent viscosity rη = 108 and various model geometries together with the values of a rk = 1, b rk = 10 and c rk = 100. Notethat for the cases with 3D spherical geometry we showed the flow structures in meridian planes assuming the incipient flows occurring with only
zonal components (i.e., the spherical harmonic order m = 0). Shown in color are the distributions of infinitesimal perturbations in temperature T ′normalized by their maxima max. For the color scale for the distribution of T /�max, see the scale bar at the bottom of the figure. Shown in solidlines are the contour lines of a “potential” ψ of the mass flux ρv ′ . The contour lines are drawn for 0.1 ≤ |ψ/ψmax | ≤ 0.9 with the interval of 0.1, whereψmax is the maximum of |ψ|Kameyama ﻿Earth, Planets and Space

Fig. 3 (See legend on previous page.)(2021) 73:167

Page 10 of 17

Kameyama ﻿Earth, Planets and Space(2021) 73:167

Page 11 of 17

Fig. 4 The radial profiles of the velocity perturbations V against rη in the fluid layer of different model geometries obtained for the cases with ark = 1, b rk = 3, c rk = 10, d rk = 30 and e rk = 100. In each diagram, the vertical axis indicates the dimensionless height, while the horizontal
axis indicates log10 rη. Shown in color and thin solid contours are the amplitudes of radial velocity perturbations V normalized by their maxima
Vmax for given rη. For the color scale for the distribution of V /Vmax, see the scale bar at the bottom of the figure. The contour lines are drawn for
0.1 ≤ |V /Vmax | ≤ 0.9 with the interval of 0.1. The green lines indicate the heights which yield V = Vmax for given rη. In some cases where regions ofstable thermal stratification exist, we also show by thick gray lines the heights below which the thermal stratification is stablevalues of rk and model geometries. Indeed, all the plots in
Figure 4 show that the regions with very small V develop
near the top surface when rη is sufficiently larger than
the threshold around O ­(103–104). In addition, for thecases with rk = 100 in 3D and 2D spherical geometries,
the regions with very small V develop near the bottom
boundary when rη is sufficiently large. The occurrence ofsuch regions reflects the emergence of thick regions ofstable thermal stratification above the bottom boundaries, and is consistent with the results of the Cartesiancases where such regions develop for rk ≥ 30.
Moreover, from a detailed comparison of the locationof the maxima in vertical velocity V between the caseswith different rk in Fig. 4, we can see that the static stability of thermal stratification in the fluid layer exerts strong
controls on the structures of incipient flows. For the cases
where regions of stable thermal stratification are absent
in the fluid layer (rk ≤ 10 for 3D and 2D spherical geometries, and rk ≤ 3 for Cartesian geometry), the locationsof the maxima in V shift downward with increasing rη,which approach to the bottom boundary for asymptotically large rη. For the other cases where regions of stablestratification exist, by contrast, the locations of the maxima in V, which shift downward with rη for small rη, do
not reach the bottom boundary even when rη is asymptotically large. Instead, the asymptotic height of theKameyama ﻿Earth, Planets and Space(2021) 73:167

Page 12 of 17

Fig. 5 An overall structure of the stagnant-lid convection of Ra = 109, rk ≡ kbot /ktop = 30 and rη ≡ ηtop /ηbot = 107 obtained in the numericalexperiments of Kameyama and Yamamoto (2018). a The distributions of temperature T. The contour interval is 0.05. The color scale is indicatedat the bottom of the figure. b Plots against dimensionless height z of the horizontally averaged temperature T  (thick red line) at the height ofz. Also shown for comparison by thin black lines are the plots of several adiabats against z, while by dotted black lines are the vertical profiles oftemperature T given by steady one-dimensional heat conduction in the vertical direction. c Plots against z of the magnitude of vertical gradient ofT  (red solid line) at the height of z. Also shown for comparison by blue dashed line is the magnitude of adiabatic temperature gradients given by
T  at z. d The distributions of potential temperature Tpot. The color scale is indicated at the bottom of the figure. The contour lines for Tpot are shownat the interval of 0.005 only in the range of Tpot ≥ 0.49. e Plots against z of the horizontally
averaged potential temperature Tpot at the height of z. fPlots against z of the root-mean-squares of the magnitude of the velocity vector v 2  (red solid line) and the vertical velocity vz2  (blue dashedline) at the height of z

maxima in V for very large rη is equal to that at the top
surface of the stable stratification indicated by the thickgray lines in the figure. In short, the incipient convection
cells tends to shrink in the radial direction around theheight of the top surface of the stable stratification withincreasing rη, in response to the thickening of stagnant
lids along the top surface and that of regions of insignificant fluid motion above the bottom boundary.
We will conclude this section by demonstrating thesimilarity between the results of linear analysis and thoseof nonlinear one. In Fig. 5 we show the overall structureof stagnant-lid convection with large rk (= 30) obtainedin our earlier numerical experiments (Kameyama andYamamoto 2018). Because of the large increase in thermal
conductivity k with depth, the temperature gradient in
the vertical direction is greatly reduced at depth (Fig. 5a,b), and it eventually becomes smaller than the adiabatic
gradient near the bottom boundary (Fig. 5c). This leads toan emergence of a region of stable thermal stratification

at the bottom of the layer (Fig. 5d, e), although its dimensionless thickness (∼ 0.25) is smaller than that expectedfrom the static stability of the reference conductive state(∼ 0.37; Fig. 1c). In this region, the fluid motion is indeed
insignificant particularly in the vertical direction (Fig. 5f ).
From this figure, we can confirm that the insignificant
fluid motion in the basal region is largely due to the thestable thermal stratification there.

Conclusion
In this paper, we carried out a series of linear analyses
on the onset of thermal convection of highly compressible fluids whose physical properties strongly vary inspace in convecting vessels of either a three-dimensional(3D) spherical shell or a two-dimensional (2D) sphericalannulus, by extending our earlier work in the Cartesiangeometry (Kameyama 2016): the variations in thermodynamic properties (thermal expansivity and reference
density) with depth are taken to be relevant for theKameyama ﻿Earth, Planets and Space(2021) 73:167

super-Earths with ten times the Earth’s mass, while thethermal conductivity and viscosity are assumed to exponentially depend on depth and temperature, respectively.
One of the major finding in our analyses is that, for thecases with strong temperature dependence in viscosity and strong depth dependence in thermal conductivity, the critical Rayleigh number is on the order of 108
to 109, which is close to the values of Rayleigh numberof the thermal convection in the mantle of super-Earths
with ten times the Earth’s mass extrapolated from that in
the Earth’s mantle. This suggests that the mantle convection of massive super-Earths is most likely to fall in thestagnant-lid regime very close to the critical condition, if
the properties of their mantle materials are quite similar
to the Earth’s.
Our analyses also demonstrated that the structuresof incipient flows of stagnant-lid convection in thepresence of strong adiabatic compression are strongly
affected by the depth dependence in thermal conductivity and the geometries of convecting vessels, throughthe changes in the static stability of thermal stratification of the reference (conductive) state. For example,
the increase in thermal conductivity with depth tends
to enhance the stability at depth in the fluid layer, by
raising the temperature and reducing the gradient in
temperature in the radial direction at depth. Whenthe increase is sufficiently large ( rk ≡ kbot /ktop ≃ 100 ),regions of stable thermal stratification can occur in the
deep mantles of super-Earths with a mass ten times
that of the Earth. The occurrence of stable thermal
stratification in the deep mantle strongly prevents theincipient cells of stagnant-lid convection from extending down to the bottom boundary, leading to the formation of “deep stratospheres” of insignificant fluid
motions above the bottom hot boundaries in addition
to the stagnant lids along the top cold surfaces. The
occurrence of such motionless regions above the bottom boundaries is consistent with our earlier fully
dynamic numerical experiments in Cartesian geometry(Kameyama and Yamamoto 2018). On the other hand,the spherical geometries of convecting vessels tend to
suppress the occurrence of “deep stratospheres” above
the bottom boundaries. This is because the reductionof surface area with depth in the spherical geometries
steepen the gradient in temperature in the radial direction at depth. Furthermore, the increase in thermal conductivity k may not necessarily be so large as assumed
here in silicate mantles of real super-Earths (e.g., Dalton et al. 2013; Hsieh et al. 2017, 2018; Deschamps
and Hsieh 2019), and the temperature dependence inPage 13 of 17

k is likely to alter the nature of stagnant-lid convection(Deschamps 2021). Our analyses showed, however, thatthe regions of insignificant fluid motion can be formed
above the bottom boundaries even in the three-dimensional spherical geometry in the presence of a factor
of about 100 increase in the thermal conductivity withdepth. We can, therefore, speculate that the stagnantlid convection in the mantles of massive super-Earths
is accompanied by another motionless regions at thebase of the mantles if the thermal conductivity of theirmantle materials strongly increases with depth (orpressure).
The occurrence of stable thermal stratification in the
mantles of massive super-Earths may strongly influence the mantle dynamics and surface environmentsof the particular terrestrial planets. For example, thestable region is most likely to greatly suppress the hot
ascending plumes originating from the core–mantle
boundary, resulting in the lack of hotspot volcanism
on super-Earths. In addition, the stable region is most
likely to slow down the cooling of underlying metallic
cores, which may further imply the absence of intrinsic
magnetic fields of super-Earths.
Our findings on the structures of incipient convection of highly compressible fluids may have strong
impact on the studies on the evolution of massive terrestrial planets, in particular on the theoretical ones
based on a parameterized convection (e.g. Tachinami
et al. 2011; Foley 2019). In these earlier studies, theefficiency of convective heat transfer and emission ofmasses (such as greenhouse gases) from the planetary
mantles are estimated, under the assumption that theconvection structures in the mantles of massive superEarths are similar to those of smaller terrestrial planets such as the Earth’s. However, the occurrence of“deep stratospheres” presented in this study indicatesthat this assumption can be significantly violated in
the mantles of massive terrestrial planets owing to theeffect of adiabatic compression, in particular when thetransport properties (i.e., viscosity and thermal conductivity) of mantle materials strongly vary in space.
In addition, from the estimate of Rayleigh number andcomparison with the critical value, the mantle convection of super-Earths is expected to have flow structuresquite similar to the incipient ones presented here. We,
therefore, conclude that the development of convection
models which appropriately take into account the effectof strong adiabatic compression is of ultimate importance in deepening the understanding of the evolutionof massive super-Earths.

Kameyama ﻿Earth, Planets and Space(2021) 73:167

Appendix
Detail of linear analysis in two‑ and three‑dimensional
spherical geometries

According to a general course of linear analysis (e.g.,
Chandrasekhar 1961; Kameyama et al. 2013, 2015;
Kameyama 2016), the evolutionary equations for perturbations are solved by separation of variables into temporal, radial, and horizontal directions. Suppose that thetemporal and spatial dependence of perturbations intemperature T ′ and radial velocity vr′ is given by,
 ′
 

T (r, θ , φ, t)�(r) lm=
Y (θ, φ) exp (at)(21)
vr′ (r, θ , φ, t)
V (r)Page 14 of 17

for three-dimensional (3D) spherical geometry or
 ′
 

T (r, θ , t)�(r)=exp (iK θ + at)vr′ (r, θ , t)
V (r)

(22)

for two-dimensional (2D) spherical cases. Here a is thegrowth rate, Y lm for 3D cases is the spherical harmonics
of the degree ℓ and order m, and K for 2D cases is thewavenumber of perturbation in the horizontal (θ -) direction. By substituting the above equation into (12), (13),and (14), and further eliminating pressure p′ and assuming a = 0, we obtain the equations for the radial profilesof the perturbations in temperature  and radial velocity
V for the critical state of thermal convection. The equations for 3D spherical geometry are1
d4 V
Rac L2 2 ρ α� = η 4
r
dr
�
�
� �
dη d3 V
1 d ln ρ+2+ η 8 +rdr
dr dr 3
�
�
� �
�
�
2
2
�
�1
ln
ρ
ρ
ρ
η
dη
d2 V
d
ln
d
d
ln
1
d
1
+
+3
+
2
10
+
+ η − 2L2 − 12 2 + 6
r
r dr
dr 2
dr
r
dr
dr 2 dr 2

 �
�
�
�
3
2
 dV
 η −4L2 13 − L2 − 6 12 d ln ρ + 12 1 d ln2 ρ + 3 d ln3 ρr dr �
dr
r
r
dr
�
�
�
+
�
�
2
2
 + dη − 2L2 − 6 12 + 6 1 d ln ρ + 4 d ln2 ρ + d η2 2 1 + d ln ρ  dr
r dr
r
dr
dr
r
dr
dr
 � �
�
�
�
�
3
2
4
d
ln
d
ln
ρ
ρ
d

η L2 L2 − 2 r14 − 2L2 r13 dr − L2 − 6 r12 dr 2 + 6 1r drln3 ρ + d drln4 ρ 






� �
�
� 1 d ln ρdηd3 ln ρ1 d2 ln ρ2V,+ + dr − 2L − 2 r 2 dr + 6 r dr 2 + 2 dr 3


��
�


�
2
2


η
d
d
ln
ρ
+

L2 − 2 r12 + dr 2
dr 2

(23)

Kameyama ﻿Earth, Planets and Spaceρ(2021) 73:167

Page 15 of 17




dT
d2 �
2 dk d�1+ Di α T V = k 2 + k +− L2 k 2 �,
dr
dr
r
dr dr
r

(24)
where
The equations for 2D sphericalgeometry are, on the other hand,L2 ≡ ℓ(ℓ + 1).

temperature dependence in η ), we seek for a critical
Rayleigh number Rac of perturbations characterized byℓ or K. In this procedure, we first recast the set of Eqs.
(23) and (24) [or that of (25) and (26)] into an eigenequation of the form

1
d4 V
Rac K 2 2 ρ α � = η 4
r
dr
�
�
� �
dη d3 V
1 d ln ρ+2+ η 6 +rdr
dr dr 3
�
�
� �
�
�
�
�1
dη
d2 ln ρd ln ρ
1d 2 η d2 V
1 d ln ρ
2
+
+3
7 +2
+ 2
+ η − 2K − 5 2 + 5
r
r dr
dr 2
dr
r
dr
dr
dr 2

 � �
�
�
�
�
3
2
ln ρ

η − 2K 2 + 1 r13 − K 2 − 3 r12 d dr
+ 10 1r d drln2 ρ + 3 d drln3 ρ 




 dV

� �
�
�1
dη
d2 ln ρ1 d ln ρ2+ + dr − 2K − 1 r 2 + 5 r dr + 4 dr 2


dr
�
�





 + d2 η 1 + d ln ρ
2
r
dr
dr
 ��
�
�2 1
�
� 2
4
3
ln ρ
2

− K 2 − 3 r12 d drln2 ρ + 5 1r d drln3 ρ + d drln4 ρ 
η K − 1 r 4 − 2K 2 r13 d dr






�
�
�1
dη � 2d3 ln ρ1 d2 ln ρ
1 d ln ρ2+ + dr K − 1 r 3 − 2K r 2 dr + 5 r dr 2 + 2 dr 3
V,


��
�


�1
2
2


d
η
d
ln
ρ
2

+
K − 1 r 2 + dr 2
dr 2



dT
d2 �
1 dk d�1+ Diα T V = k 2 + k +− K 2 k 2 �.
ρ
dr
dr
r
dr dr
r


(26)
In (23) and (25) we rewrote Ra with Rac.
The boundary conditions for  are=0

at r = rmax , rmin ,(27)

where r = rmax and r = rmin denote the top and bottom
boundaries in nondimensional units, respectively. On theother hand, the boundary conditions for V are
1 d ln ρ dV
d2 V
+ 2 +
0=V =
dr 2
r
dr
dr

for 3D spherical geometry, and
1 d ln ρ dV
d2 V
+
+
0=V =
dr 2
r
dr
dr

at r = rmax , rmin

(28)

at r = rmax , rmin(29)

for 2D spherical geometry.
We study the variation in the nature of the infinitesimal perturbations which become unstable for thesmallest Rayleigh number by a following strategy.
First, for given conditions (depth dependence in k and

A·x =

1
x
Rac

(25)

(30)

and then solved for an eigenvalue (1/Rac ) together witheigenfunctions x (radial profiles of perturbations of temperature  and radial velocity V) describing a criticalstate (i.e., growth rate a = 0). The numerical details on
how to construct the matrix A can be found in Kameyama et al. (2015). The eigensystem developed thus was
solved for all eigenvalues and eigenvectors by invoking
the LAPACK subroutine “dgeev” designed for real nonsymmetric matrices. Among the series of eigenvalues andeigenfunctions in hand, we then search the ones whichyield the absolute minimum of the critical Rayleigh number Rac0 and the degree ℓc0 or wavenumber Kc0 corresponding to Rac0.
To conclude this appendix, we note that a special caremust be taken in the final stage of the solution of the present eigensystem. Because the matrix A is in general a nonsymmetric matrix, the eigenequation (30) has to be solvedin a field of complex numbers, even if the matrix A is real(e.g., Zebib et al. 1983; Zebib 1993). In fact, we obtained aset of complex eigenvalues and their conjugation. In addition, the distributions of eigenvalues of A significantly
vary depending on whether a region of stable stratification

Kameyama ﻿Earth, Planets and Space(2021) 73:167

occurs or not. Indeed, all the eigenvalues have positive realparts (i.e. A is positive definite) when a thermal stratification is unstable in the entire layer. However, when a regionof stable stratification occurs in the fluid layer, the eigenvalues have both positive and negative real parts (i.e. A is
indefinite). ( Kameyama et al. (2015) further showed that
all the eigenvalues have negative real parts (i.e. A is negative definite) when a thermal stratification is stable in theentire layer.) Since the eigenvalues of A whose real parts
are negative are physically unsound (see Kameyama et al.
(2015) for detail), we can only focus on the eigenvalues
1/Rac with positive real parts.
Acknowledgements
We thank Masaki Ogawa, Takehiro Miyagoshi, Takatoshi Yanagisawa andHaruhiko Dekura for fruitful discussion.
Authors’ contributions
MK developed the numerical codes, performed numerical simulations,analyzed simulation data, and prepared the manuscript. The author read andapproved the final manuscript.
Funding
This work was supported by JSPS KAKENHI Grant Number 18K03724, a
joint research program at Geodynamics Research Center, Ehime University,“Exploratory Challenge on Post-K computer” (Elucidation of the Birth of Exoplanets [Second Earth] and the Environmental Variations of Planets in the Solar
System), and “Program for Promoting Research on the Supercomputer Fugaku”(Toward a unified view of the universe: from large scale structures to planets).
Availability of data and materials
The data that support the findings of this study and the numerical programs
used in this study are available upon reasonable request to M.K. (email: kameyama@sci.ehime-u.ac.jp).

Declarations
Competing interests
The author declares no competing interests.
Author details
1
Geodynamics Research Center (GRC), Ehime University, Matsuyama, Ehime
790‑8577, Japan. 2 Japan Agency for Marine-Earth Science and Technology(JAMSTEC), Yokohama, Kanagawa 236‑0001, Japan.
Received: 13 May 2021 Accepted: 15 August 2021

References
Chandrasekhar S (1961) Hydrodynamic and hydromagnetic stability. The
Clarendon Press, Oxford, p 652
Čížková H, van den Berg A, Jacobs M (2017) Impact of compressibility on heattransport characteristics of large terrestrial planets. Phys Earth Planet Int
268:65–77. https://​doi.​org/​10.​1016/j.​pepi.​2017.​04.​007
Dalton DA, Hsieh WP, Hohensee GT, Cahill DG, Goncharov AF (2013) Effect ofmass disorder on the lattice thermal conductivity of MgO periclase underpressure. Sci Rep. https://​doi.​org/​10.​1038/​srep0​2400
Dekura H, Tsuchiya T (2017) Ab initio lattice thermal conductivity of MgO froma complete solution of the linearized Boltzmann transport equation. Phys
Rev B 95:184303. https://​doi.​org/​10.​1103/​PhysR​evB.​95.​184303
Dekura H, Tsuchiya T (2019) Lattice thermal conductivity of MgSiO3 postperovskite under the lowermost mantle conditions from ab initio anharmonic lattice dynamics. Geophys Res Lett 46(22):12919–12926. https://​doi.​org/​10.​1029/​2019G​L0852​73

Page 16 of 17

Dekura H, Tsuchiya T, Tsuchiya J (2013) Ab initio lattice thermal conductivity of MgSiO3 perovskite as found in Earth’s lower mantle. Phys Rev Lett
110(2):025904. https://​doi.​org/​10.​1103/​PhysR​evLett.​110.​025904
Deschamps F (2021) Stagnant lid convection with temperature-dependent
thermal conductivity and the thermal evolution of icy worlds. Geophys J
Int 225(3):1870–1889. https://​doi.​org/​10.​1093/​gji/​ggaa5​40
Deschamps F, Hsieh WP (2019) Lowermost mantle thermal conductivity constrained from experimental data and tomographic models. Geophys J Int
219(Supplement 1):S115–S136. https://​doi.​org/​10.​1093/​gji/​ggz231
Foley BJ (2019) Habitability of Earth-like stagnant lid planets: climate evolution
and recovery from snowball states. Astrophys J 875(1):72. https://​doi.​org/​
10.​3847/​1538-​4357/​ab0f31
Gaidos E, Conrad CP, Manga M, Hernlund J (2010) Thermodynamic limits onmagnetodynamos in rocky exoplanets. Astrophys J 718(2):596–609.
https://​doi.​org/​10.​1088/​0004-​637X/​718/2/​596
Guerrero JM, Lowman JP, Deschamps F, Tackley PJ (2018) The influence of
curvature on convection in a temperature-dependent viscosity fluid:
implications for the 2-D and 3-D modeling of Moons. J Geophys Res
Planets 123(7):1863–1880. https://​doi.​org/​10.​1029/​2017J​E0054​97
Haigis V, Salanne M, Jahn S (2012) Thermal conductivity of MgO, MgSiO3 perovskite and post-perovskite in the Earth’s deep mantle. Earth Planet Sci
Lett 355–356:102–108. https://​doi.​org/​10.​1016/j.​epsl.​2012.​09.​002
Howard AW, Marcy GW, Johnson JA, Fischer DA, Wright JT, Isaacson H,
Valenti JA, Anderson J, Lin DNC, Ida S (2010) The occurrence and mass
distribution of close-in super-Earths, Neptunes, and Jupiters. Science33(6004):653–655. https://​doi.​org/​10.​1126/​scien​ce.​11948​54
Hsieh WP, Deschamps F, Okuchi T, Lin JF (2017) Reduced lattice thermal conductivity of Fe-bearing bridgmanite in Earth’s deep mantle. J Geophys
Res Solid Earth 122(7):4900–4917. https://​doi.​org/​10.​1002/​2017J​B0143​39
Hsieh WP, Deschamps F, Okuchi T, Lin JF (2018) Effects of iron on the lattice
thermal conductivity of Earth’s deep mantle and implications for mantle
dynamics. Proc Natl Acad Sci 115(16):4099–4104. https://​doi.​org/​10.​1073/​
pnas.​17185​57115
Jarvis GT, McKenzie DP (1980) Convection in a compressible fluid with infinite
Prandtl number. J Fluid Mech 96(3):515–583. https://​doi.​org/​10.​1017/​
S0022​11208​00022​5X
Kameyama M (2016) Linear analysis on the onset of thermal convection ofhighly compressible fluids with variable physical properties: implications
for the mantle convection of super-Earths. Geophys J Int 204(2):1164–
1178. https://​doi.​org/​10.​1093/​gji/​ggv507
Kameyama M, Kinoshita Y (2013) On the stability of thermal stratification ofhighly compressible fluids with depth-dependent physical properties:
implications for the mantle convection of super-Earths. Geophys J Int
195(3):1443–1454. https://​doi.​org/​10.​1093/​gji/​ggt321
Kameyama M, Ogawa M (2000) Transitions in thermal convection withstrongly temperature-dependent viscosity in a wide box. Earth Planet Sci
Lett 180(3–4):355–367. https://​doi.​org/​10.​1016/​S0012-​821X(00)​00171-0
Kameyama M, Yamamoto M (2018) Numerical experiments on thermal convection of highly compressible fluids with variable viscosity and thermal
conductivity: implications for mantle convection of super-Earths. Phys
Earth Planet Int 274:23–36. https://​doi.​org/​10.​1016/j.​pepi.​2017.​11.​001
Kameyama M, Ichikawa H, Miyauchi A (2013) A linear stability analysis on
the onset of thermal convection of a fluid with strongly temperaturedependent viscosity in a spherical shell. Theor Comput Fluid Dyn27(1–2):21–40. https://​doi.​org/​10.​1007/​s00162-​011-​0250-x
Kameyama M, Miyagoshi T, Ogawa M (2015) Linear analysis on the onset ofthermal convection of highly compressible fluids: implications for the
mantle convection of super-Earths. Geophys J Int 200(2):1066–1077.
https://​doi.​org/​10.​1093/​gji/​ggu457
Karato S (2011) Rheological structure of the mantle of a super-Earth: someinsights from mineral physics. Icarus 212(1):14–23. https://​doi.​org/​10.​1016/j.​icarus.​2010.​12.​005
King SD, Lee C, van Keken PE, Leng W, Zhong S, Tan E, Tosi N, Kameyama M(2010) A community benchmark for 2-D Cartesian compressible convection in the Earth’s mantle. Geophys J Int 180(1):73–87. https://​doi.​org/​10.​1111/j.​1365-​246X.​2009.​04413.x
Kite ES, Manga M, Gaidos E (2009) Geodynamics and rate of volcanism on
massive Earth-like planets. Astrophys J 700(2):1732–1749. https://​doi.​org/​
10.​1088/​0004-​637X/​700/2/​1732

Kameyama ﻿Earth, Planets and Space(2021) 73:167

Korenaga J (2010) On the likelihood of plate tectonics on super-earths: does
size matter? Astrophys J 725(1):L43–L46. https://​doi.​org/​10.​1088/​2041-​
8205/​725/1/​L43
Liu X, Zhong S (2013) Analyses of marginal stability, heat transfer andboundary layer properties for thermal convection in a compressible fluid with infinite Prandtl number. Geophys J Int 194(1):125–144.
https://​doi.​org/​10.​1093/​gji/​ggt117
Lobanov SS, Holtgrewe N, Lind JF, Goncharov AF (2017) Radiative conductivity and abundance of post-perovskite in the lowermost mantle.
Earth Planet Sci Lett 479:43–49. https://​doi.​org/​10.​1016/j.​epsl.​2017.​09.​016
Miyagoshi T, Tachinami C, Kameyama M, Ogawa M (2014) On the vigor ofmantle convection in super-Earths. Astrophys J. https://​doi.​org/​10.​1088/​2041-​8205/​780/1/​L8
Miyagoshi T, Kameyama M, Ogawa M (2015) Thermal convection and theconvective regime diagram in super-Earths. J Geophys Res Planets
120(7):1267–1278. https://​doi.​org/​10.​1002/​2015J​E0047​93
Miyagoshi T, Kameyama M, Ogawa M (2018) Effects of adiabatic compression on thermal convection in super-earths of various sizes. Earth
Planets Space 70:200. https://​doi.​org/​10.​1186/​s40623-​018-​0975-5
Miyauchi A, Kameyama M (2013) Influences of the depth-dependence ofthermal conductivity and expansivity on thermal convection withtemperature-dependent viscosity. Phys Earth Planet Int 223:86–95.
https://​doi.​org/​10.​1016/j.​pepi.​2013.​08.​001
Noack L, Breuer D (2014) Plate tectonics on rocky exoplanets: influenceof initial conditions and mantle rheology. Planet Space Sci 98:41–49.
https://​doi.​org/​10.​1016/j.​pss.​2013.​06.​020
Noack L, Rivoldini A, Van Hoolst T (2017) Volcanism and outgassing of
stagnant-lid planets: implications for the habitable zone. Phys Earth
Planet Int 269:40–57. https://​doi.​org/​10.​1016/j.​pepi.​2017.​05.​010
Ohta K, Yagi T, Taketoshi N, Hirose K, Komabayashi T, Baba T, Ohishi Y,
Hernlund J (2012) Lattice thermal conductivity of MgSiO3 perovskite
and post-perovskite at the core-mantle boundary. Earth Planet Sci Lett
349–350:109–115. https://​doi.​org/​10.​1016/j.​epsl.​2012.​06.​043
Ohta K, Yagi T, Hirose K, Ohishi Y (2017) Thermal conductivity of ferropericlase in the Earth’s lower mantle. Earth Planet Sci Lett 465:29–37.
https://​doi.​org/​10.​1016/j.​epsl.​2017.​02.​030
Papuc AM, Davies GF (2008) The internal activity and thermal evolution ofEarth-like planets. Icarus 195(1):447–458. https://​doi.​org/​10.​1016/j.​
icarus.​2007.​12.​016
Rivera EJ, Lissauer JJ, Butler RP, Marcy GW, Vogt SS, Fischer DA, Brown TM,
Laughlin G, Henry GW (2005) A sim7.5 extMo plus planet orbiting the
nearby star, GJ 876. Astrophys J 634(1):625–640. https://​doi.​org/​10.​1086/​491669
Seager S, Kuchner M, Hier-Majumder CA, Militzer B (2007) Mass-radius relationships for solid exoplanets. Astrophys J 669(2):1279–1297. https://​doi.​org/​10.​1086/​521346
Schubert G, Turcotte DL, Olson P (2001) Mantle convection in the Earth andPlanets. Cambridge University Press, Cambridge, p 940
Shahnas MH, Pysklywec RN, Yuen DA (2018) Penetrative convection insuper-Earth planets: consequences of MgSiO3 postperovskite dissociation transition and implications for super-Earth GJ 876 d. J Geophys
Res Planets 123(8):2162–2177. https://​doi.​org/​10.​1029/​2018J​E0056​33
Sotin C, Grasset O, Mocquet A (2007) Mass-radius curve for extrasolar Earthlike planets and ocean planets. Icarus 191(1):337–351. https://​doi.​org/​
10.​1016/j.​icarus.​2007.​04.​006
Stamenković V, Breuer D, Spohn T (2011) Thermal and transport propertiesof mantle rock at high pressure: applications to super-Earths. Icarus
216(2):572–596. https://​doi.​org/​10.​1016/j.​icarus.​2011.​09.​030
Stein C, Finnenkötter A, Lowman JP, Hansen U (2011) The pressure-weakening effect in super-Earths: consequences of a decrease in lower mantle

Page 17 of 17

viscosity on surface dynamics. Geophys Res Lett. https://​doi.​org/​10.​1029/​2011G​L0493​41
Stengel KC, Oliver DS, Booker JR (1982) Onset of convection in a variableviscosity fluid. J Fluid Mech 120:411–431. https://​doi.​org/​10.​1017/​
S0022​11208​20028​21
Tachinami C, Senshu H, Ida S (2011) Thermal evolution and lifetime of
intrinsic magnetic fields of super-earths in habitable zones. Astrophys J
726(2):70. https://​doi.​org/​10.​1088/​0004-​637X/​726/2/​70
Tachinami C, Ogawa M, Kameyama M (2014) Thermal convection of compressible fluid in the mantle of super-Earths. Icarus 231:377–384. https://​doi.​
org/​10.​1016/j.​icarus.​2013.​12.​022
Tackley PJ, Ammann M, Brodholt JP, Dobson DP, Valencia D (2013) Mantle
dynamics in super-Earths: post-perovskite rheology and self-regulation ofviscosity. Icarus 225(1):50–61. https://​doi.​org/​10.​1016/j.​icarus.​2013.​03.​013
Tosi N, Godolt M, Stracke B, Ruedas T, Grenfell JL, Höning D, Nikolaou A, Plesa
AC, Da Breuer, Spohn T (2017) The habitability of a stagnant-lid Earth.
Astoron Astrophys 605:A71. https://​doi.​org/​10.​1051/​0004-​6361/​20173​
0728
Tsuchiya T, Tsuchiya J (2011) Prediction of a hexagonal SiO2 phase affecting
stabilities of MgSiO3 and CaSiO3 at multimegabar pressures. Proc Nat
Acad Sci USA 108(4):1252–1255. https://​doi.​org/​10.​1073/​pnas.​10135​
94108
Valencia D, O’Connell RJ, Sasselov DD (2006) Internal structure of massive
terrestrial planets. Icarus 181(2):545–554. https://​doi.​org/​10.​1016/j.​icarus.​
2005.​11.​021
Valencia D, O’Connell RJ, Sasselov DD (2007a) Inevitability of plate tectonics onsuper-earths. Astrophys J 670(1):L45–L48. https://​doi.​org/​10.​1086/​524012
Valencia D, Sasselov DD, O’Connell RJ (2007b) Radius and structure models of
the first super-earth planet. Astrophys J 656(1):545–551. https://​doi.​org/​
10.​1086/​509800
van den Berg AP, Yuen DA, Umemoto K, Jacobs MHG, Wentzcovitch RM (2018)
Mass-dependent dynamics of terrestrial exoplanets using ab initio
mineral properties. Icarus 317:412–426. https://​doi.​org/​10.​1016/j.​icarus.​
2018.​08.​016
Wagner FW, Sohl F, Hussmann H, Grott M, Rauer H (2011) Interior structure
models of solid exoplanets using material laws in the infinite pressure
limit. Icarus 214(2):366–376. https://​doi.​org/​10.​1016/j.​icarus.​2011.​05.​027
Wagner FW, Tosi N, Sohl F, Rauer H, Spohn T (2012) Rocky super-Earth interiors:
structure and internal dynamics of CoRoT-7b and Kepler-10b. AstronAstrophys 541:A103. https://​doi.​org/​10.​1051/​0004-​6361/​20111​8441
Yanagisawa T, Kameyama M, Ogawa M (2016) Numerical studies on convective stability and flow pattern in three-dimensional spherical mantle ofterrestrial planets. Geophys J Int 206(3):1526–1538. https://​doi.​org/​10.​1093/​gji/​ggw226
Yao C, Deschamps F, Lowman JP, Sanchez-Valle C, Tackley PJ (2014) Stagnant
lid convection in bottom-heated thin 3-D spherical shells: influence of
curvature and implications for dwarf planets and icy moons. J Geophys
Res Planets 119(8):1895–1913. https://​doi.​org/​10.​1002/​2014J​E0046​53
Zebib A (1993) Linear and weakly nonlinear variable viscosity convection in
spherical shells. Theor Comput Fluid Dyn 4(5):241–253. https://​doi.​org/​10.​
1007/​BF004​17930
Zebib A, Schubert G, Dein JL, Paliwal RC (1983) Character and stability ofaxisymmetric thermal convection in spheres and spherical shells.
Geophys Astrophys Fluid Dyn 23(1):1–42. https://​doi.​org/​10.​1080/​03091​
92830​82090​38

Publisher’s NoteSpringer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.