Geophys. J. Int. (2022) 231, 1457–1469
Advance Access publication 2022 July 9
GJI Geodynamics and Tectonics

https://doi.org/10.1093/gji/ggac259

Numerical experiments on thermal convection of highly compressible
fluids with variable viscosity and thermal conductivity in 2-D
cylindrical geometry: implications for mantle convection ofsuper-Earths
1,2

1 Geodynamics Research Center (GRC), Ehime University, Matsuyama, Ehime 790–8577, Japan. E-mail: kameyama@sci.ehime-u.ac.jp
2 Japan Agency for Marine-Earth and Technology (JAMSTEC), Yokohama, Kanagawa 236–0001, Japan

Accepted 2022 July 7. Received 2022 July 1; in original form 2021 December 3

SUMMARY
We conduct a series of numerical experiments of thermal convection of highly compressible
fluids in 2-D cylindrical annulus, in order to study the mantle convection on super-Earths.
The variations in thermodynamic properties (thermal expansivity and reference density) withdepth are taken to be relevant for the super-Earths with 10 times the Earth’s mass, while thosein transport properties (viscosity and thermal conductivity) are modelled by an exponential
dependence on temperature and/or depth. From our experiments we identified a distinct regimeof convecting flow patterns induced by the interplay between the adiabatic temperature change
and the spatial variations in viscosity and thermal conductivity. That is, for the cases withstrong temperature-dependent viscosity and large increase in thermal conductivity with depth,a ‘deep stratosphere’ of stable thermal stratification is formed at the base of the mantle, in
addition to thick stagnant lids at their top surfaces. In the ‘deep stratosphere’, the fluid motion
is insignificant particularly in the vertical direction in spite of smallest viscosity owing to
its strong dependence on temperature. From the comparison with the experiments with the
Cartesian geometry, we also found that the occurrence of ‘deep stratosphere’ tends to be
suppressed for the cases with cylindrical geometry, owing to the reduction of the surface
area with depth which helps increase the temperature gradient in the lowermost mantle. Ourfinding may further imply that both the effects of adiabatic compression and those of spherical(or cylindrical) geometry of mantle are of crucial importance in understanding the mantle
dynamics of massive super-Earths in the presence of spatial variations in physical properties.
Key words: High-pressure behaviour; Numerical modelling; Planetary interiors; Dynamics:
convection currents, and mantle plumes; Heat generation and transport.

1 I N T RO D U C T I O N
Since the first discovery of extra-solar ‘super-Earths’ (Rivera et al.
2005), a growing effort has been devoted to the studies in the internal(radial) structure (e.g. Seager et al. 2007; Valencia et al. 2007b;
Sotin et al. 2007; Wagner et al. 2011), the thermal evolution (e.g.
Papuc & Davies 2008; Kite et al. 2009; Korenaga 2010; Gaidos
et al. 2010; Tachinami et al. 2011; Čı́žková et al. 2017) and thehabitability (e.g. Noack et al. 2017; Tosi et al. 2017; Foley 2019)of such massive terrestrial bodies. During the last decades, thestudies on mantle dynamics have been playing an important role inaddressing these issues (e.g. van den Berg et al. 2018; Shahnas &Pysklywec 2021).
It can be intuitively understood that, because of their large 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 andsmaller terrestrial planets. For example, the variations of physical properties of mantle materials are expected to be larger inthe much thicker 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 theconvecting flow patterns have been widely studied so far, under
the assumption that the effects of other physical mechanisms suchas compressibility can be minor (e.g. Kameyama & Ogawa 2000;
Miyauchi & Kameyama 2013). In particular, the effects of strong
temperature-dependence in viscosity have been intensively studied
by numerical models of convection in 3-D spherical geometry (e.g.
Yanagisawa et al. 2016; Guerrero et al. 2018).
On the other hand, it has been well acknowledged that the effectof adiabatic compression affects the nature of thermal convection


C The Author(s) 2022. Published by Oxford University Press on behalf of The Royal Astronomical Society. This is an Open Accessarticle distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which
permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.

1457

Downloaded from https://academic.oup.com/gji/article/231/2/1457/6638987 by Ehime University Library user on 05 August 2022

Masanori Kameyama

1458

M. Kameyama

