Theor. Comput. Fluid Dyn. (2013) 27:21–40
DOI 10.1007/s00162-011-0250-x

O R I G I NA L A RT I C L E

Masanori Kameyama · Hiroki Ichikawa · Arata Miyauchi

A linear stability analysis on the onset of thermal convectionof a fluid with strongly temperature-dependent viscosity
in a spherical shell

Received: 2 September 2010 / Accepted: 3 December 2011 / Published online: 29 December 2011
© The Author(s) 2011. This article is published with open access at Springerlink.com

Abstract A linear stability analysis was performed in order to study the onset of thermal convection in thepresence of a strong viscosity variation, with a special emphasis on the condition for the stagnant-lid (ST)convection where a convection takes place only in a sublayer beneath a highly viscous lid of cold fluid. We
consider the temporal evolution (growth or decay) of an infinitesimal perturbation superimposed to a Boussinesq fluid with an infinite Prandtl number which is in a static (motionless) and conductive state in a basallyheated planar layer or spherical shell. The viscosity of the fluid is assumed to be exponentially dependent on
temperature. The linearized equations for conservations of mass, momentum, and internal (thermal) energy
are numerically solved for the critical Rayleigh number, Rac , as well as the radial profiles of eigenfunctions
for infinitesimal perturbations. The above calculations are repeatedly carried out by systematically varying (i)the magnitude of the temperature dependence of viscosity, E, and (ii) the ratio of the inner and outer radiiof the spherical shell, γ . A careful analysis of the vertical structure of incipient flows demonstrated that thedominance of the ST convection can be quantitatively identified by the vertical profile of Δh (a measure ofconversion between horizontal and vertical flows), regardless of the model geometries. We also found that,in the spherical shell relevant to the Earth’s mantle (γ = 0.55), the transition into ST convection takes placeat the viscosity contrast across the layer rη  104 . Taken together with the fact that the threshold value ofrη falls in the range of rη for a so-called sluggish-lid convection, our finding suggests that the ST-mode of
convection with horizontally elongated convection cells is likely to arise in the Earth’s mantle solely from thetemperature-dependent viscosity.
Keywords Mantle convection · Temperature-dependent viscosity · Spherical shell · Linear stability analysis ·Stagnant-lid convection

1 Introduction
Temperature-dependent viscosity is one of the most important properties that control the convective flow patterns in the mantle of terrestrial planets [4,26]. For example, the viscosity varies by several to ten orders ofmagnitude associated with the temperature change within the Earth’s mantle [9,25]. Owing to the strongly temperature-dependent viscosity of mantle materials, a highly viscous lid of cold rocks is expected to develop along
the top surface of the mantle, which is quite analogous to the plates or lithospheres of the planetary surfaces[3,9,10,15,16, for example]. From this viewpoint, thermal convection of the fluid with strongly temperaturedependent viscosity has attracted much interest, both experimentally and numerically, from researchers on
mantle convection [14, for a review].
Communicated by: Malik.
M. Kameyama (B) · H. Ichikawa · A. Miyauchi
Geodynamics Research Center, Ehime University, Matsuyama 790-8577, Japan
E-mail: kameyama@sci.ehime-u.ac.jp

22

M. Kameyama et al.

The influence of the temperature-dependent viscosity on the pattern of the thermal convection has been
investigated mainly by theoretical [19] and numerical [3,9,10,15,16] studies in Cartesian domains. From thechanges in vertical flow structures, three convective regimes have been identified depending on the magnitude of temperature dependence of viscosity; the small-viscosity-contrast or whole-layer (WH) regime whereconvection involves all of the fluid in the vessel, the stagnant-lid (ST) regime where convection occurs only
beneath a cold stagnant lid that develops along the surface boundary, and the transitional (TR) regime between
the WH and ST regimes. Convection of WH-, TR-, and ST-modes arises when the temperature dependence ofviscosity is weak, moderate, and strong, respectively, and the transitions between the three convective regimes
are induced by the changes in the stiffness of the cold thermal boundary layer coming from the changes in themagnitude of temperature dependence of viscosity. In particular, a theoretical study by [19] for the onset of
convection in a planar layer estimated that the ST-mode of convection takes place when the viscosity contrast
across the convecting layer exceeds a threshold value of around exp(8)  3,000, which has been regarded as
the critical point for the classification of the convective regimes.
On the other hand, mainly from the numerical models of three-dimensional vessels, it has also beenrecognized that the temperature-dependent viscosity affects the horizontal length scales of convection cells[11,13,20,24,29]. Besides the three regimes described above in terms of the changes in the vertical structures,
another suite of three convective regimes have been separately identified; (a) a regime with narrow convection
cells for weak temperature dependence, (b) a regime sometimes called “sluggish-lid” (SL) mode for moderatetemperature dependence, which is characterized by wide convection cells, and (c) a regime with narrow convection cells for strong temperature dependence. Taken together with the two suites of changes in convective
flow regimes, it has often been conjectured that the convection with strongly temperature-dependent viscosity
is always characterized by narrow convection cells beneath stagnant lids of cold and highly viscous fluid.
Moreover, in almost all of the recent numerical models of convection in three-dimensional spherical shells[18,27,28], the reduction in the horizontal length scales of convection cells has been used as the criterion of
the transition into the ST regimes, even though in some of these studies [18,28], the relations between the twoseries of changes were slightly discussed. Implicit in the above conjecture is the assumption that the two seriesof changes in convective flow regimes are identical. However, as has been demonstrated from a systematic
numerical study for thermal convection with strongly temperature-dependent viscosity in a two-dimensional
rectangular box [3], the changes in horizontal and vertical flow structures belong to different series of transitions for convection. In addition, recent numerical studies by using two-dimensional cylindrical [5,7] andthree-dimensional spherical shell geometries [17] suggested the existence of convection patterns with large
horizontal length scales beneath thick cold thermal boundary layers. It is, therefore, very unlikely that, in athree-dimensional spherical shell geometry, the transition into the ST regime always coincides with the changes
from wide to narrow convection cells. This indicates that in order to precisely classify the convective regimes
in a spherical shell geometry, the transition into the ST-mode of convection must be identified not indirectlyfrom the horizontal but directly from the vertical flow structures.
In this study, we carry out a linear stability analysis on the onset of thermal convection of a fluid withstrongly temperature-dependent viscosity in a spherical shell geometry, in order to classify the convective flow
patterns in terms of the changes in the vertical flow structures and, in particular, to find the condition for thedominance of the ST convection. We will begin by revisiting the analysis by [19] where the transition into
the ST regimes has been identified in a theoretical manner for a convection in a planar layer. By carefully
studying the vertical structures of incipient flows, we will first develop an empirical but quantitative criterion
for the dominance of the ST convection, which is valid regardless of the model geometries. We will then apply
the criterion to the convection in a spherical shell geometry in order to theoretically locate the transition into
the ST regime. By further comparing the conditions for the changes in vertical and horizontal flow structures,
we aim at deepening the insights into convective planforms in spherical shells in the presence of strongly
temperature-dependent viscosity.

2 Model descriptions
In Fig. 1, we schematically show the conceptual models employed here. The onset of thermal convection is considered for a basally heated Boussinesq fluid of infinite Prandtl number with strongly temperature-dependent
viscosity either in a planar layer or in a spherical shell whose temperatures are kept to be Ts and Tb at the
top (inner) and bottom (outer) surfaces, respectively. The aspect ratio of the spherical shell is characterizedby γ ≡ ri /ro , where ri and ro are the inner and outer radii of the shell, respectively. In order to focus on theA linear stability analysis

23

(a)

z

(b)

T=Ts
r

T=Ts
ro

L

T=Tb

0

T=Tb

ri

Fig. 1 Schematic illustration of the numerical model in this study

effects of temperature-dependent viscosity, the spatial variations in other physical quantities (such as thermal
expansivity and thermal conductivity) are ignored.
In this study, we solved for the temporal evolution of an infinitesimal perturbation superimposed to a reference state described by a stationary (motionless) state with steady one-dimensional heat conduction in thevertical (radial) direction. In the following, the quantities with overbars and primes denote those of reference
state and perturbation, respectively. The general forms of the dimensionless equations for the reference state
are written as
1
0 = ∇ 2 T , 0 = v, η = exp −E T −,(1)
2
while those for infinitesimal perturbations are given by
0 = ∇ · v, 
0 = −∇ p  + ∇ · η ∇ ⊗ v  + v  ⊗ ∇ + RaT  eg ,∂T = −v  · ∇T + ∇ 2 T  .
∂t

