1

2-D numerical experiments of thermal convection of highly vis-

2

cous fluids under strong adiabatic compression: Implications on

3mantle convection of super-Earths with various sizes

4

Masanori Kameyama, Geodynamics Research Center (GRC), Ehime University, Matsuyama,
Ehime 790-8577, Japan, kameyama.masanori.my@ehime-u.ac.jp
Also at Japan Agency for Marine-Earth Science and Technology (JAMSTEC), Yokohama,
Kanagawa 236-0001, Japan.

5

Preprint submitted to Earth, Planets and Space on August 1, 2024
1

6

Abstract

7

We conduct a series of numerical experiments of thermal convection of compressible fluids with8

temperature-dependent viscosity, in order to study how the adiabatic compression and model

9

geometries affect the mantle convection on super-Earths. A two-dimensional basally-heated convection

10

is considered under the truncated anelastic liquid approximation (TALA), either in a rectangular box or11

in a cylindrical annulus. We varied the magnitude of adiabatic heating and the Rayleigh number as well12

as the depth profile of thermodynamic properties (thermal expansivity and reference density) in13

accordance with the planetary sizes. From our calculations by varying the planetary sizes up to 10

14times the Earth’s mass, we confirmed that the adiabatic compression affects the thermal convection

15

more strongly for larger planets. The activity of hot plumes originating from the core-mantle boundary

16

is significantly suppressed in the terrestrial planets whose mass is larger than the Earth’s by a factor of

17

about 3 regardless of the model geometries. We also developed scaling relationships between the vigor

18of thermal convection and the planetary mass by appropriately incorporating the effect of adiabatic

19

compression into those of Boussinesq (or incompressible) cases. Our scaling relationships suggest that20

the stress level in the top cold thermal boundary layers is almost independent of the planetary mass,

21

which may further imply that the emergence of plate tectonics is not likely to be enhanced for massive

22

terrestrial planets larger than our Earth.

23

Keywords

24

mantle convection, super-Earths, adiabatic compression

25

Introduction

26

Since the first discovery of extra-solar “super-Earths” (Rivera et al. 2005), a growing effort has been

27

devoted to the studies in the internal (radial) structure (e.g., Seager et al. 2007; Sotin et al. 2007; Wagner

28

et al. 2011), the thermal evolution (e.g., Papuc and Davies 2008; Kite et al. 2009; Korenaga 2010; Gaidos

29

et al. 2010; Tachinami et al. 2011; Čı́žková et al. 2017) and the habitability (e.g., Noack et al. 2017; Tosi30

et al. 2017; Foley 2019) of such massive terrestrial bodies. During the last decades, the studies on mantle

31

dynamics have been playing an important role in addressing these issues (e.g., Kameyama and Yamamoto2

32

2018; van den Berg et al. 2019; Shahnas and Pysklywec 2021; Kameyama 2022; Ricard and Alboussière332023).

34

It has been well acknowledged that the effect of adiabatic compression affects the nature of thermal

35

convection (Jarvis and McKenzie 1980; King et al. 2010; Liu and Zhong 2013). This effect is expected to36

be of greater importance in the convection in the thicker mantles of larger super-Earths, considering that37

the pressure can reach on the order of terapascals in their interiors (Tsuchiya and Tsuchiya 2011). Indeed,38

our earlier studies (Miyagoshi et al. 2014, 2015) demonstrated that adiabatic compression substantially39

reduces the vigor of thermal convection in a massive super-Earth whose mass is ten times the Earth’s,40

in contrast to the predictions of many of earlier studies that are based on the Boussinesq approximation

41

where adiabatic compression is neglected (e.g., Valencia et al. 2007; Valencia and O’Connell 2009; van

42

Heck and Tackley 2011; Foley et al. 2012). In particular, the heart of the suppressing effect of adiabatic

43

compression lies in the adiabatic temperature change which tends to decelerate the vertical motion of44

fluid parcels (Tachinami et al. 2014; Kameyama and Kinoshita 2013; Kameyama 2016, 2021).

45

In order to discuss how strongly adiabatic compression affects the thermal convection depending on

46

planetary mass, Miyagoshi et al. (2018) carried out numerical experiments of thermal convection of47

a compressible fluid in the two-dimensional rectangular box by systematically varying the physical48

parameters with increasing the mass of modeled planets. They qualitatively confirmed that the adiabatic

49

compression affects more strongly the thermal convection in the mantles of planets with larger mass.

50

Based on somewhat simple analysis and visual comparison of their numerical results, the threshold mass51

of modeled terrestrial planets was estimated above which the effects of adiabatic compression become52

dominant on their mantle convection.

53

This paper is an extension or revision of the work by Miyagoshi et al. (2018), aiming at detailed and54

improved discussion on the effects of planetary sizes on the thermal convection in the mantles. From

55

an analogy of the mantle convection of terrestrial planets in our solar system whose dynamic nature is56

well characterized by the lithosphere and hot ascending plumes (e.g., Schubert et al. 2001; Davies 2011),

57

the changes in the vigor of thermal convection with planetary mass will be discussed in a much more58

quantitative manner than in Miyagoshi et al. (2018). We will also discuss the effects of model geometries

59

on the compressible convection, by carrying out calculations with identical physical parameters both in60

two-dimensional rectangular box and cylindrical annulus.

361

Model descriptions

62

We consider a time-dependent thermal convection of a highly compressible fluid of infinite Prandtl

63

number, as a model of mantle convection of massive super-Earths. The numerical model is almost64

the same as those of our earlier ones (Kameyama and Yamamoto 2018; Kameyama 2022) except that65∗
the mass M ∗ of modeled planets varies up to 10 times the Earth’s mass MEarth(Miyagoshi et al. 2018).

66(The asterisks stand for dimensional quantities throughout this paper.) In this study, we assume that67∗
the thickness d∗M of the planetary mantle and the gravitational acceleration gM
depend on M ∗ as

0.28

d∗M /d∗Earth = (M/MEarth )

0.5

∗
∗
gM
/gEarth= (M/MEarth )

,

,

(1)

68

∗
∗
(where d∗Earth and gEarth
are the values for M ∗ = MEarth) according to the scaling relationships proposed69

by Valencia et al. (2006).

70

The convecting vessel is taken to be either a two-dimensional rectangular or cylindrical geometries. For71

the rectangular models, a Cartesian box of height d∗M and an aspect ratio (width/height) 4 is employed.

72

For the cylindrical models, on the other hand, a half (azimuth θmin ≤ θ ≤ θmax where θmax − θmin = π) of

73

∗
∗
∗
∗
∗
∗
cylindrical annulus (radius rmin
≤ r∗ ≤ rmax
where rmax
− rmin= d∗M and rmin
/rmax
= 0.5) is considered.

74