2 MODEL DESCRIPTIONS
We consider a time-dependent thermal convection of a highly compressible fluid of infinite Prandtl number, as a model of mantle
convection of a super-Earth with 10 times the Earth’s mass. The
numerical model is almost the same as those of our earlier one(Kameyama & Yamamoto 2018) except that the convecting vessel
is taken to be a half (azimuth θ min ≤ θ ≤ θ max where θ max − θ min= π ) of 2-D cylindrical annulus whose inner and outer radii arermin and rmax , respectively, and rmin /rmax = 0.5. The temperature Tof the layer is fixed to be T = Ttop and Tbot (>Ttop ) at the top andbottom boundaries, respectively, while the vertical side walls are
insulating (zero heat flux). The impermeable and shear-stress-free
conditions are adopted along all boundaries. In order to keep thepresent analysis as simple as possible, we ignored the effects ofinternal heating and spatial variations in specific heat Cp of the fluid
and the gravity g throughout the layer. In addition we fixed Ttop /T= 0.1 (where T ≡ Tbot − Ttop is the temperature difference acrossthe layer) as is expected in the Earth’s mantle, in order to reduce thenumbers of varying parameters. We here admit that the choice ofTtop /T = 0.1 is somewhat inappropriate for the mantle convectionof super-Earths whose surface temperature is expected to be muchhigher than the Earth’s owing to their smaller orbital radii. In otherwords, our results may overestimate the vigour of thermal convection in the presence of adiabatic compression, because of lower Ttop
and colder interiors than the actual ones which in turn reduces theeffects of adiabatic temperature changes (see also Tachinami et al.
2014).
2.1 Physical properties of modelled compressible fluid
As a model of pressure-dependence in thermodynamic propertiesof mantle materials of massive super-Earths, we use the same variations in thermal expansivity α and reference density ρ with depths
as in our earlier works (Miyagoshi et al. 2015; Kameyama & Yamamoto 2018). The thermal expansivity α decreases with depth,
whose radial profile is taken to be similar to that estimated for MgOfrom an ab initio calculation (Tsuchiya & Tsuchiya 2011) for the
mantles of super-Earths with 10 times the Earth’s mass. That is, thevalue of α(r = rmin ) at the bottom boundary is as small as about
1/14 of the value α(r = rmax ) ≡ α top (=4 × 10−5 K−1 ) at the top
surface, and the value of its average α ave is taken to be r =rmax
α(r )
αave
1
=
dr  0.1574.
(1)
αtoprmax − rmin r =rmin αtop
On the other hand, the reference density ρ is assumed to linearly
increase with depth, yielding the value of ρ(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 increase in ρ
assumed here is not likely to spoil the significance of the present
analysis despite its geophysical unsoundness, since the vertical profile of (reference) density has a very 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 in transport properties. In
this study the thermal conductivity k is taken to be exponentiallydependent on depth as,
rmax − r,(2)k(r ) = ktop exp ln(rk ) ×rmax − rmin
where ktop is the thermal conductivity at r = rmax , and rk is a dimensionless constant equal to the ratio of k(r = rmin ) ≡ kbot to ktop . Fora geophysical significance, we restrict ourselves to the cases whererk ≥ 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 & 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 to two orders of magnitude acrossthe mantles of super-Earths. In this study we varied rk up to 100.
On the other hand, the viscosity η is assumed to be exponentiallydependent on temperature T and depth as,Tbot − T− ln(r p )η(T, r ) = ηbot exp ln(r T ) ×
Tbot − Ttopr − rmin,(3)
×rmax − rmin
where ηbot is the viscosity at the bottom boundary (T = Tbot and r= rmin ), and rT and rp are dimensionless constants controlling thetemperature- and depth-dependence in viscosity, respectively. From(3), the ratio rη of viscosity ηtop at the top surface (T = Ttop and r =rmax ) to ηbot can be written as


ηtop
rη ≡
= exp ln(r T ) − ln(r p ) = r T /r p .
(4)
ηbot

Downloaded from https://academic.oup.com/gji/article/231/2/1457/6638987 by Ehime University Library user on 05 August 2022

(Jarvis & McKenzie 1980; King et al. 2010; Liu & Zhong 2013).
This effect is expected to be of great importance in the convection in the thick mantles of massive super-Earths, considering thatthe pressure can reach on the order of terapascals in their interiors(Tsuchiya & Tsuchiya 2011). Indeed, our earlier studies (Tachinami
et al. 2014; Miyagoshi et al. 2015; Kameyama & Yamamoto 2018)
showed that the nature of thermal convection in the thick mantles ofmassive super-Earths is strongly affected by the interplay between
the effects of adiabatic compression and those of spatial variations
in physical properties of mantle materials. However, in these studies the effects of spherical geometries of planetary mantles havebeen thoroughly ignored, since the analysis 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 in physicalproperties of mantle materials on the convecting flow structures in
spherical vessels, in order to gain the fundamental insights into themantle convection of massive super-Earths.
Recently we have carried out a series of linear analysis on thethermal convection in spherical domains of a highly compressible
fluid with the spatial variation in physical properties (Kameyama2021), and demonstrated that the spherical geometries of planetary
mantles can significantly affect the critical state and condition ofthermal convection of a highly compressible fluid. Motivated by thisanalysis, we will perform numerical experiments of thermal convection in a polar coordinate system in the presence of the interplay
between the adiabatic compression and spatial variations in physicalproperties such as viscosity and thermal conductivity. In this paper,
however, we restrict our numerical experiments to the convection
in a 2-D cylindrical annulus, in order to avoid severe computational costs required for numerical solutions with strong variations
in physical properties in large domains. By further comparing thepresent results with our earlier one in the Cartesian (planar) model(Miyagoshi et al. 2015; Kameyama & Yamamoto 2018), we aim at
deepening the insights into the thermal convection in the mantles ofmassive terrestrial super-Earths both from the thermodynamic andfluid-dynamic aspects.

Highly compressible 2D cylindrical convection
In this study, we assume rT ≥ 1 and rp ≥ 1; the viscosity η decreases
with temperature for rT > 1 and decreases with depth for rp > 1.
We set rp = 1 unless otherwise noted. In some cases where rp > 1
is used, we restricted ourselves to the cases with rp ≤ 103 and rη
≥ 105 . By using this assumption, we will exclude the cases with ahuge increase in viscosity in the deep mantle of the super-Earths
proposed in several earlier studies (Stamenković et al. 2011, 2012).
Indeed, a recent theoretical work by Ritterbex et al. (2018) predicts
a substantial decrease in viscosity in the mantles of massive superEarths associated with the B1-B2 transition of MgO at around 500
GPa.

balance is expected to be improved by using better approximations
such as ALA (Leng & Zhong 2008; Liu & Zhong 2013; Ricard et al.
2022).
The present formulation includes two non-dimensional parameters. First one is the dissipation number Di of the modelled fluid
defined by,Di ≡1(∇ ⊗ v + v ⊗ ∇) ,
2

 
1
0 = −∇ p + 2η ε̇ − (∇ · v)I · ∇ + Raρ αT er ,3∂TDi
= ∇ · (k∇T − ρvT ) − Diρ α T vr +2η (ε̇ II )2 ,
ρ
∂t
Ra

ε̇ =

(5)
(6)
(7)
(8)

where the non-dimensionalization is carried out with the lengthscale of d ≡ rmax − rmin , time scale of d2 /κ top (where κ top =ktop /ρ top Cp is the value of thermal diffusivity at the top surface),and temperature scale of T ≡ Tbot − Ttop . Here v is the velocity
vector, ε̇ is the strain rate tensor, ⊗ stands for the tensor product ofvectors, v ⊗ ∇ is the transpose of ∇ ⊗ v, I is the identity tensor, er
is the unit vector in radial direction positive outward, p is pressure,vr ≡ (v · er ) is the velocity in the radial direction, and ε̇ II is the second invariant of strain rate tensor. In the right-hand side of eq. (8),the first, second, and third terms represent the effects of apparent
heat transport (advection plus conduction), adiabatic temperaturechange and viscous dissipation (frictional heating), respectively.
We note here the validity of the truncated anelastic liquid approximation (TALA) in the present numerical experiments. By using
TALA instead of its untruncated counterpart (ALA), we neglected
the effects on buoyancy of dynamic pressure (the deviation of pressure from the hydrostatic one), by assuming that the density of modelled fluids is independent of dynamic pressure (Jarvis & McKenzie1980). Since the formulation based on TALA is known to introduce
an imbalance between viscous and adiabatic heating (Leng & Zhong2008), one may think that it can severely violate the energy budget
in the entire domain particularly under the strong compression characterized by a large dissipation number (King et al. 2010). However,as had been already demonstrated in our earlier work in the Cartesian geometry (Kameyama & Yamamoto 2018), it is not likely thecase for the mantle convection of massive terrestrial planets. For
example, at a statistically steady state of the case with Ra = 109 , rk
= 30 and rη = 104 (rT = 104 and rp = 1), the heat flux at the top surface Qtop defined by (16) is about 7.14 per unit arc length, while thatat the bottom boundary Qbot is about 13.93, indicating that the mismatch between the input (rmin Qbot ) and output (rmax Qtop ) of heat at
the top and bottom boundaries of cylindrical annulus with rmin /rmax= 0.5 is only up to a few per cent. We therefore believe that thepresent modelling based on TALA can be sufficiently valid for theconvecting flows even in the mantles of massive super-Earths since
the energy balance almost holds in the entire domain, although the(9)ρtop αtop T gd 3.
ηbot κtop(10)

Note here that, as indicated by eq. (3), the viscosity is scaled by thevalue ηbot at the bottom boundary while other properties are doneby those at top surface. The values of Di and Ra for the mantles ofterrestrial planets are expected to increase with the mass M of theplanets, coming from the increase in the gravity g and the mantle
thickness d 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∝M0.5 andd∝M0.28 ; (Valencia et al. 2007a)] together with Di ∼ 0.7 and Ra= O(107 ) for the Earth’s mantle (Schubert et al. 2001), one canestimate Di ∼ 5 and Ra = O(108 − 109 ) for the mantles of superEarths whose mass is 10 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 is taken to be 108 ≤ Ra ≤1010 .
The used initial condition for temperature T(r, θ, t = 0) is taken
to be
T (r, θ, t = 0) = Tcond (r ) − δT (r ) cos πθ − θminθmax − θmin

,

(11)

where Tcond (r) represents the steady 1-D conductive state in theradial direction and δT(r) is an amplitude of perturbation which istaken to be O(10−2 ) in the dimensionless unit. The radial profile ofTcond (r) is given bydTcond1 drk, Tcond (r = rmax ) = Ttop ,r drdrTcond (r = rmin ) = Tbot = Ttop + 1,

0=

(12)

and the values of its radial gradient (dTcond /dr) at the top (r =rmax ) and bottom (r = rmin ) boundaries for several rk are shown in
Table 1. We continue the calculation till the average temperaturein the convecting vessel and the heat flow on both the surface andbottom boundaries reach their statistically steady state.

2.3 Numerical method
The basic equations are discretized by the finite-difference method
based on the control-volume method (Patankar 1980). A uniform
mesh is used with 256 and 1024 mesh points in radial and horizontal directions, respectively. We confirmed that the thermal boundary
layer (TBL) along the bottom boundary is spanned by at least threemesh points excluding the boundary point for all the cases where
the bottom TBLs are clearly formed (see below). The energy eq. (8)is discretized by the Crank-Nicholson method in time. An upwind
scheme, called power-law scheme (Patankar 1980), is used to evaluate the advection and conduction terms in the energy equation. The
eqs (5) and (7) for flow field are solved for a ‘potential’ ψ of the

Downloaded from https://academic.oup.com/gji/article/231/2/1457/6638987 by Ehime University Library user on 05 August 2022

2.2 Basic equations

0 = ∇ · (ρv) ,αtop gd.
Cp

Second one is the Rayleigh number Ra of thermal convection defined by,Ra ≡Under the truncated anelastic liquid approximation (TALA), thebasic equations can be written in a dimensionless form as:

1459

1460

M. Kameyama
Table 1. The dimensionless values of the radial gradientof temperature Tcond given by steady-state conduction at
the top (r = rmax ) and bottom (r = rmin ) surfaces of the2-D cylindrical annulus with rmin /rmax = 0.5 and several
values of rk ≡ kbot /ktop controlling the depth-dependence
in thermal conductivity. Also shown for comparison are
the values of dTcond /dz at top (z = 1) and bottom (z = 0)
boundaries of the Cartesian model.
Cylindrical|dTcond /dr|rk

0.72131.26492.08682.99994.11181.44270.84330.41740.20000.0822at z = 1

at z = 0

1.00001.64792.55843.51854.65171.00000.54930.25580.11730.0465

mass flux ρv defined by,
ρvr =

1 ∂ψ
,
r ∂θ

ρvθ = −

∂ψ
,
∂r

(13)