(2)
(3)
(4)

Here, T is temperature, v is velocity vector, η is viscosity, p is pressure, and eg is the unit vector in vertical(radial) direction positive upward (outward). In deriving these equations, the second-order terms of infinitesimal (primed) quantities are ignored, and the non-dimensionalization is carried out with the length scale ofL (≡ ro − ri for spherical shell) and temperature scale of Tb − Ts . Note also that the value of viscosity η isscaled by ηref , the value at non-dimensional temperature T = 21 . In these equations, the rheological propertiesof the system are expressed by two non-dimensional parameters. One is the Rayleigh number Ra of thermal
convection defined byRa ≡ρref α(Tb − Ts )gL 3,
ηref κ

(5)where ρref , g, α, and κ are the reference values of density, gravity, thermal expansivity, and thermal diffusivity,
respectively. Note that the Rayleigh number Ra is defined with the viscosity ηref . Another is E that describes
the magnitude of temperature dependence in viscosity. At the top and bottom boundaries, we kept the perturbations in temperature, T  , and the normal velocity, v  · eg , to be zero. For the shear velocity at the boundaries,
on the other hand, we assumed either no-slip or free-slip conditions.
2.1 Explicit forms of basic equations
According to a general course of linear analysis [1, for example], the evolutionary equations for perturbations
are solved by separation of variables into temporal, vertical (radial), and horizontal directions. In order tosimplify the present analysis, we assume that the infinitesimal velocity vector, v  , consists only of poloidal
field as:v  = ∇ × ∇ × (W  eg ),(6)where W  is a potential for poloidal field. Note that the velocity field is free from the toroidal components,
since there is no body force acting in the horizontal directions.

24

M. Kameyama et al.

Suppose, for a planar layer, the temporal and spatial dependence of perturbations is



 

Θ(z)
T  (x, y, z, t)
exp i(k x x + k y y) + at .
=W (z)W  (x, y, z, t)(7)Here, a is the growth rate and k x and k y are the wavenumbers of perturbation in the x and y directions,
respectively. By substituting the above equation into (2), (3), and (4), and eliminating pressure p  , we obtain
the equations for the vertical profiles of the perturbations in temperature Θ and poloidal potential of velocity
W as


dT 2
d2
2
aΘ = −
K W + −K + 2 Θ,
dz
dz


2



2
d
dη
d2 dW
d2
d2 η
2
2
2
−K + 2+ 2 K + 2 W,
RaΘ = η −K + 2 W + 2
dz
dz
dz
dz
dz
dz

(8)
(9)

where
K ≡

kx 2 + k y 2(10)is the horizontal wavenumber of perturbations, and the vertical profile of vertical velocity, V , with wavenumber
K is given by V = K 2 W . The boundary conditions for Θ and W are
Θ=W =0d2 W
dW=0= 0 ordz
dz 2

at z = 0, 1

(11)

at z = 0, 1

(12)where z = 1 and z = 0 denote the top and bottom boundaries in non-dimensional units, respectively.
Similarly, for a spherical shell, we will obtain the equations for the radial profiles of the perturbations intemperature Θ and poloidal potential of velocity, W . By expanding the perturbations into spherical harmonics
in the horizontal directions, we obtain the equations for Θ and W of the spherical harmonics degree, , as,
1 dT
2 dΘ
d2 Θ
1
aΘ = −L 2+W+− L 2Θ
2
r dr
dr
r dr
r
 4

2 d2 W
4 dW
d W
L(L − 6)
− 2L 2 + 3L
W+
RaΘ = η
dr 4
r
dr
r
dr
r4


1 d2 W1 dW
dη d3 W
1
−
−
(L
−
1)
W
+2
+
2L
dr dr 3r dr 2r 2 drr3
d2 η d2 W
1
2 dW
+ 2
+ L 2W
−
2
dr
dr
r dr
r

(13)

(14)

where L ≡ ( + 1), and the radial profile of radial velocity V with spherical harmonic degree is given by
V = rL2 W . Note also that there is no dependence on the order of the spherical harmonics in (13) or (14). The
boundary conditions for Θ and W are
Θ=W =0at r = r0 , r1


1 dW
dW
d
= 0 or= 0 at r = r0 , r1
dr
dr r 2 dr

(15)
(16)

where r = r0 and r = r1 denote the inner and outer surfaces of spherical shell in non-dimensional units,
respectively.

A linear stability analysis

25

2.2 Numerical techniques
In this article, we study the variation in the nature of the infinitesimal perturbations that become unstable forthe smallest Rayleigh number by a following strategy. First, for given conditions (model geometry, boundary conditions, and temperature dependence of viscosity), we seek for a critical Rayleigh number, Rac , of
perturbations with wavenumber K or spherical harmonic degree, . In this procedure, we iteratively solved
the eigenequations describing a critical state (i.e., growth rate a = 0) for an eigenvalue (Rac ) together witheigenfunctions (vertical profiles of perturbations of temperature Θ and vertical velocity V ). The details of
the numerical technique can be found in Appendix. Among the series of eigenvalues and eigenfunctions inhand, we then search the one that yields the absolute minimum of the critical Rayleigh number, Rac0 , and thewavenumber, K c0 , or the spherical harmonic degree, c0 , corresponding to Rac0 .
The differential equations for Θ and W (Eqs. 8 and 9 for a planar layer, 13 and 14 for a spherical shell) are
discretized with a sixth-order accuracy in space using equally spaced 257 grid points. As will be demonstrated
in later sections, the reliability of the present numerical technique is verified by the comparison with the earlier
results of [12,19].

3 Results
3.1 Control experiments: onset of convection in a planar layer
We first carried out control experiments similar to the work by [19]: the onset of convection in a planar layer.
In the following, we denote no-slip boundaries by R (rigid) and stress-free boundaries by F (free-slip) for
convenience. We then use a fractional notation with the top of the fraction being the condition at the top
boundary and the bottom of the fraction being the condition at the bottom boundary. Thus, F/R means that the
top surface is free-slip and the bottom surface is rigid. We consider four cases: F/F, R/R, F/R, and R/F.
In Fig. 2, we show the plots of (a) the absolute minimum value of critical Rayleigh number, Rac0 , and(b) the wavenumber, K c0 , corresponding to Rac0 , against E for various boundary conditions indicated in the
figure. From these experiments, we confirmed the earlier finding [19] that three regimes can be distinguished:(i) the regime with weakly temperature-dependent viscosity (E  1.5) where Rac0 and K c0 are nearly constant, (ii) that with moderate dependence (1.5  E  8) where Rac0 increases and K c0 is nearly constant ordecreases moderately depending on the boundary conditions, and (iii) that with strong dependence (E  8)
where Rac0 reaches a peak and then decreases, while K c0 rises rapidly. In particular, for the cases with free-slip
top boundaries (F/F and F/R), the plots of Rac0 have maxima at around E = 8, which has been commonly
regarded as the transition into the “Stagnant-Lid” (ST) regime. Indeed, as can be seen from the figure, the
(b)
6

2500

5

2000

R/R

1500

R/F

critical wavenumber

critical Rayleigh number

(a)
3000

F/R

1000

4

R/R
3

R/F
F/R
2

F/F

F/F
500

1
0

2

4

6

8

E=ln(ηmax/ηmin)

10

12

14

0

2

4

6

8

10

12

14

E=ln(ηmax/ηmin)Fig. 2 The plots of a the absolute minimum value of critical Rayleigh number, Rac0 , and b the wavenumber, K c0 , corresponding
to Rac0 , against E for various boundary conditions indicated in the figure (e.g., “R/F” for a rigid top and a free-slip bottom
boundaries). Note that the Rayleigh number is defined with the viscosity for non-dimensional temperature T = 0.5

26

M. Kameyama et al.

(a) F/ F(b) R/ F
1.0

1.0

E=0

E=0
0.8

normalized heightnormalized height

0.8

0.6

E=5
0.4

E=9
E=11

0.2

E=5
0.6

0.4

E=9
0.2

E=7

E=7
0.0

E=11

0.0

0

1

2

3

4

0

1

normalized vertical velocity

2

3

4

normalized vertical velocity

Fig. 3 Plots of vertical profiles of eigenfunctions for vertical velocity V for a F/F and b R/F boundaries with various values of
E indicated in the figure. In each plot, the values of V are normalized in order to set their volumetric average to unity