∗
∗
∗
) at the top and bottom
and Tbot(> Ttop
The temperature T ∗ of the layer is fixed to be T ∗ = Ttop

75boundaries, respectively, while the vertical side walls are insulating (zero heat flux). The impermeable76

and shear-stress-free conditions are adopted along all boundaries. In order to keep the present analysis

77

as simple as possible, we ignored the effects of internal heating and spatial variations in specific heat78∗/∆T ∗ = 0.1 (whereCp∗ of the fluid and the gravity g ∗ throughout the layer. In addition we fixed Ttop

79∗∗
∆T ∗ ≡ Tbot
− Ttopis the temperature difference across the layer) as is expected in the Earth’s mantle,

80

in order to reduce the numbers of varying parameters.

81

Physical properties of modeled compressible fluid82

As a model of pressure-dependence in thermodynamic properties of mantle materials of super-Earths

83

whose mass is M ∗ , we employ the similar variations in thermal expansivity α∗ and reference density ρ∗

84with depth as in our earlier work (Miyagoshi et al. 2018). In Figure 1 we schematically show their vertical85

profiles. We first assume that the hydrostatic pressure p∗ and the reference density ρ∗ linearly increase86

with the depth y ∗ from the top surface asp∗ = p∗bot (y ∗ /d∗M ) ,ρ∗ = ρ∗top + ρ∗bot (y ∗ /d∗M ) .

4(2)

Table 1. The pressure dependence of the thermal expansivity α∗ assumed in the model.
Thermal expansivity (10−5 K−1 )Pressure (GPa)
0

4.0

10

3.5

30

2.6

60

1.86

100

1.38

150

1.10

300

0.701

600

0.37641000

0.3034

1500

0.258

(a) p∗

(b) α
0

−1

−1

−1

depth [×d∗Earth ]

0

−2

−2
0

Figure 1.

(c) ρ

0

[GPa]

1500

0

∗
α∗ /αtop

−2
1

1.0

ρ∗ /ρ∗top

3.5

Plots of the variations in (a) hydrostatic pressure p∗ , (b) thermal expansivity α and (c)reference density ρ as functions of the depth from the top surfaces for the ratio of the planetary mass∗
M ∗ to the Earth’s MEarthin the range of 1 ≤ M/MEarth ≤ 10.

587Here p∗bot is the hydrostatic pressure at the bottom boundary (y ∗ = d∗M ) and given byp∗bot [GPa] = 7.241 + 127.4 × (M/MEarth ),

(3)

88

while the values of ρ∗top and ρ∗bot are the values of reference density at the top (y ∗ = 0) and bottom

89boundaries and taken to satisfy
ρ∗bot /ρ∗top = 1.75 × (M/MEarth )0.258 .

(4)

90

For M/MEarth = 10, for example, p∗bot is about 1.3 terapascals and ρ∗bot is larger than ρ∗top by about 3.17

91

(Valencia et al. 2006). On the other hand, the thermal expansivity α∗ decreases with depth from the top

92

surface along with the increase in the hydrostatic pressure p∗ , whose depth profile is calculated by a linear

93

interpolation from the variation in α∗ with pressure shown in Table 1. The estimate of α∗ was made94

for MgO from an ab-initio calculation (Tsuchiya and Tsuchiya 2011). We note that the linear increase95

in ρ∗ assumed here is not likely to spoil the significance of the present analysis despite its geophysical96

unsoundness, since the vertical profile of (reference) density has a very minor effect on the nature of97

critical states of thermal convection of highly compressible fluids (Kameyama et al. 2015).

98

In addition to the variations in thermodynamic properties (α∗ and ρ∗ ), we also assume that the viscosity99η ∗ depends on temperature T ∗ asT∗ − T∗∗η ∗ (T ∗ ) = ηbot
exp ln(rT ) × ∗bot,∗
Tbot − Ttop

(5)

100

∗
∗
where ηbotis the viscosity at the bottom boundary (T ∗ = Tbot), and rT is a dimensionless constant101∗
controlling the temperature-dependence in viscosity. This equation means that the viscosity ηtop
at the

102

∗
∗
∗
∗
top surface (T ∗ = Ttop
) satisfies rT = ηtop/ηbot
. In this study we assumed that ηbotand rT = 105 are kept

103

unchanged regardless of the planetary mass M ∗ . The value of rT is chosen so that the stiff lithosphere104∗
develops along the cold surface of the planets. By assuming that ηbotis independent of M ∗ , we will105

exclude the cases with a huge increase in viscosity in the deep mantle of the super-Earths proposed in106

several earlier studies (Stamenković et al. 2011, 2012). Indeed, a recent theoretical work by Ritterbex

107

et al. (2018) predicts a substantial decrease in viscosity in the mantles of massive super-Earths associated108with the B1-B2 transition of MgO at around 500 GPa. On the other hand, the thermal conductivity k ∗109

is assumed be uniform in the entire layer and independent of M .

6

110

Basic equations

111

Under the truncated anelastic liquid approximation (TALA), the basic equations can be written in a112

dimensionless form as:
0 = ∇ · (ρv) ,(6)1(∇ ⊗ v + v ⊗ ∇) ,
2
 

1
0 = −∇p + 2η ε̇ − (∇ · v)I· ∇ + Raρ αT eup ,
3

ε̇ =

ρ

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

(7)
(8)
(9)

113

where the nondimensionalization is carried out with the length scale of d∗M (thickness of planetary mantle),

114

time scale of (d∗M )2 /κ∗top (where κ∗top = k ∗ /(ρ∗top Cp∗ ) is the value of thermal diffusivity at the top surface),

115

∗
∗
velocity scale of κ∗top /d∗M and temperature scale of ∆T ∗ ≡ Tbot
− Ttop. (Note that the scales of time116

and velocity depend on the planetary mass M ∗ because the length scale d∗M depends on M ∗ .) Here v is117

the velocity vector, ε̇ is the strain rate tensor, ⊗ stands for the tensor product of vectors, v ⊗ ∇ is the118

transpose of ∇ ⊗ v, I is the identity tensor, eup is the unit vector in the upward direction (eup = ez for

119

rectangular box, eup = er for cylindrical annulus), p is pressure, vup ≡ (v · eup ) is the velocity in the120

upward direction, and ε̇II is the second invariant of deviatoric strain rate tensor. In the right-hand side121of (9), the first, second, and third terms represent the effects of apparent heat transport (advection plus

122

conduction), adiabatic temperature change, and viscous dissipation (frictional heating), respectively.

123

We note here the validity of the truncated anelastic liquid approximation (TALA) in the present numerical124

experiments where the effects on buoyancy of dynamic pressure (the deviation of pressure from the125hydrostatic one) are neglected. Since the formulation based on TALA is known to introduce an imbalance

126

between viscous and adiabatic heating (Leng and Zhong 2008), one may think that it can severely violate

127

the energy budget in the entire domain particularly under the strong compression (King et al. 2010).

128

However, it is not likely the case for the mantle convection of massive terrestrial planets. As had been129

already demonstrated in our earlier works (Kameyama and Yamamoto 2018; Kameyama 2022), the130

mismatch between the input and output of heat at the top and bottom boundaries is only up to a few

131

per cent. We therefore believe that the present modeling based on TALA can be sufficiently valid for the132

convecting flows even in the mantles of massive super-Earths since the energy balance almost holds in133

the entire domain, although the balance is expected to be improved by using better approximations such134

as its untruncated counterpart (ALA) (Leng and Zhong 2008; Liu and Zhong 2013) or by solving fully7

135

compressible equations without approximations (Ricard et al. 2022; Ricard and Alboussière 2023).

136

The present formulation includes two nondimensional parameters which vary depending on the planetary

137

mass M ∗ . First one is the dissipation number Di of the modeled fluid defined by,
Di ≡

138

∗
∗ ∗
αtop
gM
dM
.
∗
Cp

(10)

Second one is the Rayleigh number Ra of thermal convection defined by,
Ra ≡

∗
∗
∆T ∗ gM(d∗M )3
ρ∗top αtop.
∗∗ηbot κtop

(11)

139

∗
Note here that, as indicated by (5), the viscosity is scaled by the value ηbotat the bottom boundary

140

while other properties are done by those at top surface. The values of Di and Ra for the mantles of141

terrestrial planets are expected to increase with the mass M ∗ of the planets, coming from the increase142∗
in the gravity gMand the mantle thickness d∗M of the terrestrial planets whose composition is similar to143

the Earth’s (e.g., Seager et al. 2007; Valencia et al. 2007; Wagner et al. 2011). According to the scaling144

relationship given by (1), Di and Ra increase with M ∗ as
0.78

Di/DiEarth = (M/MEarth ),Ra/RaEarth = (M/MEarth )

1.34

.

(12)

145

In this study we assumed Di = DiEarth = 0.83 and Ra = RaEarth = 4.6 × 107 for M/MEarth = 1

146while Di = 5.0 and Ra = 1.0 × 109 for M/MEarth = 10. (Tachinami et al. 2014; Miyagoshi et al. 2014,

147

2015; Kameyama and Yamamoto 2018; Kameyama 2022). (Note that the values of Di and Ra for given

148

M employed here are slightly different from those in Miyagoshi et al. (2018), because we stick to the149

consistency with the choices for M/MEarth = 10 in our earlier studies.) In the later section, however, we

150

will consider “effective” dissipation number and Rayleigh number in addition to Di and Ra defined here,151

as appropriate measures of the vigor of thermal convection of compressible fluids.

152

The employed initial condition for temperature T is taken to be equal to the steady one-dimensional

153

conductive state in the vertical direction. We added a small perturbation in temperature on the order154

of O(10−2 ) in the dimensionless unit which induces a single convection cell in the model domain. We

155

continue the calculation till the average temperature in the convecting vessel and the heat flow on both156

the surface and bottom boundaries reach their statistically steady state.

157

Numerical method158

The basic equations are discretized by the finite-difference method based on the control-volume method

159(Patankar 1980). A uniform mesh is employed with 256 and 1024 mesh points in vertical and horizontal
8

160

directions, respectively. We confirmed that the thermal boundary layer (TBL) along the bottom boundary161

is spanned by at least three mesh points excluding the boundary point.

162

is discretized by the Crank-Nicholson method in time. An upwind scheme, called power-law scheme

163(Patankar 1980), is used to evaluate the advection and conduction terms in the energy equation. The

164

The energy equation (9)equations (6) and (8) for flow field are solved for a “potential” ψ of the mass flux ρv defined by,


∂ψ
∂ψ


ρvx = −, ρvz =
for 2-D rectangular box
∂z
∂x



ρvr = 1 ∂ψ , ρvθ = − ∂ψ for 2-D cylindrical annulus
r ∂θ
∂r

(13)

165

instead of directly solving for primitive variables (velocity v and pressure p). Note that the function166

ψ is a generalization of the stream function for the Boussinesq or extended Boussinesq approximation

167

where ρ = 1 throughout the layer. As in the case with incompressible fluids, the velocity field v defined168

above automatically satisfies the equation of continuity (6). The reliability of the numerical codes has

169

been already verified (Kameyama and Harada 2017; Kameyama and Yamamoto 2018; Tsuchida and170

Kameyama 2020; Kameyama 2022). See Kameyama and Yamamoto (2018) and Kameyama (2022) for171

the explicit forms of basic equations and our numerical technique for the solution of flow field.

172

Results and Discussion

173

We carried out numerical experiments of thermal convection of highly compressible fluids both in the174

two-dimensional rectangular box and cylindrical annulus by varying the mass of the modeled planets M

175

in the range of 1 ≤ M/MEarth ≤ 10 where MEarth is the Earth’s mass. The snapshots of the overall

176

thermal convection for all the cases calculated in this study can be found in Figs S1 to S4. For the177sake of conciseness, we will discuss in the following the convecting flow patterns mainly for the cases in

178

two-dimensional cylindrical annulus.

179

Estimates on the effects of depth-dependent thermal exansivity on compressible convection180

As had been introduced earlier, we assumed that the thermal expansivity α depends on depth y from181

the top surface, and that its variation with depth also depends on the mass M of modeled planets (see182

Figure 1b). Before presenting our numerical results with varying M , we will discuss in this subsection

183

how the change in α with depth can affect the thermal convection of compressible fluids.

184

Let us first estimate the average thermal expansivity αave whose dimensional value is calculated by
∗
αave
=

Z y∗ =d∗My ∗ =09α∗ (y ∗ )dy ∗ ,(14)(a) αbot /αtop and αave /αtop
0.6
0.5
0.4

∝(

(b) ∆Tpot /∆T
0.7

M/
M

0.3

∝(

0.6

Ea
rth ) −

0.5

2

M/
M

Ea
rt

0.5

h ) − 0.

0.2

19

0.4
0.1
0.09
0.08
0.07
0.06

● αbot /αtop★ αave /αtop

1

0.3

2

5

10

∆Ttop /∆T = 0.1
12f(c) Di and Di
5
4
3

109

Di
f
★ Di
●

●

∝

1
0.9
0.8
0.7
0.6
0.5
0.4

(M

/M

∝ (M

rt
Ea

10

Ra

f
★ Ra

8
0.7

2

5

f
(d) Ra and Ra
4

1 .3

)
h

t

ar

108

th
/M Ear

)
h

)0