instead of directly solving for primitive variables (velocity v andpressure p). Note that the function ψ is a generalization of thestream function for the Boussinesq or extended Boussinesq approximation where ρ = 1. As in the case with incompressible fluids,
the velocity field v defined above automatically satisfies the equation of continuity (eq. 5). In Appendix we describe in detail the
explicit forms of basic equations and our numerical technique forthe solution of flow field. The reliability of the numerical code
has been already verified (Kameyama & Harada 2017; Tsuchida &Kameyama 2020).

3 R E S U LT S
We carried out numerical experiments of thermal convection ofhighly compressible fluids by systematically varying the depthdependence in thermal conductivity (1 ≤ rk ≤ 100) together withthe temperature-dependence in viscosity (102 ≤ rT ≤ 109 and rp =1) and the Rayleigh number Ra (108 , 109 and 1010 ) as indicated by
the symbols in Fig. 1. The snapshots of the overall thermal convection in 2-D cylindrical annulus for all the cases calculated inthis study can be found in Figs S2–S14. (We also conducted some
additional calculations by incorporating the depth-dependence inviscosity with 100 ≤ rp ≤ 103 and rη ≡ rT /rp ≥ 105 for Ra = 109
and rk ≥ 30, whose results can be found in Figs S15–S20.) In theregime diagram of Fig. 1 we classified the obtained flow patterns ofthermal convection as in our earlier works with the Cartesian geometry (Miyagoshi et al. 2015; Kameyama & Yamamoto 2018). Herethe regime of convective flow patterns is identified mainly based theratio of the average horizontal velocity Vtop at the top surface (r =rmax ) to the root-mean-square velocity Vrms of the entire layer shownin Fig. S1. The open circles and solid squares in the figure indicate
the sets of parameters Ra and rη where thermal convection takes
place in the ‘small viscosity contrast’ (SVC) and ‘stagnant-lid’ (ST)regimes, respectively. The solid triangle for rη = 105 and Ra ≥ 109 ,on the other hand, indicates the occurrence of so-called ‘sluggishlid’ (SL) regime (Kameyama & Ogawa 2000; Miyagoshi et al.
2015; Kameyama & Yamamoto 2018). The changes in convecting
flow patterns induced by the changes in the control parameters (Ra,rη and rk ) obtained here are in fairly good agreement with thoseobtained in our earlier Cartesian models (Miyagoshi et al. 2015;

Figure 1. The regime diagram of convection flow patterns in the plane of Ra
versus rη for the cases with purely temperature-dependent viscosity (rp =100 ). The open circles, solid triangle and solid squares indicate that thermal
convection occurs in the small-viscosity-contrast (SVC), the sluggish-lid(SL) and the stagnant-lid (ST) regimes, respectively. The lines show theplots against rη of the critical Rayleigh numbers Rac estimated by the linear
analysis (Kameyama 2021) for the values of rk indicated in the figure.

Kameyama & Yamamoto 2018). We also show by several lines theplots against rη of the critical Rayleigh numbers Rac estimated by
the linear analysis (Kameyama 2021) for the values of rk indicatedin the figure. These plots show that in general the values of Rac become larger for larger rη and rk , which can be intuitively understoodfrom an overall increase in viscosity and thermal conductivity withrη and/or rk . In particular, the values of Rac exceed 108 for the caseswith rk = 100 and rη  106 , which is consistent with the present
numerical results that small perturbations of thermal convection did
not grow with time under these conditions. This is the reason why
we showed the snapshots of the cases with Ra = 108 and rk = 100
only for rη ≤ 105 in Fig. S2.
In the following we will display the distributions of the potential
temperature Tpot as well as those of the ‘raw’ temperature T, in orderto clearly show the structures of thermal convection. Here, Tpot at
the dimensionless position (r, θ) is calculated by, r =rmax
α(r )dr .
(14)
Tpot (r, θ ) = T (r, θ ) exp −Dir =r

As will be seen later, the distributions of potential temperature Tpot
is more appropriate than that of temperature T for discussing thebuoyancy force that drives thermal convection, since the temperature change due to adiabatic (de)compression is subtracted. By
using (1) and (14), we obtain the dimensionless values of Tpot (r =rmax , θ) = Ttop = 0.1 andTpot (r = rmin , θ) = (1 + Ttop ) exp −Di

αave
αtop

 0.5

(15)

for the present numerical model.
We will next focus on the effects of the depth-dependence inthermal conductivity k on the thermal convection of highly compressible fluids with strongly temperature-dependent viscosity η. In
Figs 2 and 3 we present the overall convection structures for several
values of rk together with Ra = 109 and rη = 107 . In Fig. 2 we show

Downloaded from https://academic.oup.com/gji/article/231/2/1457/6638987 by Ehime University Library user on 05 August 2022

1
3
10
30
100

at r = rmax at r = rmin

Cartesian|dTcond /dz|

Highly compressible 2D cylindrical convection

(a)

(b)

(c)

(b)

(a)

1461

(c)