values of Rac0 are the largest at E = 8.3 for both cases. Figure 2 also shows that for sufficiently large E, theplots for the case with free-slip top boundaries (F/F and F/R) become close to those for the corresponding casewith rigid top boundaries (R/F and R/R). This is consistent with the nature of ST mode of convection, where
the motion of the highly viscous fluid is negligibly small along the top cold surface boundary.
The transition into ST regime with increasing E can also be observed from Fig. 3, where we compared
the vertical profiles of eigenfunctions for vertical velocity V for the critical states with various values of E
obtained with (a) F/F and (b) R/F boundary conditions. The plots of V in the figure show that, as E becomeslarger, the values of V become smaller and the positions of their local maxima become lower. In particular,
for sufficiently large E, the profiles of eigenfunctions for vertical velocity, V , become almost identical bothfor the F/F and R/F cases (see, for example, the profiles with E = 11), implying that the cold fluid behaves asa stagnant lid regardless of the boundary conditions along the top surface.
The fact that the transition into ST regime occurs at around E  8 has been analytically estimated by theearlier work by [19], based on an assumption that convection begins in a sublayer spanning 0 ≤ z ≤ ẑ(< 1).
ˆ given by,
They consider the “local” Rayleigh number of the sublayer Ra3ˆ ≡ ρref α[(Tb − Ts )ẑ]g(L ẑ) = Ra ηref ẑ 4 ,Raη̂ref κ
η̂ref

(17)

where (Tb − Ts )ẑ is the temperature drop in the sublayer. By further taking the reference viscosity η̂ref in theˆ to Ra can be written as
sublayer to be that at its mid-depth (T = 1 − 2ẑ at z = 2ẑ ), the ratio of Ra

 
ˆ
1
Ra
ẑ
= ẑ 4 exp E
−
.
Ra
2 2

(18)

Differentiating the above equation with respect to ẑ gives
∂
∂ ẑ

ˆ
Ra
Ra

=

ˆ
Ra
Ra



4
E
−
ẑ
2


,

(19)

ˆ
implying that for E ≥ 8, there exists a sublayer with ẑ = 8/E ≤ 1, which yields a maximum of Ra/Ra
largerthan 1. This allows the onset of convection only in a sublayer that does not extend to the upper boundary.
This estimate of threshold E agrees well with our numerical finding that the Rac0 has a maximum at aroundE = 8.3.
As has already been demonstrated in the earlier work [19], the validity of the above analysis can be also
ˆ c0 defined
confirmed from the variations of the absolute minimum value of critical “local” Rayleigh number RaA linear stability analysis

27

critical Rayleigh number

3000

2500

R/R
2000

F/R
R/F

1500

F/F

1000

500
0

10

20

30

E=ln(hmax/hmin)

ˆ c0 for the stagnant-lid (ST)Fig. 4 The plots against E of the absolute minimum value of critical “local” Rayleigh number Ra
regime (E ≥ 8; thick lines). The meanings of the labels for the boundary conditions are the same as in Fig. 2. For comparison,
we also plotted against E the values of Ra c0 by thin lines

by


 
1 ẑ c04
ˆ
Ra c0 ≡ Ra c0 ẑ c0 exp E
−
,
2
2

(20)

ˆ c0 is the “local” Rayleigh number, corresponding to Rac0 , of the basal sublayer
where ẑ c0 = 8/E. Namely, Raˆwhose Ra is maximum. We show in Fig. 4 the plots against E the values of Ra c0 over a broad range of E ≤ 30
(or, the viscosity contrast of exp(30)  1013 of the entire layer). As can be seen from the plots, the values ofˆ c0 are always larger than those of Ra c0 for the ST convection (see also 18). In addition, according to theRaˆ c0 are expected toassumption that the onset of convection occurs only in a basal sublayer, the values of Ra
be constant in the stagnant-lid regime. The plots in the figure clearly show that the above conjecture is valid
ˆ c0 are around 2,038 and 1,568 when the lower
for the range of E studied here. The asymptotic values of Ra
boundary is rigid and free-slip, respectively.
In order to quantitatively study the changes in the vertical flow structures toward the dominance of the ST
convection, we consider a quantity Δh defined by
∂v y
∂vx
∂vz
dV
+
=−
=−
exp[i(k x x + k y y) + at]
∂x
∂y
∂z
dz
≡ Δh (z) exp[i(k x x + k y y) + at].

(21)

By definition, the quantity Δh (z) is a measure of the divergence of horizontal velocity (vx and v y ) in thehorizontal plane at a non-dimensional height z. Note that the value of Δh (z) is also a measure of the extent
of conversion between horizontal and vertical flows at a height z. This indicates that the source regions ofascending and descending flows can be estimated from its vertical profiles. In Fig. 5, we show the verticalprofiles of Δh (z) for various values of E obtained for the case with F/F boundary condition. The vertical profile
of Δh for an isoviscous case (E = 0) is symmetric with respect to the mid-depth of the layer (z = 0.5), withtwo maxima at top and bottom boundaries (z = 0 and z = 1). In addition, the minimum at z = 0.5 comesfrom the fact that the eigenfunction of vertical velocity V has a maximum at the depth.
As can be seen from the changes in the vertical profiles of Δh with increasing E in Fig. 5, the transition in
the convective regimes can be identified by the difference in the positions and magnitudes of the local maxima
of Δh occurring in the upper and colder part of the layer (i.e., the maxima other than that at z = 0). The verticalprofiles of Δh (z) for E  3 have maxima at z = 1. This is consistent with the feature of so-called whole-layer(WH) mode of convection [3,15] occurring with a weakly temperature-dependent viscosity, where the coldest
fluid along the top surface is involved in an overall convective motion. For E = 5, in contrast, the profile of Δhhas a local maximum in the mid-depth at around z = 0.6. This indicates that the motion of the highly viscous

28

M. Kameyama et al.

1.0

E=0
E=3
E=5
E=7
E=8
E=8.3
E=9

normalized height

0.8

0.6

0.4

0.2

0.0
10

-2

10

-1

10

0

10

1

10

2

|Δh(z)/Δh(z=1)|

Fig. 5 Plots of the vertical profiles of Δh , a measure of conversion between horizontal/vertical flows, for various values of E
indicated in the figure. Shown are the results of the cases with F/F boundaries. Note that plotted are the absolute values of theratios of Δh to the values at the top surface (z = 1)

fluid in the uppermost part of the layer is slightly decoupled from that of less viscous fluid in the interior. As
E increases further, the local maximum occurs at a deeper part of the layer, and the value of local maximum
Δhmax becomes larger. When the ST convection becomes dominant (E = 8.3; see also Fig. 2), the value ofΔhmax becomes more than 10 times larger than that of Δh (z = 1). This indicates that the fluid motion becomes
negligible in the uppermost part of the layer, and, in other words, a stagnant lid of cold and highly viscous
fluid develops along the top surface boundary.
In summary, the present numerical method has successfully reproduced the change in the style of convection
depending on the changes in temperature-dependent viscosity. The transition can be identified by the quantity
Δh , which is a measure of the extent of conversion between horizontal and vertical flow at a particular height.
In addition, we found that, when the stagnant-lid (ST) convection becomes dominant, the vertical profiles ofΔh (z) have local maxima in the mid-depth whose magnitude is more than 10 times larger than that at the top
surface (z = 1). From this finding, we will identify in the next section, the transition in convective patterns inthe spherical shell geometry based on the vertical profiles of Δh .
3.2 Onset of convection in spherical shell geometries
In this section, we will study the onset of convection in a spherical shell of a fluid with temperature-dependent
viscosity.
First we verified the validity of the present technique, by comparing our numerical results with those ofan earlier study by [12]. Here we calculated the values of the critical Rayleigh number, Rac , in a sphericalshell with γ = 0.55, by varying the spherical harmonic degrees, , up to 6 and the viscosity contrast due tothe temperature dependence rη ≡ exp(E) up to 103 . In Table 1, we summarize the comparison of the present
results with those by [12]. As can be seen in Table 1, the values of Rac obtained from two different studies
show an excellent agreement for the ranges of and rη considered here. From this comparison, we conclude
that our technique successfully captures the critical state of convection in the spherical shell geometry.
In the rest of this section, we will study the condition for the stagnant-lid (ST) convection for the sphericalshell geometry by varying the magnitude of the temperature-dependent viscosity, E. In particular, we will
focus on how the transition is affected by the model geometry, including the ratio of the inner and outer radiiof the spherical shell, γ .
3.2.1 Empirical estimates for the condition for the ST-mode for spherical shell geometry
To clearly show the differences from the results of the control experiments, we will start with a comparisonof the cases with thin (aspect ratio γ = 0.95) and thick (γ = 0.55) spherical shells. Note that a model with