∝

.26

(M

E
/M

∝ (M

M Earth

/

)0

.63

107
1

2

5

10

1

M/MEarth2

5

10

M/MEarth

Figure 2. (a) The plots of the ratio of the thermal expansivity αbot at the bottom boundary to αtop (redcircles) and the ratio of the average thermal expansivity αave in the entire layer to αtop (black stars)
against the ratio of the planetary mass M to the Earth’s one MEarth . (b) The plot against M/MEarthof the difference ∆Tpot in the potential temperature across the layer. (c) The plots against M/MEarthf (black stars). (d) The plots against M/MEarth of theof the dissipation number Di (red circles) and Di
f (black stars). Shown by the solid gray lines in (a) to (d) are
Rayleigh number Ra (red circles) and Ra
the relationships fitted by the powers of M/MEarth .

10

185

where the spatial integration is taken from the top surface (y ∗ = 0) to the bottom boundary (y ∗ = d∗M ).

186(Note that the dimensional depth d∗M of the mantle depends on M ∗ .) In Figure 2a we plot the ratio αave

187

to αtop against the ratio of planetary mass M to the Earth’s one MEarth calculated from the profiles of α

188

shown in Figure 1b, together with those of αbot to αtop . As can be seen in the plot, αave /αtop is smaller189

for larger planets, and decreases approximately in proportion to (M/MEarth )−0.52 .

190

We will next estimate how the adiabatic temperature change is affected by the variation in α with depth.

191

To this purpose we consider the potential temperature Tpot whose dimensionless value at dimensionless

192

depth y is calculated by
Z y′ =y"
Tpot = T exp −Di

y ′ =0#α(y ′ ) ′
dy ,
αtop

(15)

193

where y = 1 − z for rectangular and y = rmax − r for cylindrical cases. The physical meaning of Tpot

194

can be clearly understood from the adiabatic temperature gradient (dT /dy)s whose dimensionless form

195

is written as


dT
dy


=
s

∗ ∗
dM
α∗ gM
α
T = Di
T.
Cp∗
αtop

(16)

196

That is, the potential temperature Tpot of a fluid parcel located at depth y is the temperature of the197

parcel after it would ascend to the top surface (y = 0) in an adiabatic manner. In particular, the potential

198

temperature Tpot,bot at the bottom boundary (y = 1) can be given by"

Z y′ =1

Tpot,bot = T (y = 1) exp −Di
y ′ =0#αave
α(y ′ ) ′dy = Tbot exp −Di
.
αtop
αtop

(17)

199

f ≡ Di(αave /αtop ). In Figure 2c we plot theWe will hereafter denote an “effective” dissipation number Di

200

f against M/MEarth . By combining Di ∝ (M/MEarth )0.78 and αave ∝ (M/MEarth )−0.52
values of Di and Di

201

f ∝ (M/MEarth )0.26 . That is, owing to the depth-dependence in thermal expansivity α∗ , thewe obtain Di

202

f in the entire layer is much smaller than Di given by (10) using the“effective” dissipation number Di

203

∗
f
value of thermal expansivity αtop
at the top surface. For M/MEarth = 10, for example, the value of Di

204

is approximately 0.79, which is close to the value assumed in Tackley et al. (2013).

205

The physical meaning of the potential temperature Tpot further implies that the distribution of Tpot is

206

more appropriate than that of the “raw” temperature T for discussing the buoyancy force that drives

207

thermal convection, since the temperature change due to adiabatic (de)compression is subtracted. Let

208

us therefore consider a “super-adiabatic” temperature difference ∆Tpot across the entire layer, which

209

dictates the difference in Tpot between the top and bottom boundaries. By definition, ∆Tpot is calculated

11

210

by


h

i
f − Ttop = ∆T − (Ttop + ∆T ) 1 − exp −Di
f .
∆Tpot ≡ Tpot,bot − Ttop = Tbot exp −Di

(18)

211

f and ∆T but also on Ttop . In particular, ∆Tpot becomes(Note here that ∆Tpot depends not only on Di

212

smaller for larger Ttop .) In Figure 2b we plot the ratio of ∆Tpot to ∆T = Tbot − Ttop , calculated with the213

variation of α with depth shown in Figure 1b and Ttop /∆T = 0.1 assumed in this study. As can be seen214

in the plot, ∆Tpot /∆T is smaller for larger planets, and its value decreases approximately in proportion215

to (M/MEarth )−0.19 .

216

f defined
From the discussion so far in this subsection, we can consider an “effective” Rayleigh number Ra

217

by
f ≡
Ra

∗
∗
ρ∗top αave
∆Tpot ∗ gM
(d∗M )3
αave ∆Tpot
= Ra
,
∗
∗
ηbot κtop
αtop ∆T(19)

218

in addition to Ra given in (11). Taken the relationships of αave /αtop and ∆Tpot /∆T with M into account,219

f varies with the planetary mass in proportion to (M/MEarth )0.63 (see also Figure 2d). As will be seenRa

220

later, the change in the vigor of thermal convection of compressible fluid can be estimated much more

221

f rather than by Ra.
appropriately by Ra

222

Changes in the structures of thermal convection with planetary mass

223

In Figures 3 and 4 we present the overall convection structures calculated in two-dimensional cylindrical224

annulus using several values of M/MEarth . Note in these figures that the thickness dM of convecting layer

225

differs between the cases depending on the mass M of modeled planets (see Equation (1)). Similarly, both226

the Rayleigh number Ra and the dissipation number Di increase with M as indicated in the figures. In

227

Figure 3 we show (a) the distributions of temperature T and (b) the plots against the dimensionless height

228

from the bottom boundary (r − rmin ) of the horizontally averaged temperature ⟨T ⟩ at r. Also shown by

229

thin solid lines in Figure 3b are the radial profiles of temperatures given by several adiabats (i.e., the lines

230

of constant Tpot ). In Figure 4, on the other hand, we show (a) the distributions of potential temperature

231

Tpot and (b) the plots of the horizontally averaged potential temperature ⟨Tpot ⟩ at r, calculated for232

the same snapshots in Figure 3. Also indicated by the shaded regions in Figure 4b are the ranges of233|Tpot − ⟨Tpot ⟩| ≤ δTpot (r), where δTpot is given bys
Z θ=θmax
1
2[Tpot (r, θ) − ⟨Tpot ⟩(r)] dθ.
δTpot ≡
θmax − θmin θ=θmin

(20)

234

That is, δTpot is the standard deviation of the distribution of potential temperature in the horizontal

235

plane at r.
12

r − rmin(a) temperature T(b) ⟨T ⟩M/MEarth = 1 (Ra = 4.6 × 107 , Di = 0.83)
1
1

0
0.1
1.1
0.1
1.1
M/MEarth = 2 (Ra = 1.2 × 108 , Di = 1.42)
1
1r − rmin

0

0
0.1
1.1
0.1
1.1
M/MEarth = 4 (Ra = 2.9 × 108 , Di = 2.45)
1
1

r − rmin

0

0
0.1
1.1
0.1
1.1
M/MEarth = 6 (Ra = 5.0 × 108 , Di = 3.36)
1
1

r − rmin

0

0
0.1
1.1
0.1
1.1
M/MEarth = 8 (Ra = 7.4 × 108 , Di = 4.20)
1
1

r − rmin

0

0
0.1
0.1
0.1

Figure 3.

1.1
1.1

0

1.1

Comparison of overall structures of the thermal convection calculated for the cases in

two-dimensional cylindrical annulus with various values of the ratio of planetary mass M to the Earth’s
mass MEarth . (a) The distributions of dimensionless temperature T . The contour interval is 0.1. The
color scale is indicated at the bottom of the figure. (b) Plots against dimensionless height from thebottom boundary of the horizontally averaged temperature ⟨T ⟩ (thick red lines) at the height of r.
Indicated by the thin solid lines for comparison are the plots of several adiabats against r.

13

(a) Tpot(b) ⟨Tpot ⟩M/MEarth = 1 (Ra = 4.6 × 107 , Di = 0.83)
1
1

r − rmin

0.2
0.3
0.4
0.5
0.

6

0.7

0
0.1
1.1
0.1
1.1
M/MEarth = 2 (Ra = 1.2 × 108 , Di = 1.42)
1
10

r − rmin

0.2
0.3
0.4
0.5

0.6

0.6

0.6

0
0.1
1.1
0.1
1.1
M/MEarth = 4 (Ra = 2.9 × 108 , Di = 2.45)
1
1
0

0.5

r − rmin

0.2
0.3
0.4
0.5

0
0.1
1.1
0.1
1.1
M/MEarth = 6 (Ra = 5.0 × 108 , Di = 3.36)
1
1
0

4
0.5

r − rmin

0.2
0.3
0.

0.5

0.5

0
0.1
1.1
0.1
1.1
M/MEarth = 8 (Ra = 7.4 × 108 , Di = 4.20)
1
1
0

r − rmin

0.2
0.3
0.4

0.5

0
0.1
0.1
0.10

Figure 4.

1.1
1.1

0

0.72

Comparison of overall structures of the thermal convection calculated for the cases in

two-dimensional cylindrical annulus with various values of the ratio of planetary mass M to the Earth’s
mass MEarth . (a) The distributions of dimensionless potential temperature Tpot . The contour interval
is 0.05. The color scale is indicated at the bottom of the figure. (b) Plots against dimensionless heightfrom the bottom boundary of the horizontally averaged potential temperature ⟨Tpot ⟩ (thick red lines)
at the height of r. Delineated by the shaded regions in the plots are the ranges where |Tpot − ⟨Tpot ⟩| issmaller than the standard deviation of Tpot at given r. The vertical dashed lines indicate the values ofpotential temperature Tpot,ic in the isentropic core, while the horizontal dotted lines indicate the levelsof r at the bottom of the cold thermal boundary layers.
14

236

From the snapshots in Figures 3 and 4 we can clearly see that the thermal convection in the layer has

237

a very wide convection cell characterized by (i) a cold thermal boundary layer (TBL) whose thickness

238

increases from the left to the right corners of the box along the top surface (see the slope of the isotherms239

in the cold TBL), (ii) a major descending plume along the right-side wall, and (iii) a thin hot TBL along240

the bottom boundary from which minor ascending plumes are originating (see also Figures S1 to S4). The

241

fluid in the cold TBL has a much higher viscosity than the underlying fluid owing to the temperature-242

dependence in viscosity (rT ≡ ηtop /ηbot = 105 ), indicating that the cold TBL acts as a stiff lid similar243

to the lithosphere of terrestrial planets. In addition, there occurs a region at mid-depth between the top244

and bottom TBLs where the temperature T increases with depth along with an adiabat (see Figure 3a245

and 3b). This region, hereafter denoted by “isentropic core”, can be also seen as a broad region where the246potential temperature Tpot is almost uniform (see Figure 4a and 4b). A similar convecting flow structure247

is obtained for all the cases calculated in this study (see Figures S1 to S4).

248

However, from a detailed comparison between the cases, we can clearly notice the differences in the249

convection structures induced by the changes in mass M of the modeled planets. Indeed, the potential

250

temperature becomes lower in the deep mantle of larger planets with increasing M . Of course, this comes251

from the effect of adiabatic compression. A stronger adiabatic temperature change occurs during the252

vertical motion of fluid parcels, and it reduces the differences in the potential temperature across the253

mantle to a larger extent.

254

In order to discuss in a much more quantitative manner the effects of adiabatic compression on the255

thermal structure at depth, we plot against the planetary mass M in Figure 5a the values of the potential

256

temperature Tpot,ic of the “isentropic cores” for all the cases calculated in this study. Here the values

257

of Tpot,ic is estimated from those of ⟨Tpot ⟩ at the height where the standard deviation δTpot of potential

258

temperature is the local minimum. (In Figures 4b and S1d to S4d, we show the derived values of Tpot,ic259

by the vertical dashed lines.) The red solid circles and blue solid squares stand for the values obtained

260

for the cases in two-dimensional cylindrical and rectangular cases, respectively. Also shown by the black

261

stars for comparison are the values of the potential temperature Tpot,bot at the bottom boundaries given262

by Equation (17). From these plots we can first confirm the earlier finding that both the potential

263

temperatures at the bottom boundary and in the isentropic regions decrease with the planetary mass M .

264

In addition, the comparison between the different model geometries shows that for given M the values of265

Tpot,ic are higher in the rectangular cases than those in cylindrical cases. The difference in Tpot,ic due to15

(a) Tpot,ic

0.75(b) ∆Tpot,btbl

0.06

cylindrical■ rectangular★ Tpot,bot
●

0.70
0.65

●
■

cylindrical
rectangular

0.60
0.55
0.50
0.45

0.00
0

2

4 6 8
M/MEarth

10

0

2

4 6 8
M/MEarth

10

Figure 5. (a) Plots against the ratio of the planetary mass M to the Earth’s mass MEarth of the potential
temperature Tpot,ic of the isentropic cores calculated for the cases in two-dimensional cylindrical annulus(red solid circles) and in two-dimensional rectangular box (blue solid squares). Also plotted by solid
black stars for comparison are the values of potential temperature Tpot,bot at the bottom boundary. (b)Plots against M/MEarth of the difference in potential temperature ∆Tpot,btbl ≡ Tpot,bot − Tpot,ic acrossthe bottom thermal boundary layers.

16

Table 2. Summaries of the scaling relationships with the planetary mass M ∗ of the (a) dimensionlessroot-mean-square velocity vrms of convection and (b) dimensionless thickness ℓttbl of the top cold TBL
indicated by gray solid lines or black dotted lines in Figure 6.

f ∝ (M/MEarth )0.63
Ra(a) vrms ∝ Ra2/3(b) ℓttbl ∝ Ra−1/3

vrms ∝ (M/MEarth )−0.42ℓttbl ∝ (M/MEarth )−0.21

vrms ∝ (M/MEarth )−0.893ℓttbl ∝ (M/MEarth )−0.447(gray solid lines)
Ra ∝ (M/MEarth )1.34(black dotted lines)266

the model geometries comes from the difference in the changes in the surface area with depth. The surface

267

area decreases with depth in proportion to r−1 for the cylindrical cases, while it is independent of depth

268

for the rectangular cases. Owing to this effect, the magnitude of the vertical gradient in temperature T

269

in the shallow part of fluid layer is smaller in the cylindrical cases cases than in rectangular cases and, at270

the same time, the magnitude in the deep part is larger. This results in lower Tpot,ic in the entire layer271

for cylindrical cases than in rectangular cases.

272

What should be also noted in Figure 5a is that the difference between the values of Tpot,ic and Tpot,bot

273

decreases with increasing M for both geometries. To see this more clearly, we plot in Figure 5b the values274of ∆Tpot,btbl ≡ Tpot,bot − Tpot,ic against M . The plots in Figure 5b clearly show that in general the values

275

of ∆Tpot,btbl decreases with M for both geometries. Considering that ∆Tpot,btbl represents the difference276

in the potential temperatures which is supported by the bottom TBLs, the decrease in ∆Tpot,btbl with277

M is a clear manifestation that the vigor of ascending plumes is reduced with M . In addition, the values

278of ∆Tpot,btbl for given M are larger in the cylindrical cases than those in rectangular cases, which can be279

intuitively understood from the effect of model geometries. Moreover, the plots in Figure 5b clearly show

280

that the value of ∆Tpot,btbl rapidly decreases with M for M/MEarth ≲ 3 while it decreases very mildly

281

for M/MEarth ≳ 3. From this we can conclude that the ascending plumes become weak from the bottom282

boundary when the mass M of modeled terrestrial planets is larger than about 3MEarth .

283

Scaling relationships of the vigor of thermal convection with planetary mass

284

In order to quantitatively see the effects of planetary mass on the vigor of overall thermal convection, we

285

next plot against M/MEarth in Figure 6a the dimensionless root-mean-square velocity vrms of convection286

in 2-D cylindrical annulus (red solid circles) and in 2-D rectangular box (blue solid squares). Also17

(a) vrms(b) N u500(c) ℓttbl

10

0.6

400
300
200
100

●
■

0
0

2

cylindrical
rectangular

4 6 8
M/MEarth

10

●
■

0
0

2

cylindrical
rectangular

4 6 8
M/MEarth

10

●
■

0.0
0

2

cylindrical
rectangular

4 6 8
M/MEarth

Figure 6. (a) Plots against the ratio of the planetary mass M to the Earth’s mass MEarth of the dimensionless root-mean-square velocity vrms of convection calculated for the cases in two-dimensional cylindricalannulus (red solid circles) and in two-dimensional rectangular box (blue solid squares). The gray solidand black dotted lines indicate the relationships of vrms ∝ (M/MEarth )0.42 and vrms ∝ (M/MEarth )0.893 ,respectively. (see main text and Table 2 for detail). (b) Plots against M/MEarth of the Nusselt number
N u. The gray solid and black dotted lines indicate the relationships of N u ∝ (M/MEarth )0.21 andN u ∝ (M/MEarth )0.447 , respectively. (c) Plots against M/MEarth of the dimensionless thickness ℓttbl
of the top thermal boundary layer. The gray solid and black dotted lines indicate the relationshipsof ℓttbl ∝ (M/MEarth )−0.21 and ℓttbl ∝ (M/MEarth )−0.447 , respectively. Also shown by the black solidlines in (b) and (c) are the relationships of N u ∝ (M/MEarth )0.163 and ℓttbl ∝ (M/MEarth )−0.163 , respectively.

18

10

287

shown by the gray solid and black dotted lines in Figure 6a for comparison are the relationships between

288

the planetary mass M and vrms based on the well-known Boussinesq results where vrms is proportional289

to the Rayleigh number to the power of two-thirds (e.g., Turcotte and Schubert 2014). In order to290further derive the relationships between these quantities with the planetary mass M , we substituted291

the relationships between the Rayleigh number and M obtained earlier in Figure 2d. The gray solid292

f givenline in Figure 6a indicates the expected relationship using the “effective” Rayleigh number Ra

293

by Equation (19), while the gray dotted line indicates that using the original Ra defined by Equation

294(11). (In Table 2 we summarize the relationships obtained thus and shown in Figure 6.) The plots in295

Figure 6a clearly show that the variation in vrms with M obtained in numerical experiments follows very296

f indicated by gray solid lines.
closely the expected relationship using the effective Rayleigh number Ra

297

In contrast, the relationship using the original Rayleigh number Ra indicated by black dotted lines can298

hardly explain the numerical results. From this comparison we can conclude that the overall vigor of299

fthermal convection of compressible fluids is precisely measured with the effective Rayleigh number Ra,

300

which appropriately incorporates the effects of the decrease in thermal expansivity with depth as well as

301

those of the superadiabatic temperature difference to drive the convection.

302

Similar conclusions can be drawn also from the plots in Figure 6b and 6c where we show the variations in

303

the Nusselt number N u and the dimensionless thickness ℓttbl of the top cold TBL with M , respectively.

304

Here N u is defined by





∂T


−

∂z z=1
Nu =




1
∂T


−
rmax ln(rmax /rmin ) ∂r r=rmax

for 2-D rectangular box(21)
for 2-D cylindrical annulus

305

in dimensionless units. On the other hand, the base of the top cold TBL is taken to be the depth level306

where ⟨Tpot ⟩ − Ttop becomes 98% of Tpot,ic − Ttop , and is indicated by the horizontal dotted lines in Figure307

4b (see also Figures S1b to S4b). The plots in Figure 6b and 6c show that N u increases with M while308ℓttbl decreases with M . This is again in general agreement with the increase in the Rayleigh number with309

M . We also show by the gray solid lines and black dotted lines for comparison the relationships between

310

the planetary mass M and N u or ℓttbl expected from the relationships where N u is proportional to the311

Rayleigh number to the power of one-third (e.g., Turcotte and Schubert 2014) together with ℓttbl ∝ N u−1 .

312

f while the black dotted lines indicate
The gray solid lines indicate the relationships expected from Ra,

313

those from original Ra. The plots in Figure 6b and 6c clearly show that the variations in N u and ℓttbl
19

314

with M obtained for the cylindrical cases follow much more closely the expected relationship using the315

f indicated by gray solid lines than those using the original Rayleigh number
effective Rayleigh number Ra

316

Ra indicated by black dotted lines.

317

We also show by the thin black solid lines in Figure 6b and 6c the relationships between the planetary

318

mass M and N u or ℓttbl expected from the relationships where N u is proportional to the Rayleigh319

number to the power of 0.258 based on the compilation of the experiments of Boussinesq convection with320

temperature-dependent viscosity in 2-D Cartesian domain (e.g., Deschamps and Sotin 2000). From the321plots in Figure 6b and 6c we can clearly see that for 2-D Cartesian cases the variations in N u and ℓttbl with322f
M follow very closely the relationship based on N u ∝ Ra

323

indicated by the gray solid lines. The better agreement with the former relationship is probably due to324

the effect of model geometries which affects the overall thermal states and, in turn, an effective viscosity

325

in the convecting layers.

326

We finally note that the present results are not very consistent with those of our earlier study (Miyagoshi

327

et al. 2018) where ℓttbl was almost independent with M . The discrepancy between the present study

328

and earlier one is probably due to the difference in the numbers of convection cells in the vessels. The

329

number of convection cell is always one for all the cases in this study (see also Figures S1 to S4), while330

it varies with M in Miyagoshi et al. (2018) by allowing multiple cells for small M . In other words, the331

decrease in ℓttbl with increasing M is most likely to be overshadowed in Miyagoshi et al. (2018), owing

332

to the occurrence of narrow convection cells for the cases with small M .

0.258

333

1/3

f
, rather than that based on N u ∝ Ra

Feasibility of plate tectonics on larger terrestrial planets?

334

In this subsection we will revisit the issue of plate tectonics on super-Earths larger than the Earth. In

335

earlier literature (e.g., Valencia et al. 2007; Valencia and O’Connell 2009; van Heck and Tackley 2011;

336

Miyagoshi et al. 2018), the stress induced in the lithosphere by mantle convection has been estimated and337

compared with the rupture strength of the lithosphere to see if plate tectonics operates in super-Earths.

338

The key in their estimates lies in the assumption that the generated stress is inversely proportional to339

the thickness of the lithosphere along the top surface and is proportional to the convecting velocity in340

the underlying mantle. In the following we will make similar estimates based on our numerical results.

341∗
In Figure 7 we plot against M/MEarth (a) the dimensional root-mean-square velocity vems
in the entire

342

layer and (b) the dimensional thickness ℓ∗ttbl of the cold TBL along the top surface. In this figure we343

also show by the gray solid lines the relationships between M/MEarth and vrms or ℓttbl expected from the20

rectangular[mm/yr]■(b) ℓ∗ttbl (dimensional)
2● cylindrical2ℓ∗ttbl /ℓ∗ttbl,Earth∗(a) vrms(dimensional)
4● cylindrical0

rectangular100Figure 7.

■

2

4 6 8
M/MEarth

10

0

2

4 6 8
M/MEarth

10(a) Plots against the ratio of the planetary mass M to the Earth’s mass MEarth of the∗of convection calculated for the cases in two-dimensional
dimensional root-mean-square velocity vrms

cylindrical annulus (red solid circles) and in two-dimensional rectangular box (blue solid squares). (b)Plots against M/MEarth of the ratio of the dimensional thickness ℓ∗ttbl of the top thermal boundary
layer to ℓ∗ttbl,Earth ≡ ℓ∗ttbl (M = MEarth ) calculated for the cases in two-dimensional cylindrical annulus(red solid circles) and in two-dimensional rectangular box (blue solid squares).

In (a) and (b),∗the gray solid lines indicate the relationships based on our estimates (vrms∝ (M/MEarth )0.14 andℓ∗ttbl ∝ (M/MEarth )0.07 ), while the black dashed lines indicate the relationships proposed by Valencia∗
et al. (2007) (vrms∝ (M/MEarth )1.19 and ℓ∗ttbl ∝ (M/MEarth )−0.45 ).

21

344

relationships developed in the previous subsection (see Figure 6 and Table 2) together with the increase345∗
in the length scale d∗M with M/MEarth given by Equation (1). These plots clearly show that both vrms

346

and ℓ∗ttbl gradually increase with M/MEarth , whose changes can be approximately given by∗vrms= vrms (κ∗top /d∗M ) ∝ (M/MEarth )0.14 ,ℓ∗ttbl = ℓttbl d∗M ∝ (M/MEarth )0.07 .

(22)

347

∗
These estimates lead to (vrms/ℓ∗ttbl ) ∝ (M/MEarth )0.07 , indicating that the stress level in the lithosphere348

is almost independent of the planetary mass M for M/MEarth ≥ 1. This conclusion is in a stark contrast349

to that of Valencia et al. (2007); by employing the Boussinesq approximation Valencia et al. (2007)350∗
derived the relationships indicated by the gray dashed lines in Figure 7 (vrms∝ (M/MEarth )1.19 and351ℓ∗ttbl ∝ (M/MEarth )−0.45 ) and concluded that plate tectonics takes place more readily in larger terrestrial

352

planets. The difference between our conclusion and earlier ones (e.g., Valencia et al. 2007; van Heck and353

Tackley 2011) is a clear manifestation that the effects of adiabatic compression is of ultimate importance354

in understanding the tectonics and dynamics of the mantles of massive super-Earths.

355

We finally note that the difference in our conclusion and that of Miyagoshi et al. (2018) simply comes from356

the difference in the dependence of the dimensionless thickness ℓ∗ttbl of the cold TBL on the planetary

357

mass M reflecting the difference in the horizontal length scales of convection cells.

358

Conclusions

359

In this paper we carried out a series of numerical experiments of thermal convection of compressible fluids

360

with temperature-dependent viscosity either in two-dimensional rectangular box and cylindrical annulus,361

in order to study how the adiabatic compression and model geometries affect the mantle convection362

on super-Earths whose mass is up to 10 times the Earth’s. We have demonstrated that the adiabatic

363

compression affects more strongly the thermal convection in the mantles of larger terrestrial planets. This364

is because the adiabatic temperature change, which occurs during the ascent or descent of fluid parcels,

365

reduces the vigor of thermal convection. In addition, regardless of the model geometries, ascending plumes366

from the bottom boundary become significantly weaker for super-Earths whose mass is larger than the367

Earth’s by a factor of about 3, resulting in ascending plumes associated with very small temperature368

difference from the surroundings. We also found that the adiabatic compression also affects the overall369

vigor of thermal convection measured by the dimensionless convecting velocity and the efficiency of heat

370

transport. In particular, both the increase in the convecting velocity and the decrease in the thickness
22

371

of cold thermal boundary layers (TBLs) are more mild than those expected from the increase in the372

nominal Rayleigh number under the Boussinesq approximation. This in turn makes the stress levels373

in the cold TBLs almost constant regardless of the planetary mass and/or the model geometries. Our374

finding suggests that the occurrence of plate tectonics is not likely to be enhanced for massive terrestrial

375

planets larger than our Earth, as opposed to the prediction of incompressible convection (e.g., Valencia

376

et al. 2007; van Heck and Tackley 2011). Taken together with very rapid cooling of small planets, our377

speculation from the viewpoint of mantle dynamics is that Earth-sized terrestrial planets are most likely

378

to be “most habitable”.

379

The fact that qualitatively similar results were obtained both for two-dimensional rectangular and380

cylindrical geometries suggests that the overall nature of thermal convection of compressible fluids is381

not significantly affected by the model geometries. This may further imply that the findings in the382

present numerical experiments can be directly applied to the mantles of massive super-Earths in nature

383

which takes place in three-dimensional spherical shell geometry. In particular, ascending plumes from384

the bottom boundary become significantly weaker in the three-dimensional mantles of super-Earths with385∗
M ∗ ≳ 3MEarth, although larger temperature difference between the plumes and the surroundings can be386

expected owing to the effect of spherical geometry which reduces the surface area in proportion to r−2 .

387

What is also worth being noted in our findings is that the variation in the overall vigor of mantle

388

convection with the planetary mass is precisely estimated from the changes in the effective Rayleigh389

f defined by (19) rather than the nominal Ra. This is a strong manifestation of an ultimate
number Ra

390

importance of the adiabatic temperature change in the thermal convection of compressible fluids on the391

one hand. This indicates, on the other hand, that the scaling relationship developed for the Boussinesq

392(or incompressible) convection (e.g., Deschamps and Sotin 2000; Turcotte and Schubert 2014) can be also393

f which incorporates
used for compressible cases simply by replacing the nominal Ra with the effective Ra394

the effects of adiabatic temperature change (and variation in thermal expansivity with depth). We thus

395

expect that the extensive archives of mantle convection studies under the Boussinesq approximation

396

conducted so far can help enrich the understanding the mantle dynamics of super-Earths under the397

non-Boussinesq conditions.

398

Although we do believe that our numerical model captures the fundamental nature of the convecting

399

motions expected in the mantles of super-Earths larger than the Earth, the significance of our finding

400

should be verified in future by incorporating various complexities which may affect the dynamics of the
23

401

convecting mantles of massive terrestrial planets. For example, we should take into account the effects of

402

internal heating either by radiogenic elements or by tidal dissipation, which are expected to increase the403

temperature in the interior and in turn to enhance the effects of adiabatic temperature change. This will404

be particularly true for the thick mantles of massive super-Earths, whose thermal evolution (or cooling) is

405

expected to be extremely slow (Miyagoshi et al. 2017). In addition, we will need to incorporate more likely

406

configurations for the mantle convection of super-Earths, such as the influences of internal structures,

407

physical properties, and phase changes of mantle materials. Among them, the transport properties

408(viscosity and thermal conductivity) under extreme conditions is of great importance in the dynamics of409

the lowermost mantle, since these properties may significantly affect the thermal state at depth. Indeed,410

a stable thermal stratification can be formed above the core-mantle boundary of massive super-Earths

411

if the thermal conductivity is very high at depth (Kameyama and Yamamoto 2018; Kameyama 2022),

412

resulting in an absence in ascending plumes there.

413

We finally note that our rheology model needs to be improved particularly under shallow and cold414

conditions, in order to address the occurrence of plate tectonics in a much more appropriate manner.

415

Among the various models employed so far to reproduce the plate tectonics in the framework of mantle

416

convection (e.g., Moresi and Solomatov 1998; Bercovici 1998; Tackley 2000), our candidate is the stress-417

history-dependent rheology proposed by Ogawa (2003). As had been demonstrated in our earlier work418

with incompressible mantle convection in three-dimensional rectangular box (Miyagoshi et al. 2020), this419

rheology model induces a hysteresis in viscosity under certain conditions, which can be the keys to the420

coexistence of a “damaged” state for weak narrow plate boundaries and an “intact” state for stiff broad

421

plate interiors. We thus speculate that an extension of this model to compressible cases may help greatly

422

deepen the understanding of plate tectonics on super-Earths.

423

Declarations

424

Availability of data and materials

425

The data that support the findings of this study and the numerical programs used in this study the are426

available upon reasonable request to M.K. (email: kameyama.masanori.my@ehime-u.ac.jp).

427

Competing interests

428

The author declares no competing interests.
24

429

Funding

430

This work was supported by JSPS KAKENHI Grant Number 22K03702, a joint research program at

431

Geodynamics Research Center, Ehime University, the Earth Simulator of the Japan Agency for Marine-

432

Earth Science and Technology (JAMSTEC), “Exploratory Challenge on Post-K computer” (Elucidation433

of the Birth of Exoplanets [Second Earth] and the Environmental Variations of Planets in the Solar

434

System), and “Program for Promoting Research on the Supercomputer Fugaku” (Toward a unified view435

of the universe: from large scale structures to planets).

436

Authors’ contributions

437

M.K. developed the numerical codes, performed numerical simulations, analysed simulation data, and438

prepared the manuscript.

439

Acknowledgments

440

We thank Masaki Ogawa, Takehiro Miyagoshi, Takatoshi Yanagisawa, Haruhiko Dekura, and Shin-ichi

441

Takehiro for fruitful discussion.

442

References

443

Bercovici D (1998) Generation of plate tectonics from lithosphere-mantle flow and void-volatile444

self-lubrication. Earth Planet Sci Lett 154(1-4):139–151, DOI 10.1016/S0012-821X(97)00182-9445

van den Berg AP, Yuen DA, Umemoto K, Jacobs MHG, Wentzcovitch RM (2019) Mass-dependent446

dynamics of terrestrial exoplanets using ab initio mineral properties. Icarus 317:412–426,

447

DOI 10.1016/j.icarus.2018.08.016

448

Davies GF (2011) Mantle convection for geologists. Cambridge University Press, 232 pp.

449

Deschamps F, Sotin C (2000) Inversion of two-dimensional numerical convection experiments for a fluid

450with a strongly temperature-dependent viscosity. Geophys J Int 143(1):204–218,

451

DOI 10.1046/j.1365-246x.2000.00228.x

452

453

Foley BJ (2019) Habitability of Earth-like stagnant lid planets: Climate evolution and recovery fromsnowball states. Astrophys J 875(1):72, DOI 10.3847/1538-4357/ab0f31454

Foley BJ, Bercovici D, Landuyt W (2012) The conditions for plate tectonics on super-Earths:

455

Inferences from convection models with damage. Earth Planet Sci Lett 331-332:281–290,
25

456

457

458

459

460

461

462

463

DOI 10.1016/j.epsl.2012.03.028
Gaidos E, Conrad CP, Manga M, Hernlund J (2010) Thermodynamic limits on magnetodynamos inrocky exoplanets. Astrophys J 718(2):596–609, DOI 10.1088/0004-637X/718/2/596
van Heck HJ, Tackley PJ (2011) Plate tectonics on super-Earths: Equally or more likely than on Earth.
Earth Planet Sci Lett 310(3-4):252–261, DOI 10.1016/j.epsl.2011.07.029
Jarvis GT, McKenzie DP (1980) Convection in a compressible fluid with infinite Prandtl number. J
Fluid Mech 96(3):515–583, DOI 10.1017/S002211208000225X
Kameyama M (2016) Linear analysis on the onset of thermal convection of highly compressible fluids464with variable physical properties: Implications for the mantle convection of super-Earths. Geophys J465

Int 204(2):1164–1178, DOI 10.1093/gji/ggv507

466

Kameyama M (2021) Linear analysis on the onset of thermal convection of highly compressible fluids467with variable viscosity and thermal conductivity in spherical geometry: Implications for the mantle

468

convection of super-Earths. Earth, Planets and Space 73:67, DOI 10.1186/s40623-021-01499-w

469

Kameyama M (2022) Numerical experiments on thermal convection of highly compressible fluids with

470

variable viscosity and thermal conductivity in 2-D cylindrical geometry: implications for mantle

471

convection of super-Earths. Geophys J Int 231(2):1457–1469, DOI 10.1093/gji/ggac259

472

Kameyama M, Harada A (2017) Supercontinent cycle and thermochemical structure in the mantle:

473

Inference from two-dimensional numerical simulations of mantle convection. Geosciences 7(4):126,

474

DOI 10.3390/geosciences7040126, 25 pages

475

Kameyama M, Kinoshita Y (2013) On the stability of thermal stratification of highly compressible

476

fluids with depth-dependent physical properties: Implications for the mantle convection of477

super-Earths. Geophys J Int 195(3):1443–1454, DOI 10.1093/gji/ggt321

478

Kameyama M, Yamamoto M (2018) Numerical experiments on thermal convection of highly

479

compressible fluids with variable viscosity and thermal conductivity: Implications for mantle

480

convection of super-Earths. Phys Earth Planet Inter 274:23–36, DOI 10.1016/j.pepi.2017.11.001

481

Kameyama M, Miyagoshi T, Ogawa M (2015) Linear analysis on the onset of thermal convection of482

highly compressible fluids: implications for the mantle convection of super-Earths. Geophys J Int

483

200(2):1066–1077, DOI 10.1093/gji/ggu457

484

485

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
26

486

487

488

489

490

491

492

493

180(1):73–87, DOI 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, DOI 10.1088/0004-637X/700/2/1732
Korenaga J (2010) On the likelihood of plate tectonics on super-earths: Does size matter? Astrophys J
725(1):L43–L46, DOI 10.1088/2041-8205/725/1/L43
Leng W, Zhong S (2008) Viscous heating, adiabatic heating and energetic consistency in compressible
mantle convection. Geophys J Int 173(2):693–702, DOI 10.1111/j.1365-246X.2008.03745.x
Liu X, Zhong S (2013) Analyses of marginal stability, heat transfer and boundary layer properties for494

thermal convection in a compressible fluid with infinite Prandtl number. Geophys J Int

495

194(1):125–144, DOI 10.1093/gji/ggt117

496

497

498

499

500

Miyagoshi T, Tachinami C, Kameyama M, Ogawa M (2014) On the vigor of mantle convection insuper-Earths. Astrophys J 780(L8), DOI 10.1088/2041-8205/780/1/L8
Miyagoshi T, Kameyama M, Ogawa M (2015) Thermal convection and the convective regime diagramin super-Earths. J Geophys Res Planets 120(7):1267–1278, DOI 10.1002/2015JE004793
Miyagoshi T, Kameyama M, Ogawa M (2017) Extremely long transition phase of thermal convection in501

the mantle of massive super-earths. Earth, Planets and Space 69(1):46,

502

DOI 10.1186/s40623-017-0630-6

503

Miyagoshi T, Kameyama M, Ogawa M (2018) Effects of adiabatic compression on thermal convection in

504super-earths of various sizes. Earth, Planets and Space 70:200, DOI 10.1186/s40623-018-0975-5

505

Miyagoshi T, Kameyama M, Ogawa M (2020) Tectonic plates in 3D mantle convection model with

506

507

stress-history-dependent rheology. Earth, Planets and Space 72:70, DOI 10.1186/s40623-020-01195-1
Moresi LN, Solomatov V (1998) Mantle convection with a brittle lithosphere: thoughts on the global

508

tectonic styles of the Earth and Venus. Geophys J Int 133(3):669–682,

509

DOI 10.1046/j.1365-246X.1998.00521.x

510

Noack L, Rivoldini A, Van Hoolst T (2017) Volcanism and outgassing of stagnant-lid planets:

511

Implications for the habitable zone. Phys Earth Planet Inter 269:40–57,

512

DOI 10.1016/j.pepi.2017.05.010

513

Ogawa M (2003) Plate-like regime of a numerically modeled thermal convection in a fluid with

514

temperature-, pressure-, and stress-history-dependent viscosity. J Geophys Res 108(B2):2067,

515

DOI 10.1029/2000JB000069
27

516

517

Papuc AM, Davies GF (2008) The internal activity and thermal evolution of Earth-like planets. Icarus
195(1):447–458, DOI 10.1016/j.icarus.2007.12.016

518

Patankar SV (1980) Numerical Heat Transfer and Fluid Flow. Computational Methods in Mechanics

519

and Thermal Sciences, Hemisphere Publishing Corporation (CRC Press, Taylor & Francis Group),

520

Washington D.C.

521

522

523

524

525

526

527

Ricard Y, Alboussière T (2023) Compressible convection in super-earths. Phys Earth Planet Inter
341:107062, DOI 10.1016/j.pepi.2023.107062
Ricard Y, Alboussière T, Labrosse S, Curbelo J, Dubuffet F (2022) Fully compressible convection forplanetary mantles. Geophys J Int 230(2):932–956, DOI 10.1093/gji/ggac102
Ritterbex S, Harada T, Tsuchiya T (2018) Vacancies in MgO at ultrahigh pressure: About mantle
rheology of super-Earths. Icarus 305:350–357, DOI 10.1016/j.icarus.2017.12.020
Rivera EJ, Lissauer JJ, Butler RP, Marcy GW, Vogt SS, Fischer DA, Brown TM, Laughlin G, Henry

528

GW (2005) A ∼7.5 M⊕ planet orbiting the nearby star, GJ 876. Astrophys J 634(1):625–640,

529

DOI 10.1086/491669

530

531

532

533

534

Schubert G, Turcotte DL, Olson P (2001) Mantle convection in the Earth and Planets. Cambridge
Univ. Press, Cambridge, 940 pp.
Seager S, Kuchner M, Hier-Majumder CA, Militzer B (2007) Mass-radius relationships for solid
exoplanets. Astrophys J 669(2):1279–1297, DOI 10.1086/521346
Shahnas MH, Pysklywec RN (2021) Focused penetrative plumes: A possible consequence of the

535

dissociation transition of post-perovskite at ∼0.9 TPa in massive rocky super-Earths. Geochem

536

Geophys Geosyst 22(8):e2021GC009910, DOI 10.1029/2021GC009910

537

538

Sotin C, Grasset O, Mocquet A (2007) Mass-radius curve for extrasolar Earth-like planets and ocean
planets. Icarus 191(1):337–351, DOI 10.1016/j.icarus.2007.04.006

539

Stamenković V, Breuer D, Spohn T (2011) Thermal and transport properties of mantle rock at high540

pressure: Applications to super-Earths. Icarus 216(2):572–596, DOI 10.1016/j.icarus.2011.09.030

541

Stamenković V, Noack L, Breuer D, Spohn T (2012) The influence of pressure-dependent viscosity on

542

the thermal evolution of super-earths. Astrophys J 748(41):22 pages,

543

DOI 10.1088/0004-637X/748/1/41

544

545

Tachinami C, Senshu H, Ida S (2011) Thermal evolution and lifetime of intrinsic magnetic fields ofsuper-earths in habitable zones. Astrophys J 726(2):70, DOI 10.1088/0004-637X/726/2/70
28

546

547

Tachinami C, Ogawa M, Kameyama M (2014) Thermal convection of compressible fluid in the mantleof super-Earths. Icarus 231:377–384, DOI 10.1016/j.icarus.2013.12.022

548

Tackley PJ (2000) Self-consistent generation of tectonic plates in time-dependent, three-dimensional

549

mantle convection simulations 2. strain weakening and asthenosphere. Geochem Geophys Geosyst

550

1(8):1026, DOI 10.1029/2000GC000043

551

Tackley PJ, Ammann M, Brodholt JP, Dobson DP, Valencia D (2013) Mantle dynamics in552

super-Earths: Post-perovskite rheology and self-regulation of viscosity. Icarus 225(1):50–61,

553

DOI 10.1016/j.icarus.2013.03.013

554

Tosi N, Godolt M, Stracke B, Ruedas T, Grenfell JL, Höning D, Nikolaou A, Plesa AC, Breuer Da,

555

Spohn T (2017) The habitability of a stagnant-lid Earth. Astoronomy and Astrophysics 605:A71,

556

DOI 10.1051/0004-6361/201730728

557

Tsuchida M, Kameyama M (2020) 2-D numerical simulations on formation and descent of stagnant

558

slabs: important roles of trench retreat and its temporal change. Front Earth Sci 8:117,

559

DOI 10.3389/feart.2020.00117

560

Tsuchiya T, Tsuchiya J (2011) Prediction of a hexagonal SiO2 phase affecting stabilities of MgSiO3 and561

CaSiO3 at multimegabar pressures. Proc Nat Acad Sci USA 108(4):1252–1255,

562

DOI 10.1073/pnas.1013594108

563

Turcotte DL, Schubert G (2014) Geodynamics, 3rd edn. Cambridge Univ. Press, 636 pp.

564

Valencia D, O’Connell RJ (2009) Convection scaling and subduction on Earth and super-Earths. Earth

565

566

567

568

569

570

571

572

Planet Sci Lett 286(3-4):492–502, DOI 10.1016/j.epsl.2009.07.015
Valencia D, O’Connell RJ, Sasselov DD (2006) Internal structure of massive terrestrial planets. Icarus
181(2):545–554, DOI 10.1016/j.icarus.2005.11.021
Valencia D, O’Connell RJ, Sasselov DD (2007) Inevitability of plate tectonics on super-earths.
Astrophys J 670(1):L45–L48, DOI 10.1086/524012
Čı́žková H, van den Berg A, Jacobs M (2017) Impact of compressibility on heat transport characteristics
of large terrestrial planets. Phys Earth Planet Inter 268:65–77, DOI 10.1016/j.pepi.2017.04.007
Wagner FW, Sohl F, Hussmann H, Grott M, Rauer H (2011) Interior structure models of solid573

exoplanets using material laws in the infinite pressure limit. Icarus 214(2):366–376,

574

DOI 10.1016/j.icarus.2011.05.027

29