Figure 3. Comparison of overall structures of the thermal convection calculated for various values of rk together with Ra = 109 and rη = 107 . (a) The
distributions of potential temperature Tpot . The colour scale is indicated at
the bottom of the figure. The contour lines for Tpot are shown at the intervalof 0.005 only in the range of Tpot ≥ 0.47. Also shown for comparison by
the insets for 3 ≤ rk ≤ 30 are the distributions of Tpot obtained for the calculations in 2-D Cartesian box (Kameyama & Yamamoto 2018). (b) Plotsagainst dimensionless height from the bottom boundary of the horizontally
averaged potential temperature Tpot at the height of r. Plotted by the thinblue dashed lines are those obtained for the calculations in 2-D Cartesianbox. (c) Plots against r of the root-mean-squares of the magnitude of (red)velocity vector v 2 and the radial velocity vr 2 at the height r. Plotted22
in (c) are the values
 of  v and vr normalized by the maximum val-(a) the 2-D distributions of temperature T, (b) the plots against thedimensionless height from the bottom boundary (r − rmin ) of thehorizontally averaged temperature T at r and (c) the plots of thehorizontally averaged viscosity η against r. Also shown by thinblack lines and dotted lines in Fig. 2(b0 are the radial profiles oftemperatures given by several adiabats (i.e. the lines of constant
Tpot ) and steady 1-D heat conduction Tcond (r) in the radial direction, respectively. In Fig. 3, on the other hand, we show (a) the 2-D
distributions of potential temperature Tpot , (b) the plots of the horizontally averaged potential temperature Tpot at r and (c) the plotsof the root-mean-squares of the magnitude of velocity vector v 2
and those of radial velocity vr2 against r, calculated for the same
snapshots in Fig. 2. For the cases with 3 ≤ rk ≤ 30, in addition,
we show the results of our earlier Cartesian study (Kameyama &Yamamoto 2018) by the insets of Figs 2(a) and 3(a) and by the blue
dashed lines in Figs 2(b) and 3(b).

Figs 2 and 3 show that for the cases with Ra = 109 and rη =10 the thermal convection falls in the stagnant-lid (ST) regime
regardless of the values of rk . Indeed, in the shallow portion of
the fluid layer, the overall thermal structures are characterized by a
rapid increase in temperature T with depth in the thick (∼0.3d to0.4d) cold thermal boundary layers (TBLs) along the top surface.
In the thick cold TBLs, the fluid motion is negligibly small bothin the radial and horizontal directions (see Fig. 3c). Taken together
with their much higher viscosity than that at depth (see Fig. 2c), thecold TBLs act as stiff and stagnant lids along the top cold boundary.
In addition, the regions of ‘isentropic cores’ occur beneath the
thick stagnant lids where the temperature T gradually increases
with depth along with some adiabat (see Fig. 2b) and, in otherwords, the potential temperature Tpot is almost uniform (see Fig. 3b).
On the other hand, the comparison of the present results with ourearlier ones in the Cartesian domain (Kameyama & Yamamotoues vmax = max
in the figure.

7

v2

for given conditions, whose values are indicatedDownloaded from https://academic.oup.com/gji/article/231/2/1457/6638987 by Ehime University Library user on 05 August 2022

Figure 2. Comparison of overall structures of the thermal convection calculated with various values of rk ≡ kbot /ktop together with Ra = 109 and rη ≡
ηtop /ηbot = 107 . (a) The distributions of temperature T. The contour interval
is 0.05. The colour scale is indicated at the bottom of the figure. Also shownfor comparison by the insets for 3 ≤ rk ≤ 30 are the distributions of T obtained for the calculations in 2-D Cartesian box (Kameyama & Yamamoto2018). (b) Plots against dimensionless height from the bottom boundary of
the horizontally averaged temperature T (thick red lines) at the height of r.
Plotted by the thin blue dashed lines are those obtained for the calculationsin 2-D Cartesian box. Also shown for comparison by thin black lines are
the plots of several adiabats against r, while by dotted black lines are theradial profiles of temperature Tcond given by steady 1-D heat conduction in
the radial direction. (c) Plots against r of the horizontally averaged viscosity
η , normalized by the value of viscosity ηbot at the bottom boundary (r =rmin and T = Tbot ).

1462

M. Kameyama
comparison with those of steady-state conductive profiles dT/dr c .
For example, for the cases with rk ≤ 10 and cylindrical geometry, there occur thin regions with steep radial temperature gradientdT/dr just above the bottom boundary (i.e. TBLs), while suchregions are not observed for rk ≥ 30. The disappearance of TBLs
clearly indicates the weakening of the overall thermal convection bya large rk . For rk ≥ 30, in particular, the values of dT/dr are almost
equal to those of conductive profiles dT/dr c in the regions near
the bottom boundaries (r − rmin ≤ 0.25d for rk = 30 and r − rmin
≤ 0.4d for rk = 100). In addition, the values of the vertical gradient in the conductive temperature profiles dT/dr c become smaller
near the bottom surface for rk ≥ 30. This is because a large rk (ora large thermal conductivity k at the bottom boundary) enhances
a conductive heat transport at the bottom boundary even when thevertical gradient in temperature dT/dr c is gentle there. Taken together, for rk ≥ 30 where the profiles of dT/dr are close to thoseof dT/dr c near the bottom boundaries, the values of dT/dr are
smaller than the adiabatic gradient in temperature dT/dr s , implying that the thermal stratification is stable (Kameyama & Kinoshita
2013; Kameyama 2016, 2021). Indeed, for rk ≥ 30, the verticalgradient in horizontally averaged temperature dT/dr (thick redlines) is significantly gentler than that of the adiabatic one dT/dr s(dashed blue line) near the bottom boundary, indicating an occurrence of ‘deep stratosphere’ of stable thermal stratification (see also
Figs 2 and 3).
We here note that the ‘deep stratosphere’ is not identical with the‘CMB-lid’ proposed by Stamenković et al. (2012). In the numericalmodel of Stamenković et al. (2012) based on the Boussinesq approximation, the ‘CMB-lid’ is formed by an extremely high viscosity[up to around O(1030 ) Pa s] in the deep mantle of super-Earths, assuming that the viscosity significantly increases with pressure (ordepth). In the ‘deep stratosphere’ obtained here, in contrast, the fluid
motion is strongly suppressed in spite of smaller viscosity η owing
to its strong dependence on temperature T (see also Fig. 2c). The
difference between the ‘deep stratosphere’ and ‘CMB-lid’ significantly highlights the ultimate importance of adiabatic temperaturechange.
In order to demonstrate the ultimate importance of adiabatic
temperature change in the formation of ‘deep stratosphere’, we also
carried out an auxiliary calculation where the second term in theright-hand side of eq. (8) was artificially set to be zero, under thecondition where the ‘deep stratosphere’ otherwise occurs (Ra = 109 ,rη = 108 and rk = 100). From this calculation we confirmed that‘deep stratosphere’ never occurs in the absence of adiabatic temperature change even though a large increase in thermal conductivity
k with depth is introduced. In addition, we carried out further auxiliary calculations with 10 ≤ rk ≤ 100 where the reference density
ρ is kept constant with depth as in the cases under the extended
Boussinesq approximation (EBA) rather than the truncated anelastic liquid approximation (TALA). These calculations showed thatthe overall convection structures (such as an occurrence of ‘deep
stratosphere’) are almost unchanged by the changes in the depthdependence in the reference density. These results therefore indicatethat the effects of strong compression on thermal convection largely
come from the adiabatic temperature change.
The comparison between the plots with different geometries inFig. 4 also show the reason why the occurrence of ‘deep stratosphere’ is suppressed for the cylindrical geometry compared withthat for the Cartesian geometry even when the values of control
parameters (Ra, rη and rk ) are identical. We can clearly note fromthe figure that for given rk (≤30) the values of | dT/dr | in the lowermost part of the cylindrical cases are greater than those of | dT/dz |Downloaded from https://academic.oup.com/gji/article/231/2/1457/6638987 by Ehime University Library user on 05 August 20222018) for 3 ≤ rk ≤ 30 shows that the overall thermal structuresare very similar for the two model geometries. However, from adetailed comparison, we note that at a given height the temperature
T (and accordingly Tpot ) is slightly lower in the present cylindricalcases than in the Cartesian ones, because of the effects of model
geometries. Since the surface area decreases with depth (∝r−1 ) forthe cylindrical geometry, the magnitude of temperature gradient in
the radial (or vertical) direction in the shallow part of fluid layer issmaller in cylindrical cases than in Cartesian cases and, at the sametime, the magnitude in the deep part is larger. This results in lower
temperature T in the entire layer for the cylindrical cases than thatfor the Cartesian ones.
What should be also noted from the comparison between thecases in Figs 2 and 3 is the difference in the structures of thermal
convection in the broad regions beneath the thick stagnant lids dueto the changes in rk . As can be seen from Figs 2(a) and (b), thetemperature T in the isentropic cores becomes higher for larger rk ,leading to an increase in Tpot in the regions (see also Figs 3a andb). For rk ≥ 30, in particular, the values of potential temperature
Tpot in the isentropic core is higher than that at the bottom surface(∼0.5). That is, for rk ≥ 30 there occurs a distinct layer above thebottom boundary, called ‘deep stratosphere’ in our earlier works(Kameyama & Yamamoto 2018; Kameyama 2021), where the potential temperature Tpot increases with height and, in other words,
the thermal state is said to be stably stratified. In this layer the isolines for Tpot are almost horizontal, reflecting the fact that the fluid
motion is insignificant particularly in the radial direction in the layer(see Fig. 3c), in spite of low viscosity of hot fluids there (see also
Fig. 2c). In addition, from the comparison between the cases with rk= 30 and rk = 100 we note that for larger rk the value of Tpot in theisentropic core becomes larger and, in turn, the layer of insignificant
fluid motion becomes thicker. This is also consistent with our earlier
findings that the occurrence of ‘deep stratosphere’ is enhanced forlarger rk (Kameyama & Yamamoto 2018; Kameyama 2021). On theother hand, the comparison of the present results with our earlier
ones in the Cartesian domain (Kameyama & Yamamoto 2018) for3 ≤ rk ≤ 30 shows that the values of Tpot in the isentropic cores are
smaller in the present cylindrical cases than in the Cartesian ones,
which implies that the occurrence of ‘deep stratosphere’ tends to be
suppressed in the cylindrical cases.
In order to see the reason for the occurrence of ‘deep stratosphere’
for large increase in thermal conductivity k with depth, we next showby thick red lines in Fig. 4 the plots against r of the magnitudes of thevertical gradient of the horizontally averaged temperature dT/dr(≡ d T /dr) in the lower part of the fluid layer (0 ≤ r − rmin ≤ 0.5d)
calculated for 3 ≤ rk ≤ 100 together with Ra = 109 and rη = 107 .
Also plotted by the dashed blue lines and the dotted black lines in
the figure are the profiles of the magnitude of the adiabatic gradientin temperature dT/dr s (≡ −Diα T ), and those of the verticalgradient of steady-state conductive profile dT/dr c (≡ dTcond /dr),respectively. For the cases with rk ≤ 30, in addition, we show theresults of our earlier Cartesian study (Kameyama & Yamamoto2018). From an analogy of meteorological studies, the comparison
between dT/dr and dT/dr s enables us to see whether the thermal
stratification is stable or not at given r: Thermal stratification isstable when dT/dr is gentler than dT/dr s . On the other hand, thecomparison between dT/dr and dT/dr c indicates the vigour ofthermal convection: The profiles of dT/dr deviate from those ofdT/dr c more significantly for more vigorous convection.
The reason why the ‘deep stratosphere’ occurs for large rk ≡kbot /ktop can be understood from the comparison between the profiles of dT/dr and those of adiabatic gradient dT/dr s , through the

Highly compressible 2D cylindrical convection

(a)

(b)

(c)

1463

(d)

Table 2. Summary of the values of rη ≡ ηtop /ηbot
where ‘deep stratospheres’ are formed for the values of the Rayleigh number Ra and rk ≡ kbot /ktop
indicated in the table. The symbol ‘NoDS’ indicates that a ‘deep stratosphere’ does not occur
under these cases.

rk = 100
rk = 30
rk ≤ 10

Ra = 108

Ra = 109

Ra = 1010

≥ 105

≥ 105

rη ≥ 105rη ≥ 107NoDSrη
rη ≥ 106NoDSrη
rη ≥ 106
NoDS

of the Cartesian cases. This is again due to the nature of cylindrical
domains where the surface area decreases with depth (∝r−1 ): In
order to satisfy the equilibrium in the radial heat transfer, the radial
heat flux per unit area and, in other words, the temperature gradientin the radial direction must increase with depth along with the decrease with the surface area. Indeed, for the cases with cylindricalgeometry, the magnitude of dT/dr in the deeper part of the fluidlayer is larger than that in the shallower part (see also Table 1). Thisin turn makes the values of | dT/dr | in the lowermost part greaterin the cylindrical cases than those of | dT/dz | in the Cartesian cases
and, hence, makes the thermal stratification more unstable at depths
in the cylindrical cases.
We present in Figs S2–S14 in the Supplementary data the snapshots of the overall thermal convection for all the cases calculated
in this study using the purely temperature-dependent viscosity (rp= 1). In Table 2 we summarize the conditions under which the‘deep stratosphere’ occurs. The fact that the occurrence of ‘deep
stratosphere’ is limited to the cases with very strong increase inthermal conductivity with depth (rk ≥ 30) is consistent with thestatic stability of the thermal structure Tcond given by the steadystate conduction demonstrated in Fig. 4. For rk = 30 the thresholdvalues of rη above which the ‘deep stratosphere’ occurs increase
with the Rayleigh number Ra, indicating that ‘deep stratospheres’can be formed more readily when the convection is weaker and,hence, the thermal state at depth becomes closer to the conductive
ones. For rk = 100, on the other hand, the threshold values of rη
are not only larger than those for rk = 30 but also almost independent of Ra. The comparison of the results with rk = 30 and thosewith rk = 100 demonstrates that the occurrence of ‘deep stratosphere’ is greatly enhanced when rk is as large as 100. In addition,
based on the regime diagram shown in Fig. 1, the entire sets of Ra
and rη indicated above yield the thermal convection falling in thestagnant-lid regime. Taken together, the stagnant-lid convection ofhighly compressible fluids with very strong depth-dependent thermal conductivity can be accompanied by a ‘deep stratosphere’ of
insignificant motion at its base. These results are in good agreementwith our earlier findings by the linear analysis (Kameyama 2016,2021).
In order to see how the occurrence of ‘deep stratosphere’ is affected by the depth-dependence in viscosity, we carried out several
calculations with rp > 100 (rp = 101 , 102 and 103 ) together withRa = 109 , rk ≥ 30 and 105 ≤ rη ≤ 107 as indicated by the symbols
in Fig. S15. In Figs S16–S20 we compare the overall convection
structures for the cases with different rp but with the same valuesof rk and rη ≡ rT /rp . From the comparison for the cases with rk =30 (Figs S16 and S17) we note that the occurrence of ‘deep stratosphere’ is slightly suppressed for large rp ≥ 100 . This is because
the overall viscosity in the fluid layer is reduced with increasing rp

Downloaded from https://academic.oup.com/gji/article/231/2/1457/6638987 by Ehime University Library user on 05 August 2022

Figure 4. Plots against dimensionless height r − rmin from the bottom boundary of the rate of temperature change in the radial direction in the lower partof the fluid layer (0 ≤ r − rmin ≤ 0.5), calculated for various values of rk together with Ra = 109 and rη = 107 . Plotted by thick red lines are the profilesof the absolute values of the radial gradient of the horizontally averaged temperature | dT/dr |. Also shown by dashed blue lines are the plots for adiabatic
gradient in temperature | dT/dr s | calculated for the radial profiles of horizontally averaged temperature T , while by dotted black lines are those for steady
1-D heat conduction | dT/dr c | ≡ |dTcond /dr|. For rk ≤ 30 we also show the similar plots obtained for the calculations in 2-D Cartesian geometry (Kameyama& Yamamoto 2018).

1464

M. Kameyama

(a)

(b)

By definition, Qtop is the average heat flux at the top surface (r = rmax ) normalized by the scale of heat flux throughthe fluid layer (≡ ktop T/d). On the other hand, Nu is definedbyN u ≡ Q top ×dTcond
dr

−1

,

(17)

(r =rmax )where the values of the radial gradient of the steady-state conductive profile Tcond at the top surface are shown in Table 1 for several
values of rk . By definition, Nu is a measure of the efficiency of convective heat transfer after the correction of the enhanced conductive
heat transfer for rk > 1 as well as the change in the surface areawith depth and, in other words, a measure of the vigour of overall
thermal convection. From these plots we note that the values of
Nu become smaller with increasing rη for given rk , which comesfrom the decrease in Qtop associated with the transition into thestagnant-lid regime (see also Fig. S1). Similarly, the values of Nu
are smaller for larger rk for given rη , which reflects the fact that thevigour of thermal convection is reduced with increasing rk becauseof the increase in the critical Rayleigh number Rac with rk (see also
Fig. 1).
The plots of Nu in Fig. 6 show that the changes in Nu with rη andrk are consistent with the occurrence of ‘deep stratosphere’ for thecases with large rη and rk . We can clearly note that the values of
Nu are very small for the cases of stagnant-lid regime with large rη ,reflecting weak convection. In particular, for the cases where bothrη and rk are large enough for the occurrence of ‘deep stratosphere’,
the values of Nu are very close to unity. That is, the stagnant-lid
convection with ‘deep stratosphere’ is an extremely weak convection occurring very close to the critical condition for very large rη
and rk . In fact, such flow structures are very similar to those ofincipient convection for very large rη and rk obtained by our earlier
linear analyses (Kameyama 2016, 2021), which are characterized by
two regions of insignificant fluid motion along the top and bottom
surfaces.
4 C O N C LU S I O N
In this paper we carried out a series of numerical experiments ofthermal convection of highly compressible fluids whose physical(c)

Figure 5. Plots against rη of the heat flux Qtop across the top surface for thecases with various values of rk for (a) Ra = 108 , (b) Ra = 109 and (c) Ra =1010 .

Downloaded from https://academic.oup.com/gji/article/231/2/1457/6638987 by Ehime University Library user on 05 August 2022

under the assumption of constant rη (see, for example, Fig. S16h),
resulting from the increase in rT = rη × rp in proportion to rp . The
thermal convection is enhanced for large rp , which in turn tends tosuppress the occurrence of ‘deep stratosphere’. The comparison forthe cases with rk = 100 (Figs S18–S20), on the other hand, showsthat the values of rp do not significantly affect the occurrence of‘deep stratosphere’, at least for the cases with Ra = 109 and rη ≥105 , even though the overall viscosity in the fluid layer is reducedwith increasing rp . This is consistent with the results with higher Ra= 1010 , where the occurrence of ‘deep stratosphere’ is greatly enhanced when rk is as large as 100 despite the lower overall viscosity(see Fig. S12).
In order to further see the changes in the vigour and flow patterns of thermal convection, we show in Fig. 5 the plots againstrη of the heat flux Qtop at the top surface (r = rmax ) calculated forvarious values of rk and Ra, while in Fig. 6 we show those of theNusselt number Nu. Here Qtop is defined in a dimensionless formby
 θ=θmax
∂T
1
dθ.
(16)
Q top =
θmax − θmin θ=θmin∂r (r =rmax ,θ)Highly compressible 2D cylindrical convection

(a)

(c)

Figure 6. Plots against rη of the Nusselt number Nu for the cases withvarious values of rk for (a) Ra = 108 , (b) Ra = 109 and (c) Ra = 1010 .

properties strongly vary in space in a 2-D cylindrical annulus, byextending our earlier work in the Cartesian geometry (Kameyama &Yamamoto 2018): the variations in thermodynamic properties (thermal expansivity and reference density) with depth are taken to be
relevant for the super-Earths with 10 times the Earth’s mass, whilethe viscosity and thermal conductivity are assumed to exponentiallydepend on temperature and/or depth. Our experiments demonstrated
that, for the cases with strongly temperature-dependent viscosity ηand depth-dependent thermal conductivity k, there occurs a distinct
layer at the base of the mantle in addition to thick stagnant lids attheir top cold surfaces, as had been pointed out in our earlier Cartesian study (Kameyama & Yamamoto 2018). The distinct basal layer,termed a ‘deep stratosphere’ here, is characterized by a staticallystable thermal stratification where a fluid motion is insignificant particularly in the radial direction, in spite of low viscosity of hot fluids
there. We also found that the occurrence of ‘deep stratosphere’ tends
to be suppressed for the cases with cylindrical geometry comparedwith those with Cartesian geometry. This is because the reductionof surface area with depth in the cylindrical geometries steepen thegradient in temperature in the radial direction at depth, which in turn
reduces the static stability of thermal stratification there. Our finding may further imply that both the effects of adiabatic compression
and those of cylindrical (or spherical) geometry of mantle are of
crucial importance in understanding the mantle dynamics of massive super-Earths in the presence of spatial variations in physicalproperties.
From an analogy of the effects of 2-D cylindrical geometry
demonstrated here, we can note that the effects of 3-D sphericalgeometry are most likely to further suppress the occurrence of‘deep stratosphere’ in the mantles of real super-Earths, since thedecrease in the surface area with depth is more significant in 3-D
cases (∝r−2 ) than in 2-D cases (∝r−1 ). However, our earlier linear
analysis (Kameyama 2021) showed that the stagnant-lid convectionof highly compressible fluids with strongly temperature-dependent
viscosity is accompanied by another motionless regions above thebottom boundaries even in the 3-D spherical geometry in the presence of a factor of about 100 increase in the thermal conductivity
with depth. We can therefore speculate that the ‘deep stratospheres’of insignificant fluid motion can be formed above the bottom boundaries if the thermal conductivity k of their mantle materials strongly
increases with depth (or pressure), although the increase in k may
not necessarily be so large as assumed here in silicate mantles ofreal super-Earths (e.g. Dalton et al. 2013; Hsieh et al. 2017, 2018;
Deschamps & Hsieh 2019).
We do nevertheless believe that ‘deep stratospheres’ can occur inmassive super-Earths in the universe, considering that its occurrence
can be also greatly affected by the thermal state in the fluid layer.
This is because a high temperature in the interior of the layer increases the adiabatic gradient in temperature, which in turn stabilize
the thermal stratification (Kameyama & Kinoshita 2013; Tachinami
et al. 2014; Kameyama et al. 2015). Recall that in this study thethermal state in the convecting fluid layer is given by a statisticallysteady state without internal heat sources, which in essence results
in a colder state than expected in nature. That is, taking into accountthe effects of internal heating by radiogenic elements and the courseof secular cooling of planetary interiors from their hot (and possibly
molten) initial states, a hotter thermal state is expected in the mantleof super-Earth than that assumed in this study. This will be particularly true for the thick mantles of massive super-Earths, whose
thermal evolution (or cooling) is expected to be extremely slow
(Miyagoshi et al. 2017). In addition, for super-Earths which haveclose-in orbits around their parent stars, the temperature at their top

Downloaded from https://academic.oup.com/gji/article/231/2/1457/6638987 by Ehime University Library user on 05 August 2022

(b)

1465

1466

M. Kameyama

AC K N OW L E D G M E N T S
We thank Masaki Ogawa, Takehiro Miyagoshi, Takatoshi Yanagisawa, Sebastian Ritterbex and Haruhiko Dekura for fruitful discussion.
FUNDING
This work was supported by JSPS KAKENHI Grant Number
18K03724, a joint research program at Geodynamics Research Center, Ehime University, the Earth Simulator of the Japan Agency for
Marine-Earth Science and Technology (JAMSTEC), ‘Exploratory
Challenge on Post-K computer’ (Elucidation of the Birth of Exoplanets [Second Earth] and the Environmental Variations of Planetsin 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).
D ATA AVA I L A B I L I T Y
The data underlying this paper will be shared on reasonable request
to the corresponding author.
C o n fl i c t o f I n t e r e s t
The author declares no conflict of interest.
REFERENCES
Dalton, D.A., Hsieh, W.-P., Hohensee, G.T., Cahill, D.G. & Goncharov, A.F.,
2013. Effect of mass disorder on the lattice thermal conductivity of MgO
periclase under pressure, Sci. Rep., 3(2400), doi:10.1038/srep02400.
Dekura, H. & Tsuchiya, T., 2017. Ab initio lattice thermal conductivityof MgO from a complete solution of the linearized Boltzmann transport
equation, Phys. Rev. B, 95( 18), doi:10.1103/PhysRevB.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), 12 919–12 926.
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), doi:10.1103/PhysRevLett.110.025904.
Deschamps, F. & Hsieh, W.-P., 2019. Lowermost mantle thermal conductivity constrained from experimental data and tomographic models,Geophys. J. Int., 219(Suppl. 1), S115–S136.
Foley, B.J., 2019. Habitability of Earth-like stagnant lid planets: climate
evolution and recovery from snowball states, Astrophys. J., 875(1), 72,
doi:10.3847/1538-4357/ab0f31.
Gaidos, E., Conrad, C.P., Manga, M. & Hernlund, J., 2010. Thermodynamic
limits on magnetodynamos in rocky exoplanets, Astrophys. J., 718(2),596–609.
Guerrero, J.M., Lowman, J.P., Deschamps, F. & Tackley, P.J., 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., 123(7), 1863–1880.
Haigis, V., Salanne, M. & Jahn, S., 2012. Thermal conductivity of MgO,
MgSiO3 perovskite and post-perovskite in the Earth’s deep mantle, Earthplanet. Sci. Lett., 355-356, 102–108.
Hsieh, W.-P., Deschamps, F., Okuchi, T. & Lin, J.-F., 2017. Reduced lattice
thermal conductivity of Fe-bearing bridgmanite in Earth’s deep mantle,
J. geophys. Res., 122(7), 4900–4917.
Hsieh, W.-P., Deschamps, F., Okuchi, T. & Lin, J.-F., 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.
Jarvis, G.T. & McKenzie, D.P., 1980. Convection in a compressible fluidwith infinite Prandtl number, J. Fluid Mech., 96(3), 515–583.
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.
Kameyama, M., 2021. Linear analysis on the onset of thermal convectionof highly compressible fluids with variable viscosity and thermal conductivity in spherical geometry: implications for the mantle convectionof super-Earths, Earth, Planets Space, 73, 67, doi:10.1186/s40623-02101499-w.
Kameyama, M. & Harada, A., 2017. Supercontinent cycle and thermochemical structure in the mantle: inference from two-dimensional
numerical simulations of mantle convection, Geosciences, 7(4),doi:10.3390/geosciences7040126.
Kameyama, M. & Kinoshita, Y., 2013. On the stability of thermal stratification of highly compressible fluids with depth-dependent physical properties: implications for the mantle convection of super-Earths, Geophys. J.
Int., 195(3), 1443–1454.
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.
Kameyama, M. & Yamamoto, M., 2018. Numerical experiments on thermal convection of highly compressible fluids with variable viscosity andthermal conductivity: implications for mantle convection of super-Earths,
Phys. Earth planet. Inter., 274, 23–36.
Kameyama, M., Miyagoshi, T. & Ogawa, M., 2015. Linear analysis on theonset of thermal convection of highly compressible fluids: implications
for the mantle convection of super-Earths, Geophys. J. Int., 200(2), 1066–
1077.
Karato, S.-i., 2011. Rheological structure of the mantle of a super-Earth:some insights from mineral physics, Icarus, 212(1), 14–23.
King, S.D., Lee, C., van Keken, P.E., Leng, W., Zhong, S., Tan, E., Tosi,
N. & Kameyama, M., 2010. A community benchmark for 2-D Cartesiancompressible convection in the Earth’s mantle, Geophys. J. Int., 180(1),73–87.
Kite, E.S., Manga, M. & Gaidos, E., 2009. Geodynamics and rate of volcanism on massive Earth-like planets, Astrophys. J., 700(2), 1732–1749.
Korenaga, J., 2010. On the likelihood of plate tectonics on super-Earths:
does size matter?, Astrophys. J., 725(1), L43–L46.