A linear stability analysis

29

Table 1 The values of the critical Rayleigh number of the fluid with a temperature-dependent viscosity in a spherical shell of
γ = 0.55 whose inner and outer surfaces are free-slip boundariesrη = exp(E)Spherical harmonic degree
1
2

100

This study[12]

1,864.34
1,864.34

888.84888.84711.95711.95729.39729.38851.31851.31

1,063.95
1,063.95

101

This study[12]2,136.64
2,136.65

1,085.05
1,085.05

913.73913.73963.06
963.07

1,135.81
1,135.81

1,418.09
1,418.09

102

This study[12]2,313.50
2,313.49

1,456.04
1,456.04

1,387.90
1,387.91

1,535.19
1,535.20

1,811.00
1,811.01

2,199.59
2,199.60

3

4

5

6

103

This study
2,664.01
2,196.40
2,248.43
2,457.30
2,762.52
3,147.64[12]
2,663.98
2,196.39
2,248.43
2,457.29
2,762.52
3,147.63
In the table, we compare the values obtained in this study and those obtained by [12], for various values of spherical harmonic
degrees of the infinitesimal perturbations and the viscosity contrast due to the temperature dependence rη ≡ exp(E)551
γ = 0.95 is characterized by a volumetric average with a reference temperature of Tav = 1,141
 0.4829
and inner and outer radii of r0 = 19 and r1 = 20, respectively. On the other hand, a shell with γ = 0.55 has
77
r0 = 11/9, r1 = 20/9, and Tav = 247
 0.3117. In Fig. 6, we show the plots of (a) the absolute minimum
value of critical Rayleigh number, Rac0 , and (b) the degree of spherical harmonics, c0 , of infinitesimal perturbation corresponding to Rac0 , against E obtained both for the cases with γ = 0.95 and γ = 0.55 for various
boundary conditions indicated in the figure.
A qualitative comparison between the results presented in Figs. 2 and 6 shows that the changes in theconvective regimes for both γ are similar with those obtained for a planar layer. Indeed, from the changes in
Rac0 and c0 , three regimes can be distinguished depending on E: (i) the regime with small E where Rac0 andc0 are nearly constant, (ii) that with moderate E where Rac0 increases and c0 is nearly constant or decreases
moderately depending on the boundary conditions, and (iii) that with large E where Rac0 reaches a peak andthen decreases, while c0 rises rapidly. In addition, for sufficiently large E, the plots for the case with free-slip
top boundaries (F/F and F/R) become close to those for the corresponding case with rigid top boundaries (R/Fand R/R). This indicates that the ST mode of convection emerges in spherical shell geometries for γ = 0.95
and γ = 0.55 as well as in a planar layer.
However, a detailed comparison of the variations in Rac0 presented in Fig. 6a illuminates a quantitative
difference in the changes in the convective regimes depending on the aspect ratio of the spherical shell γ .
For the case with a thin spherical shell (γ = 0.95), the variations in Rac0 are very similar to those obtainedin a planar layer throughout the range of E studied here (see also Fig. 2). In particular, for the cases with
free-slip top boundaries (F/F and F/R), the plots of Rac0 have maxima at around E  8. For the case with athick spherical shell (γ = 0.55), in contrast, the variations in Rac0 are significantly different from those for aplanar layer. For a given E, the values of Rac0 are larger than those obtained for a planar layer regardless ofthe boundary conditions. In addition, the maxima in the plots of Rac0 against E occur at larger values of E
than those for a planar layer and for a spherical shell with γ = 0.95 (Note that, as will be discussed later indetail, the maxima in the plots of Rac0 do not necessarily imply the dominance of the ST convection).
We can also easily notice from Fig. 6b that the aspect ratio of the spherical shell γ affects the smoothness ofthe changes in convective regimes. The changes in Rac0 occur very smoothly and continuously with increasing
E for γ = 0.95, while they occur in a quite discontinuous manner for γ = 0.55. This comes from the differences in the thickness of the spherical shell and, in other words, the horizontal extent of the convecting layercompared to the vertical (radial) one. Indeed, the values of c0 are large ( c0  30) for γ = 0.95, while they
are small ( c0  10) for γ = 0.55. Since the horizontal length scale of perturbation is roughly proportional to
the reciprocal of spherical harmonic degree , the effect of the change in is minor for sufficiently large .
In order to precisely locate the dominance of the ST convection for spherical shell geometry, we consider
the radial profile of Δh , a measure of the extent of conversion between horizontal and vertical flows. Now inthe spherical geometry, the quantity Δh is calculated from the eigenfunction for radial velocity V (or that for
its poloidal potential W ) by

1
∂
1 ∂vφ
1 ∂  2 
1 d  2  m(sin θ vθ ) +=− 2
r vr = − 2
r V Yr sin θ ∂θr sin θ ∂φ
r ∂r
r dr
≡ Δh (r )Y m ,

(22)

30

M. Kameyama et al.

(a)

(b)
110

3000

γ=0.95100

critical harmonic degree

critical Rayleigh number

γ=0.95
2500

R/R

2000

1500

R/F
F/R

1000

F/F

90
80
70

R/R
60

R/F
F/R

50
40

F/F
500

30
0

2

4

6

8

10

12

14

0

2

4

E=ln(ηmax/ηmin)
5000

8

10

12

14

10

12

14

10

γ=0.55γ=0.55

R/R

4000

critical harmonic degree

critical Rayleigh number

6

E=ln(ηmax/ηmin)

F/R

3000

F/F

R/F

2000

R/R
5

R/F

F/F

F/R

2

4

1000

0

0
0

2

4

6

8

E=ln(ηmax/ηmin)

10

12

14

0

6

8

E=ln(ηmax/ηmin)Fig. 6 The plots of a the absolute minimum value of critical Rayleigh number Rac0 and b the degree of spherical harmonics c0of infinitesimal perturbation corresponding to Rac0 against E for various boundary conditions indicated in the figure (e.g., “R/F”
for a rigid top and a free-slip bottom boundaries). Shown are the cases with thin (γ = 0.95; top) and thick (γ = 0.55; bottom)
spherical shells. Note again that the Rayleigh number is defined with the viscosity for non-dimensional temperature T = 0.5,
which is of little significance to the dominance of the ST convection in a spherical shell geometry (see Sect. 3.2.2 for details)

where Y m is the spherical harmonics with degree and order m. Note the difference in the definition of Δhin (22) from that for a planar layer (see also 21). We show in Fig. 7 the radial profiles of Δh (r ) for various
values of E and the boundary condition of F/F obtained for γ = 0.95 and γ = 0.55, while Fig. 8 are those ofV (r ). Also shown in the figures are the values of c0 of the perturbations for the corresponding cases.
The comparison of the radial profiles of Δh (r ) in Fig. 7 clearly shows that the condition for the dominanceof the ST-mode of convection is affected by the aspect ratio of the spherical shell γ . For γ = 0.95, the transition
occurs in a quite similar manner to the case with a planar layer (see also Fig. 5): at E = 8.3 the local maximumof Δh (r ) occurs at around r − r0  0.3, and its value Δhmax becomes more than ten times larger than thevalue of Δh (r ) at the top surface (r = r1 ). For γ = 0.55, in contrast, the ratio of Δhmax to Δh (r = r1 ) do not
exceed the threshold value of ∼ 10 unless E becomes as large as around 9.2, implying that the dominance of
the ST convection takes place at around E = 9.2.
We also note that the transition is accompanied by the changes in c0 , which significantly affect the vertical profiles of Δh across the transition for small γ . This is mainly because the values of c0 are smaller for
γ = 0.55 than for γ = 0.95. As can be seen from Fig. 8, the change in c0 causes a significant change in theradial profiles of V when c0 is small.

A linear stability analysis

31

(a)

(b)
1.0

γ=0.95E=0 (l=43)
E=5 (l=36)
E=7 (l=41)
E=8 (l=53)
E=8.3 (l=57)
E=9 (l=64)

0.8

0.6

normalized height r-r0

normalized height r-r0

1.0

0.4

0.2