Downloaded from https://academic.oup.com/gji/article/231/2/1457/6638987 by Ehime University Library user on 05 August 2022

surface is expected to be much higher than that assumed here (i.e.
Ttop /T > 0.1), which is most likely to further help maintain hot
thermal states in their interiors. Taken these effects together, our results are most likely to underestimate, rather than overestimate, thecondition for the occurrence of ‘deep stratosphere’ in the mantles
of super-Earths in nature, although an exact knowledge is neededon the physical properties (thermal expansivity and conductivity)
as well as those on thermal state in order to precisely estimate forparticular planets.
What is also worth being noted in our findings is that the fluidin the ‘deep stratosphere’ is stagnant owing to the stable thermal
stratification in spite of its low viscosity. This may further imply
that such regions of insignificant fluid motion can be formed more
readily at the base of the mantles of massive super-Earths in thepresence of strong increase in viscosity with depth (or pressure),even though the viscosity is not likely to increase hugely enoughto make up ‘CMB-lids’ (Stamenković et al. 2012) according to therecent studies on mineral physics (e.g. Karato 2011; Ritterbex et al.
2018).
The occurrence of ‘deep stratosphere’ at the base of the mantles
of massive super-Earths may strongly influence the mantle dynamics and surface environments of the particular terrestrial planets.
For example, this distinct layer 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, this layer 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.

Highly compressible 2D cylindrical convection

Tosi, N. et al., 2017. The habitability of a stagnant-lid Earth, Astoron.
Astrophys., 605, A71, doi:10.1051/0004-6361/201730728.
Tsuchida, M. & Kameyama, M., 2020. 2-D numerical simulations on formation and descent of stagnant slabs: important roles of trench retreat andits temporal change, Front. Earth Sci., 8, doi:10.3389/feart.2020.00117.
Tsuchiya, T. & Tsuchiya, J., 2011. Prediction of a hexagonal SiO2 phaseaffecting stabilities of MgSiO3 and CaSiO3 at multimegabar pressures,
Proc. Nat. Acad. Sci. U.S.A., 108(4), 1252–1255.
Valencia, D., O’Connell, R.J. & Sasselov, D.D., 2006. Internal structure of
massive terrestrial planets, Icarus, 181(2), 545–554.
Valencia, D., O’Connell, R.J. & Sasselov, D.D., 2007a. Inevitability of plate
tectonics on super-Earths, Astrophys. J., 670(1), L45–L48.
Valencia, D., Sasselov, D.D. & O’Connell, R.J., 2007b. Radius and structure
models of the first super-Earth planet, Astrophys. J., 656(1), 545–551.
van den Berg, A.P., Yuen, D.A., Umemoto, K., Jacobs, M. H.G. & Wentzcovitch, R.M., 2018. Mass-dependent dynamics of terrestrial exoplanets
using ab initio mineral properties, Icarus, 317, 412–426.
Čı́žková, H., van den Berg, A. & Jacobs, M., 2017. Impact of compressibility
on heat transport characteristics of large terrestrial planets, Phys. Earthplanet. Inter., 268, 65–77.
Wagner, F.W., 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.
Yanagisawa, T., Kameyama, M. & Ogawa, M., 2016. Numerical studies
on convective stability and flow pattern in three-dimensional sphericalmantle of terrestrial planets, Geophys. J. Int., 206(3), 1526–1538.