γ=0.55E=0 (l=3)
E=3 (l=3)
E=8 (l=2)
E=9.1 (l=3)
E=9.2 (l=4)
E=10 (l=5)

0.8

0.6

0.4

0.2

0.0

0.0
-2

10

10

-1

0

10

10

1

10

2

10

-2

10

|Δh(r)/Δh(r=r1)|

-1

0

10

10

1

2

10

|Δh(r)/Δh(r=r1)|

Fig. 7 Plots of the radial profiles of Δh for various values of E indicated in the figure for a thin (γ = 0.95) and b thick (γ = 0.55)
spherical shells. Also indicated in the figure are the values of c0 of the perturbations. Shown are the results of the cases with
F/F boundaries. Note that plotted are the absolute values of the ratios of Δh to the values at the outer surface (r = r1 ) against thenon-dimensional height from the inner surface (r = r0 )

(a)

(b)
2.2

γ=0.95E=0 (l=43)
E=5 (l=36)
E=7 (l=41)
E=8 (l=53)
E=8.3 (l=57)
E=9 (l=64)

0.8

0.6

γ=0.55E=3 (l=3)
E=7 (l=2)
E=8 (l=2)
E=9.1 (l=3)
E=9.2 (l=4)
E=10 (l=5)

2.0

normalized radiusnormalized height r-r0

1.0

0.4

0.2

1.8

1.6

1.4

0.0
0

1

2

3

0

4

1

2

3

4

5

normalized vertical velocity

normalized vertical velocity

Fig. 8 Plots of vertical profiles of eigenfunctions for vertical velocity V for spherical shells with a γ = 0.95 and b γ = 0.55
with various values of E indicated in the figure. Also indicated in the figure are the values of c0 of the perturbations. Shown are
the results of the cases with F/F boundaries plotted against the non-dimensional height from the inner surface (r = r0 ). In eachplot, the values of V are normalized in order to set their volumetric averages to unity
Table 2 The values of E andof the spherical shellc0 at the dominance of the ST-mode of convection obtained for various values of the aspect ratio

γ

0.95

0.9

0.85

0.8

0.75

0.7

0.65

0.6

0.55

0.5

0.45

0.4

E

8.3
57

8.3
27

8.3
17

8.4
12

8.4
9

8.5
7

8.8
6

8.9
4

9.2
4

9.3
3

9.8
2

10.5
2

c0

We carried out similar calculations by varying the aspect ratio of the spherical shell γ in the range of
0.4 ≤ γ ≤ 0.95 in order to study how the conditions for the dominance of the ST convection are affectedby the variation in γ . We show in Table 2 the values of E and c0 for the transition with various values of γ .
Table 2 clearly shows that, for the range of γ studied here, the dominance of the ST-mode of convection occurs
at a larger E for smaller γ . We also note that, for smaller γ , the values of c0 for the perturbations become
smaller at the dominance of the ST-mode. This is because the horizontal extent of the spherical shell becomes
small compared to the vertical (radial) one for smaller γ .

32

M. Kameyama et al.

3.2.2 Analytical estimates for the condition for the transition into ST-mode for spherical shell geometry
In Sect. 3.2.1, we have numerically estimated the condition for the dominance of the ST-mode of convection
in a spherical shell geometry, based on the changes in radial structures of the incipient flows. We found thatthe transition occurs in a quite different manner from that in a planar layer. In this section, we will analytically
explore the condition for the transition into the ST regime in a spherical shell. In particular, we will demonstratethat a slight modification of the idea by [19], which has been used for a planar layer in Sect. 3.1, successfully
reproduces the values of E for the transition and their dependence on the aspect ratio of the spherical shell γ .
We consider a ST-mode of convection in a spherical shell whose inner and outer radii are r0 and r1 in
non-dimensional units, respectively. Suppose that, following the idea by [19], the convecting instability grows
only in a basal sublayer spanning r0 ≤ r ≤ r2 (where r0 < r2 ≤ r1 ) and T2 ≤ T ≤ 1 (where T2 is the valueof T at r = r2 ). In addition, we consider two “modified” Rayleigh numbers Ra∗ and Ra2∗ ; the former is forthe entire layer and the latter is local in the lower sublayer. Note that one needs to estimate Ra∗ and Ra2∗
based on a definition other than that employed so far, since the dominance of the ST convection can hardly beobserved from the values of Ra used earlier (see also Fig. 6). In other words, by appropriately choosing Ra∗
and Ra2∗ , the regime for ST-mode of convection can be characterized by the existence of r2 , which satisfies
Ra∗ < Ra2∗ .
In Fig. 9, we schematically show the key assumption employed here for an appropriate choice of the Rayleigh number for a fluid with temperature-dependent viscosity in a spherical shell geometry. The “modified”
Rayleigh number Ra∗ for the entire layer is defined with the value of viscosity of temperature T = T∗ at
r = r∗ given by ρref α(Tb − Ts )gL 3
1= Ra exp E T∗ −
Ra∗ ≡
,



2
ηref exp −E T∗ − 21 κ√r∗ ≡ r0 r1 ,

 √
γ −γ
1
r0 r1
1
T∗ ≡
−
,
=
r1 − r0 r∗ r1
1−γ

(23)
(24)
(25)

respectively. As can be seen from Fig. 9a, the value of r∗ is chosen to be the radius in the spherical shell where
the radial temperature gradient is equal to an average temperature gradient across the entire layer. Indeed, by
using
r0 r1
T =
r1 − r0



1
1
−
r
r1



r0 r1 1dT=−
dr
r1 − r0 r 2

,

(a)

(26)

(b)

r1

r1

heightheight

r2

r*

r2*

T*

r0
0

1

temperature

T2

r0

T2*

0

1

temperature

Fig. 9 a Schematic illustration of the meanings of “representative” radius r∗ and temperature T∗ used for a definition of Rayleighnumber for a convection of a fluid with temperature-dependent viscosity in a spherical shell geometry. b The same as a but forthe convection in a lower sublayer relevant to the ST-mode of convection

A linear stability analysis

33

(a)

(b)
0.5

"representative" temperature T*"representative" height r*-r0

0.5

0.4

0.3

Tav

0.4

0.3

0.2
0.2

0.4

0.6

0.8

aspect ratio of spherical shell γ

1.0

0.2

0.4

0.6

0.8

1.0

aspect ratio of spherical shell γ

Fig. 10 The plots against the aspect ratio of the spherical shell γ of the a position r∗ and b temperature T∗ which yield a“representative” viscosity. Also plotted in b is the values of volumetric average of temperature Tav of the reference state

we obtain ddTr = −1 at r = r∗ . We show in Fig. 10, the variations of r∗ and T∗ depending on the aspectratio of the spherical shell γ . As can be seen from the figure, both of the values of r∗ and T∗ decrease as γ
decreases from 1. In particular, the value of T∗ significantly deviates from 0.5 with decreasing γ , implying agreater difference between Ra and Ra∗ for smaller γ . Similarly, as indicated in Fig. 9b, let us define the local
Rayleigh number Ra2∗ for the basal sublayer with the value of viscosity of T = T2∗ at r = r2∗ given by√r2∗ ≡ r0 r2 .
(27)
We can also note that r2∗ satisfies
r0 r1 1
1 − T2dT=−.
(r = r2∗ ) = −dr
r1 − r0 r2∗ 2
r2 − r0
Here, in rewriting from the second to the third terms in the right-hand side, we used


r0 r1
1
1
r0 r1 r2 − r0
1 − T2 = 1 − T (r = r2 ) = 1 −−.
=r1 − r0 r2 r1
r1 − r0 r2∗ 2
In the following, we will use γ2 ≡ r0 /r2 (note that γ ≤ γ2 < 1) and the expression of T2∗ given by



 √
γ2 − γ
1
1
r0 r1
1
1
r0 r1
−
−
.
=
=
T2∗ =
√
r1 − r0 r2∗ r1
r1 − r0
r0 r2 r1
1−γ

(28)

(29)

Let us now consider the ratio of the Rayleigh number in the basal sublayer to that of the entire layer given by
r2 − r0 3 exp(−E T∗ )Ra2∗
= (1 − T2 )Ra∗
r1 − r0exp(−E T2∗ )3
4(1 − γ2 )γ
exp [E(T2∗ − T∗ )].
(30)=
4(1 − γ )γ2 3Differentiating (30) by γ2 gives



∂
Ra2∗
γ3
(1 − γ2 )4∂exp [E(T2∗ − T∗ )]=∂γ2 Ra∗(1 − γ )4 ∂γ2
γ2 3