S U P P O RT I N G I N F O R M AT I O N
Supplementary data are available at GJ I online.
Figure S1. Plots against rη ≡ ηtop /ηbot of the time-averaged valuesof (a, c and e) the average lateral velocity Vtop at the top surface and(b, d and f) the ratio of Vtop to the root-mean-square velocity Vrms
in the entire layer, obtained for the cases with various values of Ra
and rk indicated in the figure.
Figure S2. Comparison of overall structures of the thermal convection calculated with various values of rη ≡ ηtop /ηbot together with Ra= 108 and rk ≡ kbot /ktop = 100. (a) The distribution of temperature
T. The contour interval is 0.05. The colour scale is indicated at thebottom of the figure. (b) The distribution of potential temperature
Tpot . The colour scale is indicated at the bottom of the figure. The
contour lines are shown at the interval of 0.005 in the range ofTpot ≥ 0.47. (c) Plots against dimensionless height r − rmin fromthe bottom boundary of the horizontally averaged temperature T(thick red lines) at r. Also shown for comparison by thin black lines
are the plots of several adiabats against r, while by dotted blacklines are the radial profiles of temperature Tcond given by steady 1-D
heat conduction in the radial direction. (d) Plots against r − rmin of
the horizontally averaged potential temperature Tpot at r. (e) Plotsagainst r − rmin of the root-mean-squares of the magnitude of (red)v 2 and the radial velocity
vr 2 at r. Plotted
velocity vector22
are the values of  v andvr normalized by the maximum
values vmax = max
v 2 for given conditions, whose values areindicated in the figure.
Figure S3. Same as Fig. S2 but for Ra = 108 and rk = 30.
Figure S4. Same as Fig. S2 but for Ra = 108 and rk = 10.
Figure S5. Same as Fig. S2 but for Ra = 108 and rk = 3.
Figure S6. Same as Fig. S2 but for Ra = 108 and rk = 1.
Figure S7. Same as Fig. S2 but for Ra = 109 and rk = 100.
Figure S8. Same as Fig. S2 but for Ra = 109 and rk = 30.
Figure S9. Same as Fig. S2 but for Ra = 109 and rk = 10.
Figure S10. Same as Fig. S2 but for Ra = 109 and rk = 3.
Figure S11. Same as Fig. S2 but for Ra = 109 and rk = 1.

Downloaded from https://academic.oup.com/gji/article/231/2/1457/6638987 by Ehime University Library user on 05 August 2022

Leng, W. & Zhong, S., 2008. Viscous heating, adiabatic heating and energetic
consistency in compressible mantle convection, Geophys. J. Int., 173(2),693–702.
Liu, X. & Zhong, S., 2013. Analyses of marginal stability, heat transfer andboundary layer properties for thermal convection in a compressible fluidwith infinite Prandtl number, Geophys. J. Int., 194(1), 125–144.
Lobanov, S.S., Holtgrewe, N., Lind, J.-F. & Goncharov, A.F., 2017. Radiative
conductivity and abundance of post-perovskite in the lowermost mantle,
Earth planet. Sci. Lett., 479, 43–49.
Miyagoshi, T., Tachinami, C., Kameyama, M. & Ogawa, M., 2014. On
the vigor of mantle convection in super-Earths, Astrophys. J., 780(L8),doi:10.1088/2041-8205/780/1/L8.
Miyagoshi, T., Kameyama, M. & Ogawa, M., 2015. Thermal convection andthe convective regime diagram in super-Earths, J. geophys. Res., 120(7),1267–1278.
Miyagoshi, T., Kameyama, M. & Ogawa, M., 2017. Extremely long transition phase of thermal convection in the mantle of massive super-Earths,
Earth, Planets Space, 69(1), doi:10.1186/s40623-017-0630-6.
Miyagoshi, T., Kameyama, M. & Ogawa, M., 2018. Effects of adiabatic compression on thermal convection in super-Earths of various sizes, Earth,Planets Space, 70, doi:10.1186/s40623-018-0975-5.
Miyauchi, A. & Kameyama, M., 2013. Influences of the depth-dependenceof thermal conductivity and expansivity on thermal convection withtemperature-dependent viscosity, Phys. Earth planet. Inter., 223, 86–95.
Noack, L., Rivoldini, A. & Van Hoolst, T., 2017. Volcanism and outgassingof stagnant-lid planets: implications for the habitable zone, Phys. Earthplanet. Inter., 269, 40–57.
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.
Ohta, K., Yagi, T., Hirose, K. & Ohishi, Y., 2017. Thermal conductivity offerropericlase in the Earth’s lower mantle, Earth planet. Sci. Lett., 465,
29–37.
Papuc, A.M. & Davies, G.F., 2008. The internal activity and thermal evolution of Earth-like planets, Icarus, 195(1), 447–458.
Patankar, S.V., 1980. Numerical Heat Transfer and Fluid Flow, Computational Methods in Mechanics and Thermal Sciences, Hemisphere Publishing Corporation (CRC Press, Taylor & Francis Group).
Ricard, Y., Alboussière, T., Labrosse, S., Curbelo, J. & Dubuffet, F., 2022.
Fully compressible convection for planetary mantles, Geophys. J. Int.,
230(2), 932–956.
Ritterbex, S., Harada, T. & Tsuchiya, T., 2018. Vacancies in MgO at ultrahighpressure: about mantle rheology of super-Earths, Icarus, 305, 350–357.
Rivera, E.J. et al., 2005. A ∼7.5 M⊕ planet orbiting the nearby star, GJ 876,
Astrophys. J., 634(1), 625–640.
Schubert, G., Turcotte, D.L. & Olson, P., 2001. Mantle Convection in the
Earth and Planets, pp. 940, Cambridge Univ. Press.
Seager, S., Kuchner, M., Hier-Majumder, C.A. & Militzer, B., 2007. Massradius relationships for solid exoplanets, Astrophys. J., 669(2), 1279–
1297.
Shahnas, M.H. & Pysklywec, R.N., 2021. Focused penetrative plumes: apossible consequence of the dissociation transition of post-perovskite at∼0.9 TPa in massive rocky super-Earths, Geochem. Geophys. Geosyst.,
22(8), e2021GC009910, doi:10.1029/2021GC009910.
Sotin, C., Grasset, O. & Mocquet, A., 2007. Mass-radius curve for extrasolar
Earth-like planets and ocean planets, Icarus, 191(1), 337–351.
Stamenković, V., Breuer, D. & Spohn, T., 2011. Thermal and transport
properties of mantle rock at high pressure: applications to super-Earths,
Icarus, 216(2), 572–596.
Stamenković, V., Noack, L., Breuer, D. & Spohn, T., 2012. The influenceof pressure-dependent viscosity on the thermal evolution of super-Earths,
Astrophys. J., 748(41), doi:10.1088/0004-637X/748/1/41.
Tachinami, C., Senshu, H. & Ida, S., 2011. Thermal evolution and lifetimeof intrinsic magnetic fields of super-Earths in habitable zones, Astrophys.
J., 726(2), doi:10.1088/0004-637X/726/2/70 .
Tachinami, C., Ogawa, M. & Kameyama, M., 2014. Thermal convection ofcompressible fluid in the mantle of super-Earths, Icarus, 231, 377–384.