∂ T2∗
3 + γ2
Ra2∗
+E
−
=
Ra∗
γ2 (1 − γ2 )
∂γ2


1
1
3 + γ2
Ra2∗
+E
−
.
(31)
=
√
Ra∗
γ2 (1 − γ2 )2(1 − γ ) γ2

34

M. Kameyama et al.

1.2

γ=0.55
1.0

E=9
E=9.5
E=10
E=10.5
E=11

Ra2*/Ra*

0.8

0.6

0.4

0.2

0.0
0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00

aspect ratio of sublayer γ2Fig. 11 The plots of the ratio of the local Rayleigh number of the sublayer Ra2∗ to Ra∗ against the aspect ratio of the sublayer
γ2 obtained for a spherical shell of γ = 0.55 for various values of E indicated in the figure

threshold values of E

10.5

10.0

9.5

9.0

8.5

8.0
0.2

0.4

0.6

0.8

1.0

aspect ratio of spherical shell γ

Fig. 12 The variations against the aspect ratio of the spherical shell γ of the values of E c for the transition into the stagnant-lid
regimes given by (33). Also shown by thick dots are the values of E c by an empirical estimate described in Sect. 3.2.1

This means that the value of Ra2∗ /Ra∗ is maximum when1
3 + γ21=E√ .
γ2 (1 − γ2 )2(1 − γ ) γ2

(32)

In Fig. 11, we plotted the values of Ra2∗ /Ra∗ given by (30) against the aspect ratio of the sublayer γ2 forvarious values of E. As can be seen from the figure, for sufficiently large E, the plots of Ra2∗ /Ra∗ against γ2
have local maxima within the range of γ ≤ γ2 < 1 and, in other words, there exists a sublayer of aspect ratio
γ2 where Ra2∗ becomes larger than Ra∗ . By further assuming that the transition into the ST regime occurs at
E ≡ E c which gives γ2 = γ , we obtain the relations between γ and E c asEc =2(3 + γ )
.
√
γ

(33)

In order to verify our analytical estimate, we plotted in Fig. 12, the values of E c for the transition into the
ST-mode of convection for various values of the aspect ratio of spherical shell γ given by (33), together withthose derived empirically in Sect. 3.2.1. From the comparison of E c , we obtained a satisfactory agreement
between the empirical and analytical estimates.

A linear stability analysis

35

(b)

(a)

3000

γ=0.95
F/F
F/R
R/F
R/R

2500

modified critical Rayleigh number

modified critical Rayleigh number

3000

2000

1500

1000

γ=0.55
F/F
F/R
R/F
R/R

2500

2000

1500

1000

500

500
0

10

20

30

0

10

E=ln(ηmax/ηmin)20

30

E=ln(ηmax/ηmin)Fig. 13 The plots against E of the absolute minimum values of critical “modified” Rayleigh numbers Ra∗c0 and Ra2∗c0 obtained
for the cases with a thin (γ = 0.95) and b thick (γ = 0.55) spherical shells. Plotted by thick lines are Ra∗c0 of the entire layer,
while by thin lines are Ra2∗c0 local to the sublayer for the ST convection. The meanings of the labels for the boundary conditions
are the same as in Fig. 2

We show in Fig. 13, the plots against E for the cases with γ = 0.95 and γ = 0.55 of the absolute minimumof critical “modified” Rayleigh number Ra∗c0 calculated by (23), together with those of Ra2∗c0 for the sublayerof γ2 = γ2c0 which satisfies (32) for the ranges of E in the ST regime. The plots of Ra∗c0 clearly show thatthe values for both γ have maxima at the values of E which correspond to the transition into the ST regime(E  8.3 for γ = 0.95, E  9.2 for γ = 0.55). In addition, the plots of Ra2∗c0 show that in the ST regime,
the values of Ra2∗c0 are always larger than those of Ra∗c0 . This is another evidence for the validity of ouranalytical estimate regardless of the aspect ratio γ of the spherical shell. We also note that, on the other hand,
(a)

(b)

8

γ=0.55

7

No Convection

6

log10rh

5
4

(c)

3
2
1
0
2

3

4

5

6

7

8

log10Rab
Fig. 14 a A regime diagram of convective flow patterns in the spherical shell with γ = 0.55 in the plane of Rab (Rayleigh number
defined with the viscosity for T = Tb ) versus rη ≡ exp(E). Circles, triangles, and squares indicate that convective patterns
with narrow, wide, and again narrow cells have been reported in earlier studies with increasing rη , respectively; black symbols
from [28], blue symbols from [18], green symbols from [31], and red symbols from [2]. The thick solid and broken lines show
the locations of the threshold for the onset convection obtained in this study for the planar and spherical models, respectively.
The red crosses in the figure indicate the locations of the regime boundaries between the ST and other modes on the thick lines.
b and c Examples of convective flow patterns obtained by [2] for b Rab = 5 × 107 and rη = 106 and c Rab = 5 × 107 andrη = 104 . Shown are the three-dimensional plots of the distributions of lateral thermal anomaly δT ≡ T − T , where T is thehorizontally averaged temperature at given r . Indicated in cyan are the cold thermal anomalies with δT ≤ −0.1, while in yellow
are the hot anomalies with δT ≥ 0.1. Red spheres show the inner bottom boundary of the convecting layer

36

M. Kameyama et al.

there is a slight difference between the cases with γ = 0.95 and 0.55 in the variations of Ra2∗c0 for verylarge E. For a thin spherical shell (γ = 0.95), the values Ra 2∗c0 are almost constant for E up to 30, as in thecase within a planar layer (see Sect. 3.1). For γ = 0.55, in contrast, the values of Ra 2∗c0 moderately varydepending on E with a gradual increase for large E, since the aspect ratio γ2c0 of the sublayer significantly
increases with increasing E for a thick spherical shell.
4 Discussion and concluding remarks
We performed a linear stability analysis for the onset of thermal convection of a fluid with strongly temperature-dependent viscosity in a basally heated planar layer or a spherical shell. In this study, we focused on
the condition for the stagnant-lid (ST) mode of convection with increasing the temperature dependence ofviscosity, with particular emphasis on the changes in vertical structures of incipient flows. We first developed
an empirical but quantitative criterion for the dominance of the ST convection, based on the well-establishedresults for a planar layer. The criterion comes from an extension of the idea of [19], which postulates that in
the ST regime, the onset of convection occurs only in a sublayer that does not extend to the upper boundary.
We proposed a quantity Δh , defined by (21) or (22) depending on the model geometries, which represents theextent of conversion between horizontal and vertical flows at a particular height and, in other words, the source
regions of ascending and descending flows. From a series of calculations varying the model geometries as wellas the values of parameter E describing the temperature-dependent viscosity, we found that the dominanceof the ST convection is significantly affected by the model geometries, such as the aspect ratio of sphericalshell γ . The dominance occurs at E  8 in a thin spherical shell (large γ ), as in the case with a planar layer,
while it occurs at larger E in a thick spherical shell (small γ ). Our findings indicate the unsoundness of theclassification of the convecting regimes in a thick spherical shell relevant to the Earth’s mantle (γ  0.55),
which has been simply based on the results obtained in a planar layer.
We also developed an analytical estimate for the transition into the ST regimes in a spherical shell, following the idea of [19] where convection is assumed to begin in a sublayer with a maximum Rayleigh number. Our
analysis showed that the fact that the values of E for the dominance of the ST regimes depend on the aspectratio γ can be explained by a slight modification of the definition of the Rayleigh number from that employed
for a planar layer. As has been demonstrated in Sect. 3.2.2, it is quite important to appropriately choose thevalues of “representative” viscosity and, in other words, the “representative” temperature T∗ or the position r∗of the fluid of strongly temperature-dependent viscosity in a spherical shell. This is in a stark contrast to thecase in a planar layer, where the values of T∗ and r∗ are naturally chosen to be those at the mid-depth of thelayer (T∗ = 21 and r∗ = 21 ). Indeed, as can be seen in Fig. 10a, the values of r∗ significantly deviate from thatfor the mid-depth of the spherical shell (= 21 (r0 + r1 )) with decreasing γ . In addition, for γ < 1, the valuesof T∗ are different from either the mean of the top and bottom values (= 0.5) or the volumetric average Tav inthe entire layer. This indicates that an appropriate choice of the Rayleigh number is of crucial importance forunderstanding the convecting regimes in a thick spherical shell.
We will apply our results to a classification of the convecting regimes in a thick spherical shell relevantto the Earth’s mantle. We show with the thick solid line in Fig. 14a the variation of the absolute minimum of
the critical Rayleigh number against the temperature dependence of viscosity obtained for a spherical shell of
γ = 0.55. In the figure, we plotted, instead of E and Ra, the viscosity contrast across the layer rη ≡ exp(E)and the Rayleigh number defined with the viscosity at the bottom surface Rab ≡ exp(E/2)Ra, for a directcomparison with the earlier studies [3,9,10,15]. We also show by the dashed line the relation with Rab andrη for a planar layer. In addition, the red crosses indicate the transitions into the ST regimes on the curves. As
can be seen from the figure, the curves of rη against Rab have slight bends near the red crosses, regardless ofthe model geometries. This is also an evidence for the significance of the criterion for the transition into the
ST regimes proposed in this study.
We also show in Fig. 14a, the classification of convective regimes, inferred from the changes in the horizontal length scales of convection cells, obtained by earlier numerical experiments of thermal convection inthree-dimensional spherical shells [2,13,18,28,31]. Among the symbols in the figure, the triangles indicate
the conditions for a so-called “sluggish-lid” (SL) mode of convection which is characterized by convection
cells of large horizontal scales (see Fig. 14c), while the squares indicate those for a convection with small
horizontal scales (see Fig. 14b). In particular, the mode of convection denoted by the squares in the figure has
been considered to be identical to the ST mode in most of the earlier studies [13,20,24,28], although it has beenreported from two-dimensional convection models with a moderately wide box [3] that the changes from wide
to narrow convection cells with increasing rη do not necessarily coincide with that in vertical flow structures.