1467

1468

M. Kameyama

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by theauthors. Any queries (other than missing material) should be directed to the corresponding author for the paper.

A P P E N D I X A : B A S I C E Q UAT I O N S I N
2-D CYLINDRICAL POLARC O O R D I N AT E S
The basic equations in 2-D cylindrical polar coordinates (r, θ) canbe explicitly written as
1 ∂1 ∂(r ρvr ) +(ρvθ ) ,
r ∂r
r ∂θ


2
1 ∂
∂
[2ηr (ε̇rr − ε̇θθ )]p + η(∇ · v) +
0=−
∂r
3
r ∂r
∂
1 ∂
(2ηε̇θθ ) +(2ηε̇r θ ) + Raρ αT,
+
∂r
r ∂θ


2
1 ∂
p + η(∇ · v)
0=−
r ∂θ
3
 1 ∂
1 ∂ 
(2ηε̇θθ ) ,2ηr 2 ε̇θr +
+ 2
r ∂r
r ∂θ


∂T
∂T
1 ∂
ρ
=−
r ρT vr − k
∂t
r ∂r
∂r
1 ∂T
1 ∂
ρT vθ − k
−
r ∂θ
r ∂θ
Di
− Diρ α T vr +2η (ε̇ II )2 ,
Ra
0=

(A1)

(A2)

(A3)

(A4)

where the components of strain rate tensor ε̇ are defined by∂vr,
ε̇rr =
∂r
1 ∂vθ
∂vr
vr
1 ∂vθ
1 ∂
ε̇θθ =+=+(r vr ) −
,
r ∂θ
r
r ∂θ
r ∂r
∂r

ε̇r θ = ε̇θr =∂  vθ  1 ∂vr
1r+,
2 ∂r rr ∂θ

and its second invariant ε̇ II can be written as
 

1 2
1
2
ε̇rr + ε̇θθε̇ II =+ 2ε̇r2θ − (ε̇rr + ε̇θθ )2 .
2
3

(A6)

(A8)

APPENDIX B: NUMERICALT E C H N I Q U E F O R S O LV I N G F L O W
FIELDS IN TERMS OF A
‘GENERALIZED’ STREAM FUNCTION ψ
In order to calculate the flow field under the truncated anelastic
liquid approximation (TALA), one should simultaneously solve theequations for continuity (eq. A1) and momentum in r-direction (eq.
A2) and θ-direction (eq. A3). We solve the simultaneous equations by the use of a ‘generalized’ stream function ψ. In this section we obtain the basic equation for flow fields in terms of ψ.
In 2-D cylindrical polar coordinates, the ‘generalized’ stream
function ψ is defined by1 ∂ψ∂ψ, ρvθ = −
.
(B1)
r ∂θ
∂r
Note that ψ is chosen to make the velocity field v automatically satisfy the equation of continuity (eq. A1). On the other hand, by taking
the rotation of eqs (A2) and (A3) we can eliminate the pressure p
and isotropic part (2/3)η(∇ · v) of flow stress as
ρvr =

Raρ α

1 ∂T
1 ∂2
[2ηr (ε̇θθ − ε̇rr )]= 2r ∂θr ∂r ∂θ

1 ∂2
1 ∂ 1 ∂ 2ηr 2 ε̇θr − 2 2 (2ηε̇r θ ) (B2)+r ∂r r ∂r
r ∂θ
 

2
∂
∂vr
∂vθ
1 ∂
2η
(r vr ) +− 2r
= 2r ∂r ∂θ
∂r
∂θ
∂r
 

v 
∂v
∂
1 ∂ 1 ∂
θ
r
3
η r
+r
+
r ∂r r ∂r
∂r r
∂θ





2
1 ∂vr
1 ∂
∂ vθ
+
.
(B3)
− 2 2 η r
r ∂θ
∂r r
r ∂θ

In deriving this equation we implicitly used that the thermodynamic
properties (reference density ρ and thermal expansivity α) depend
only on r. By further substituting (eq. B1) into this equation, we
obtain a generalized biharmonic equation containing solely ψ as, 
1 ∂2∂ 1 ∂ψ
1 ∂ 2ψ
1 ∂T
= 2
2η
−
Raρ α
r ∂θ
r ∂r ∂θ∂r ρ ∂θρ ∂r ∂θ
∂ 1 1 ∂ψ− 2r
∂r ρ r ∂θ∂ 1 1 ∂ψ1 ∂ 1 ∂−
ηr 3
r ∂r r ∂r∂r r ρ ∂r
1 ∂ 2ψ
1 ∂ 1 ∂
η+r ∂r r ∂rρ ∂θ 2∂ 1 1 ∂ψ1 ∂2
η+r ∂θ 2∂r r ρ ∂r
−

(A5)

(A7)

1 1 ∂2
ρ r 4 ∂θ 2

η

∂ 2ψ
∂θ 2

.

(B4)

By discretizing the eq. (B4) we get simultaneous linear equations forψ, which can be solved by the use of a direct method such as theGaussian elimination and the LU factorization. In this study, we

Downloaded from https://academic.oup.com/gji/article/231/2/1457/6638987 by Ehime University Library user on 05 August 2022

Figure S12. Same as Fig. S2 but for Ra = 1010 and rk = 100.
Figure S13. Same as Fig. S2 but for Ra = 1010 and rk = 30.
Figure S14. Same as Fig. S2 but for Ra = 1010 and rk = 10.
Figure S15. Same as Fig. 1 not only for the case with (a) rp = 100
but also for the cases with (b) rp = 101 , (c) 102 and (d) 103 .
Figure S16. Comparison of overall structures of the thermal convection calculated with various values of rp together with Ra =109 , rk = 30 and rη = rT /rp = 106 . (a–e) Same as Fig. S2. (f) Plotsagainst dimensionless height r − rmin from the bottom boundary of
the horizontally averaged temperature T obtained for the valuesof rp indicated in the figure. Also shown for comparison by dotted black lines are the radial profiles of temperature Tcond given bysteady 1-D heat conduction in the radial direction. (g) Plots againstr − rmin of the horizontally averaged potential temperature Tpot
obtained for the values of rp indicated in the figure. (h) Plots againstr − rmin of the horizontally averaged viscosity η obtained for thevalues of rp indicated in the figure. Plotted are the values of η
normalized by the values of viscosity ηbot at the bottom boundary.
Figure S17. Same as Fig. S16 but for Ra = 109 , rk = 30 and rη =rT /rp = 107 .
Figure S18. Same as Fig. S16 but for Ra = 109 , rk = 100 and rη =rT /rp = 105 .
Figure S19. Same as Fig. S16 but for Ra = 109 , rk = 100 and rη =rT /rp = 106 .
Figure S20. Same as Fig. 20 but for Ra = 109 , rk = 100 and rη =rT /rp = 107 .

Highly compressible 2D cylindrical convectionused the LAPACK subroutines (dgbtrf and dgbtrs) designed forlinear equations whose coefficient matrix is a non-symmetric band
matrix. A major benefit of direct methods is the robustness against
the spatial variations in viscosity η: solutions can be obtained at thesame computational costs regardless of the variations in η, as long1469

as the coefficient matrix is safely factorized. Our experience showed
that, however, iterative methods like the conjugate gradient method
are not suitable to solve the linear equations coming from eq. (B4)even when η is constant, since their coefficient matrix cannot be
positive definite.

Downloaded from https://academic.oup.com/gji/article/231/2/1457/6638987 by Ehime University Library user on 05 August 2022