A linear stability analysis

37

The results from our linear stability analysis suggest that the conclusion drawn from two-dimensional
Cartesian models [3] is most likely to be true for the model with spherical shell geometry relevant to the
Earth’s mantle (γ = 0.55). Indeed, as can be seen in Fig. 14a, the changes in the vertical structures of incipient flows occur at around rη  104 (i.e., E  9.2), while, from earlier numerical studies [2,28,31], thatthe collapse of the SL-mode of convection occurs at the values of rη larger than 104 . It can, therefore, be
conjectured that, for a range of 104  rη  104.5 , there exists an ST-mode of convection with convection cellsof large horizontal length scales. Of course, the present findings are not entirely applicable to the classificationof convecting regimes from finite-amplitude flows, since our results are valid only near the critical condition(indicated by the thick lines in Fig. 14a). In addition, the transition into the ST regimes does not necessarily
occur at rη  104 when the Rayleigh number is sufficiently higher than the critical value. However, as can beseen from Fig. 14c, the SL-mode of convection obtained with rη = 104 and the Rayleigh number far larger thanthe critical value (Rab = 5 × 107 ) have a feature quite similar to that of “elongated ST-mode” of convection,
which has been identified by an earlier numerical model using two-dimensional moderately wide Cartesiangeometry [3]. The convective flow pattern of these cases is characterized by convection cells of large aspectratio superimposed by several minor ascending and descending plumes due to the instability in the top andbottom thermal boundary layers (TBLs). Moreover, the minor descending plumes originating from the baseof the cold TBL do not break the wide convection cells into narrow ones, implying that the cold TBL is stiff
enough to prevent the instability from penetrating upward into the TBL. The similarity of “SL-mode” and“elongated ST-mode” suggest that our conjecture is most likely to be valid for, at least, around rη = 104 evenwhen the Rayleigh number is sufficiently large.
One of the most important inferences from this study is, we believe, that much more careful studies areneeded for the transitions in convective flow patterns in a spherical shell geometry. This is largely because it isof crucial importance for understanding the convective planforms in the mantle of the terrestrial planets. It has
been commonly believed that the convective motion in the mantle of terrestrial planets is characterized by (i)thick stiff lids along the top cold surfaces and (ii) convection cells of large horizontal scales [6,14,17,25,30, forexample]. However, it has also be considered that a mode of thermal convection can be hardly obtained solely
from the effect of temperature-dependent viscosity, which bears these two features simultaneously. For instance,
when the viscosity is strongly temperature dependent, convection cells are thought to have small horizontal
scales beneath highly viscous and stiff lids. On the other hand, for a convection with cells of large horizontal
scales, the cold TBLs are not viscous enough to act as a stiff lid. In order to simultaneously reproduce these two
features, various agents, other than temperature-dependent viscosity, have been incorporated into the numericalmodels of mantle convection, such as the effects of depth-dependent viscosity [28], internal heating [8], adiabatic compression [23], and yielding within the cold TBLs [21,22,27]. Contrary to earlier conjectures describedabove, the present results indicate that the ST-mode of convection with elongated convection cells is most likely
to occur, under certain conditions, solely from the temperature-dependent viscosity. Certainly, such agents, if
they operate together with the temperature-dependent viscosity, would play an important role in simultaneously
reproducing convective patterns with wide convection cells beneath cold and stiff lids. It may be an oversimplified view, however, that such agents are inevitable mechanisms in mantle convection of terrestrial planets, since
the effects of temperature-dependent viscosity have not been fully resolved on the changes in vertical flow structures. We, thus, conclude that it is still an important issue to carefully study the effects of temperature-dependent
viscosity on the convective planforms, both vertically and horizontally and with a three-dimensional sphericalshell geometry, in order to deepen the insights into the nature of mantle convection of terrestrial planets.
Acknowledgments We thank Drs. Takashi Nakagawa, Shijie Zhong, and David A. Yuen for fruitful discussion. We also thankSabrina Whitaker for correcting the grammatical and typological errors in the manuscript. We thank three anonymous reviewers
for valuable comments which greatly improved the manuscript. We greatly acknowledge thorough support from the Global COE
program from the Ministry of Education, Culture, Sports, and Technology (MEXT) of Japan. MK acknowledges financial supports
from Grant-in-Aid for Scientific Research (C; #19540452) and (B; #22340127) from Japan Society for the Promotion of Science.
Open Access This article is distributed under the terms of the Creative Commons Attribution Noncommercial License which
permits any noncommercial use, distribution, and reproduction in any medium, provided the original author(s) and source arecredited.

A The details of numerical techniques
In this study, the critical Rayleigh number Rac for a given condition (i.e., boundary conditions, the values of
E and K ) is numerically calculated in two steps. First, for a given Rayleigh number Ra, we calculated the38

M. Kameyama et al.

infinitesimal perturbation which yields the largest growth rate amax . We then sought for the value of Rac which
gives amax = 0, by iteratively varying Ra. In this section, we will describe the details of the each technique inturn. In the following, for simplicity, we will write the eigenequation to be solved as
Ax = λx.

(34)

Here the matrix A represents the spatial discretization, and the eigenvalue λ and eigenvector x stand for thegrowth rate and infinitesimal perturbations of temperature and velocity, respectively.
A.1 Solution of the largest eigenvalue λmax
The numerical difficulty in solution of the largest eigenvalue λmax largely comes from the fact that λmax canbe either positive or negative depending on the conditions to be considered. This implies that basic numerical
techniques (such as the power method) which seek for the dominant eigenvalue (i.e., the one with the largest
absolute value) cannot be directly used. However, as will be shown below, a slight modification of the original
eigenequations enables us to use the power method in order to obtain λmax .
In the present numerical algorithm, we seek for the eigenvalues λ̃ of an “auxiliary” matrix Ã given byÃ = I + Δt A,(35)

instead of directly solving for the eigenvalues λ of A in (34). Here I is an identity matrix and Δt is a small
positive constant similar to a time increment. By multiplying an eigenvector x of A with (35), we getÃx = x + Δt Ax = (1 + λΔt)x.

(36)

This means that both A and Ã have the same set of eigenvectors and that the following relation holds between
the eigenvalues λ and λ̃ of the two matrices:λ̃ = 1 + λΔt.
In addition, let us assume that the parameter Δt is taken to be small enough to satisfy1,Δt <|λmin |
where λmin is the smallest eigenvalue of A. Together with (37), this leads to
λmin
= 0.
λ̃ ≥ 1 + λmin Δt > 1 +|λmin |

(37)

(38)

(39)

Namely, all of the eigenvalues λ̃ of Ã are positive. This means that the largest eigenvalue λ̃max can be easily
solved for by, for example, the power method since λ̃max = |λ̃max |. In other words, the eigenvalue λmax of thematrix A can be derived from λ̃max of the matrix Ã via
λmax =

λ̃max − 1
,
Δt

(40)

regardless of the sign of λmax .
We finally consider an appropriate choice of the positive parameter Δt. Suppose a following “auxiliary”
evolutionary equation
dy
= Ay.
(41)
dt
By discretizing the above equation in the direction of t by a first-order explicit scheme, we get
yn+1 − yn= Ayn .
Δt
By using (35), this equation can be rewritten as
yn+1 = (I + Δt A) yn = Ã yn .

(42)

(43)

This means that the multiplication of Ã is just a time-marching calculation of the (41) by a first-order explicit
scheme. In other words, the parameter Δt must be chosen to satisfy the Courant-Friedrichs-Lewy (CFL)
condition of (41). In the present study, the values of Δt are determined from the CFL condition of (4).

A linear stability analysis

39

A.2 Solution of the critical Rayleigh number
By using the technique in the previous subsection, we can obtain the value of amax for a given Rayleigh number
Ra. Next, we seek for the value of the critical Rayleigh number Rac , which yields amax = 0, by the well-known
bisection method. The bisection method is a root-finding method that repeatedly bisects an interval and thenselects a subinterval where a root must lie for further processing. Owing to the monotonous relation between
Ra and amax , the bisection method always converges to the correct Rac .
References
1. Chandrasekhar, S.: Hydrodynamic and Hydromagnetic Stability. The Clarendon Press, Oxford (1961)
2. Kameyama, M., Kageyama, A., Sato, T.: Multigrid-based simulation code for mantle convection in spherical shell using
Yin-Yang grid. Phys. Earth Planet. Inter. 171(1–4), 19–32 (2008). doi:10.1016/j.pepi.2008.06.025
3. Kameyama, M., Ogawa, M.: Transitions in thermal convection with strongly temperature-dependent viscosity in a widebox. Earth Planet. Sci. Lett. 180(3-4), 355–367 (2000). doi:10.1016/S0012-821X(00)00171-0
4. Karato, S.I., Wu, P.: Rheology of the upper mantle: a synthesis. Science 260(5109), 771–778 (1993). doi:10.1126/science.
260.5109.771
5. Kiefer, W.S., Li, Q.: Mantle convection controls the observed lateral variations in lithospheric thickness on present-day
Mars. Geophys. Res. Lett. 36, L18203 (2009). doi:10.1029/2009GL039827
6. Konopliv, A.S., Banerdt, W.B., Sjogren, W.L.: Venus gravity: 180th degree and order model. Icarus 139(1), 3–18 (1999).
doi:10.1006/icar.1999.6086
7. Li, Q., Kiefer, W.S.: Mantle convection and magma production on present-day Mars: Effects of temperature-dependentrheology. Geophys. Res. Lett. 34, L16203 (2007). doi:10.1029/2007GL030544
8. McNamara, A.K., Zhong, S.: Degree-one mantle convection: dependence on internal heating and temperature-dependentrheology. Geophys. Res. Lett. 32, L01301 (2005). doi:10.1029/2004GL021082
9. Moresi, L.N., Solomatov, V.: Numerical investigation of 2D convection with extremely large viscosity variations. Phys.
Fluids 7, 2154–2162 (1995)
10. Ogawa, M., Schubert, G., Zebib, A.: Numerical simulations of three-dimensional thermal convection in a fluid with strongly
temperature-dependent viscosity. J. Fluid Mech. 233, 299–328 (1991)
11. Ratcliff, J.T., Schubert, G., Zebib, A.: Effects of temperature-dependent viscosity on thermal convection in a sphericalshell. Physica D Nonlinear Phenomena 97(1–3), 242–252 (1996). doi:10.1016/0167-2789(96)00150-9
12. Ratcliff, J.T., Schubert, G., Zebib, A.: Steady tetrahedral and cubic patterns of spherical shell convection with temperaturedependent viscosity. J. Geophys. Res. 101(B11), 25,473–25,484 (1996). doi:10.1029/96JB02097
13. Ratcliff, J.T., Tackley, P.J., Schubert, G., Zebib, A.: Transitions in thermal convection with strongly variable viscosity. Phys.
Earth Planet. Inter. 102(3–4), 201–212 (1997). doi:10.1016/S0031-9201(97)00013-7
14. Schubert, G., Turcotte, D.L., Olson, P.: Mantle Convection in the Earth and Planets. Cambridge University Press, Cambridge (2001)
15. Solomatov, V.S.: Scaling of temperature- and stress-dependent viscosity convection. Phys. Fluids 7, 266–274 (1995)
16. Solomatov, V.S., Moresi, L.N.: Three regimes of mantle convection with non-Newtonian viscosity and stagnant-lid convection on the terrestrial planets. Geophys. Res. Lett. 24(15), 1907–1910 (1997). doi:10.1029/97GL01682
17. Šrámek, O., Zhong, S.: Long-wavelength stagnant-lid convection with hemispheric variation in lithospheric thickness: link
between Martian crustal dichotomy and Tharsis?. J. Geophys. Res. 115, E09010 (2010). doi:10.1029/2010JE003597
18. Stemmer, K., Harder, H., Hansen, U.: A new method to simulate convection with strongly temperature- and pressure-dependent viscosity in a spherical shell: applications to the Earth’s mantle. Phys. Earth Planet. Inter. 157(3–4), 223–249 (2006).
doi:10.1016/j.pepi.2006.04.007
19. Stengel, K.C., Oliver, D.S., Booker, J.R.: Onset of convection in a variable-viscosity fluid. J. Fluid Mech. 120, 411–
431 (1982). doi:10.1017/S0022112082002821
20. Tackley, P.J.: Effects of strongly variable viscosity on three-dimensional compressible convection in planetary mantles. J.
Geophys. Res. 101(B2), 3311–3332 (1996)
21. Tackley, P.J.: Self-consistent generation of tectonic plates in time-dependent, three-dimensional mantle convection simulations 1. Pseudoplastic yielding. Geochem. Geophys. Geosyst. 1(8), 1021 (2000). doi:10.1029/2000GC000036
22. Tackley, P.J.: Self-consistent generation of tectonic plates in time-dependent, three-dimensional mantle convection simulations 2. Strain weakening and asthenosphere. Geochem. Geophys. Geosyst. 1(8), 1026 (2000). doi:10.1029/2000GC000043
23. Tackley, P.J., Xie, S., Nakagawa, T., Hernlund, J.W.: Experimental studies of mantle convection: philosophy, accomplishments and thermo-chemical structure and evolution. In: van der Hilst, R.D., Bass, J., Matas, J., Trampert, J. (eds.) Earth’s
Deep Mantle: Structure, Composition, and Evolution, Geophysical Monograph, vol. 160, pp. 83–99. AGU (2005)
24. Trompert, R., Hansen, U.: On the Rayleigh number dependence of convection with a strongly temperature-dependent viscosity. Phys. Fluids 10(2), 351–364 (1998). doi:10.1063/1.869527
25. Turcotte, D.L., Schubert, G.: Geodynamics. pp. 456 2nd edn. Cambridge University Press, Cambridge (2002)
26. Weertman, J.: The creep strength of the Earth’s mantle. Rev. Geophys. 8, 145–168 (1970)
27. Yoshida, M.: Mantle convection with longest-wavelength thermal heterogeneity in a 3-D spherical model: degree one or
two?. Geophys. Res. Lett. 35, L23302 (2008). doi:10.1029/2008GL036059
28. Yoshida, M., Kageyama, A.: Low-degree mantle convection with strongly temperature- and depth-dependent viscosity in athree-dimensional spherical shell. J. Geophys. Res. 111, B03412 (2006). doi:10.1029/2005JB003905
29. Zebib, A.: Linear and weakly nonlinear variable viscosity convection in spherical shells. Theor. Comput. Fluid
Dyn. 4(5), 241–253 (1993). doi:10.1007/BF00417930

40

M. Kameyama et al.

30. Zhao, D.: Global tomographic images of mantle plumes and subducting slabs: insight into deep Earth dynamics. Phys. Earth
Planet. Inter. 146(1-2), 3–34 (2004). doi:10.1016/j.pepi.2003.07.032
31. Zhong, S., McNamara, A., Tan, E., Moresi, L., Gurnis, M.: A benchmark study on mantle convection in a 3-D sphericalshell using CitcomS. Geochem. Geophys. Geosyst. 9, Q10017 (2008). doi:10.1029/2008GC002048